Simulation of capillary infiltration into packing structures by the Lattice-Boltzmann method for the optimization of ceramic materials

In this work we want to simulate with the Lattice-Boltzmann method in 2D the capillary infiltration into porous structures obtained from the packing of particles. The experimental problem motivating our work is the densification of carbon preforms by reactive melt infiltration. The aim is to determine optimization principles for the manufacturing of high-performance ceramics. Simulations are performed for packings with varying structural properties. Our analysis suggests that the observed slow infiltrations can be ascribed to interface dynamics. Pinning represents the primary factor retarding fluid penetration. The mechanism responsible for this phenomenon is analyzed in detail. When surface growth is allowed, it is found that the phenomenon of pinning becomes stronger. Systems trying to reproduce typical experimental conditions are also investigated. It turns out that the standard for accurate simulations is challenging. The primary obstacle to overcome for enhanced accuracy seems to be the over-occurrence of pinning.


INTRODUCTION
Porous media are ubiquitous in science and engineering. Typical examples ripple across a number of disciplines such as hydrology, oil recovery, materials science, printing and solar technology, to cite but a few (Bear, 1972;Dullien, 1992;Furler et al., 2012;Goodall & Mortensen, 2013;Koponen et al., 1998). In simulations, their structures are investigated according to a variety of properties (mechanical, thermal and fluid dynamics) mainly through finite elements methods (FEM) (Alawadhi, 2010;Bohn & Garboczi, 2003). Capillarity is the spontaneous infiltration of a fluid into a porous medium (Alava, Dube, & Rost, 2004;Washburn, 1921). The adhesive forces between the liquid and solid phases are responsible for this phenomenon. Indicatively, its characteristic length scale is the micron (de Gennes, Brochard-Wyart, & Quéré, 2004). Our interest is in the role of the porous structure for capillary infiltration. Simulations are carried out by using the Lattice Boltzmann (LB) method (Benzi, Succi, & Vergassola, 1992;Chen & Doolen, 1998;Succi, 2009;Sukop & Thorne, 2010;Wolf-Gladrow, 2005). This is a numerical scheme for approximating the hydrodynamic behavior, especially suited for problems involving complex boundaries and interface phenomena. Its main advantage over other approaches resides in the discretization of the velocity space and a statistical treatment of particle motion and collisions. Applications encompass disparate critical systems (Bao et  The actual problem motivating our research is the reactive infiltration of molten silicon (Si) into carbon (C) preforms ( ). This is a complex phenomenon associated with a wetting transition and surface growth. A feature of this process is that Si reacts with C to form silicon carbide (SiC). Si wets SiC but does not wet C (Bougiouri et al., 2006). Furthermore, the reaction-formed phase causes the thickening of the surface behind the invading front (Bougiouri et al., 2006;Einset, 1996Einset, , 1998Israel et al., 2010;Messner & Chiang, 1990). The interaction with the fluid flow results in a retardation of the infiltration up to its interruption because of pore clogging. This research gains increasing attention also in industrial practice. Liquid Si infiltration (LSI) is a common manufacturing route to high-performance ceramics. Prominent advanced applications include armor and brake systems (Fan et al., 2008;Krenkel & Berndt, 2005;Roberson & Hazell, 2003). Significant engineering efforts concentrate on the control of the effects of reactivity. Indeed, SiC formation provides additional hardness to the resulting ceramic material, reducing also the amount of residual unreacted Si detrimental for high-temperature functioning. On the other hand, pore closure can stop infiltration before full densification is attained. The microstructure of the C preform is of course of importance for optimal impregnation, besides for the operating behavior of the final ceramic product in specific applications (e.g., the length of reinforced C fibers) ( In previous works we have studied the effects of surface growth for systems composed of single capillaries both uniform and structured in 2D (Sergi, Grossi, Leidi, & Ortona, 2014, 2015. The focus of the present investigation is on making more precise the pore configurations of particular relevance for the industrial application of LSI. The strategy consists in first realizing porous structures approximating real arXiv:1601.07697v2 [cond-mat.mtrl-sci] 27 Jul 2016 porous preforms. To this end, we generate packing systems by using particles differing in basic attributes (size, shape and distribution). In principle, these parameters can be controlled by the selection and processing of powders. The packing structures are assembled by means of the random sequential addition (Sherwood, 1997;Widom, 1966). This is a simple and effective algorithm allowing to reach reasonable volume fractions for random packings. The simulations for capillary infiltration aim at highlighting the function of specific pore and structural characteristics. The systems of this work display enhanced randomness with pores and channels of different size and orientation inducing complex dynamics for the invading front. The significance of this investigation resides in the possibility to translate the findings into inputs for preform preparation in order to improve the quality of ceramic materials.
Our simulations indicate that the interface dynamics can be assumed to be responsible for small effective radii. Pinning of the contact line is the primary process affecting capillary infiltration. Our work allows to give evidence for the mechanism determining pinning in random systems. Packings of particles based on bimodal size distributions turn out to provide features reducing the strong retardation of liquid invasion associated with pinning. Larger particles introduce walls creating channels that are filled faster while smaller particles make the pore structure more uniform. As a result, the interface travels longer distances before the loss of its curvature. Furthermore, the probability to encounter the surface of other particles is higher. Our data also show that the inertia effects induced by accelerations remain of secondary importance. With surface growth, it is found that pore closure can be effective and pinning turns out to be enhanced. The retardation effects are stronger for fast infiltrations (higher degree of orientation disorder), that remain nevertheless more advantageous. Also with surface growth, bimodal size distributions lead to better results. Simulations with higher resolution for the porous structure are also performed for a direct comparison with experimental expectations. In this case, the effect of surface growth on pore closure is mainly indirect because pinning is stronger. Furthermore, the relative strength of capillary forces is too weak. This means that the details of the infiltration process suffers from discrepancies with real conditions. For example, this implies that the ratio of the infiltrated distance to the thickness of the growing surface (i.e., the width of the reaction-formed phase) is too small in simulations. Progress seem possible especially by reducing the phenomenon of pinning.

MODELS AND ANALYSIS
For the formalism of the LB models the interested reader is addressed to Sergi et al. (2014). In the sequel, the same notation and terminology are employed. We recall that in experiments with molten Si and some alloys the invading front exhibits a kinetics linear with time (Israel et al., 2010;Voytovych et al., 2008). The standard Washburn law for capillary infiltration in vacuum predicts a parabolic time depen-dence. The observed experimental behavior can be reproduced with the simulations by assuming that the systems are composed by two fluid components (Chibbaro, 2008 ). This means that the effect of the reaction at the contact line is equivalent to the resistance due to the presence of a second non-wetting component as dense and viscous as the wetting one (Sergi et al., 2014). This approach does not explain the origin of the linear behavior for the infiltration kinetics generally ascribed to effects of the reactivity (Israel et al., 2010;Voytovych et al., 2008).
The differential equation describing capillary infiltration for two fluids with the same density ρ and dynamic viscosity µ into a single channel in 2D reads The first term on the left-hand side accounts for inertial forces while the second one introduces viscous forces. The term on the right-hand side is due to capillary forces. z designates the position of the invading front, r is the radius of the channel and L its length. When inertial forces are neglected, for a shrinking radius as r(t) = r 0 − kt (r 0 is the initial radius and k is the reaction-rate constant) a solution to the above equation is given by z 0 is the initial position. For completeness, in the absence of reaction, a solution to Eq. 1 taking into account also inertial forces is obtained by ) where V cap = γ/µ and t d = ρr 2 /3µ. For more general porous systems, the radius of the channel is replaced with the effective radius. Relevant studies for reactive infiltration based on dynamics can be found in Asthana v indicates a characteristic velocity for the system. Q 1 compares the reaction with surface tension, Q 2 the reaction with viscosity, Q 3 the reaction with inertia. For completeness, Q 4 = ρrv 2 /γ compares inertia with capillary forces. In our terminology, the first force becomes more important as the dimensionless number increases. In our notation, the Reynolds and capillary numbers are given by Re = ρvr/µ and Ca = vµ/γ. In the onset of pore closure, Q 1 tends toward zero faster and thus surface tension is the dominant force. Instead, inertia has the weakest effects since Q 3 diverges while Q 1 and Q 2 tend to 0. By using for the time-varying effective radius the average value r 0 /2, we arrive at the same conclusions.
For the characterization of the porous systems we use Darcy law (Sukop & Thorne, 2010): v is the interstitial velocity, ϕ the porosity, P the pressure and K the permeability. Another quantity that we take into account is the tortuosity (Bear, 1972). It is defined as the ratio of the actual length of the capillary paths to the corresponding length of the sample. For the computations we use the formula < u > / < u x > (Duda, Koza, & Matyka, 2011; Matyka & Koza, 2012), where u is the magnitude of the fluid velocity of components u x and u y . It is worth recalling that at the interface there are spurious currents (Wagner, 2003).

PACKING SYSTEMS
It is standard practice to reproduce the microstructure of composite materials with random packings (Scocchi et al., 2013; Torquato, 2002). The problem consists in filling the empty space with particles arranged in a disordered way. The realization of high volume fractions can also be quite challenging (Donev et al., 2004). The principle of the algorithm for the random sequential addition is to place one particle after the other in the simulation domain (Sherwood, 1997;Widom, 1966). In case of overlapping between a test particle and an existing one, the trial is not accepted and a new particle is considered. This procedure is iterated up to the desired filling fraction. With this algorithm it is not easy to exceed values beyond 0.6. In order to generate the packing systems for the LB simulations we employ the algorithm proposed in the article by  Table III: Data for the packing structures obtained with fibers and rhombs. For the side and the volume, the first value refers to the sticks and the second to the rhombs. For fibers, the side of the square base is 0.06 long. The height of the rhombs is h = 0.2, 0.16, 0.12, 0.08, that is, similar to the length of the four equal sides. The reported filling fractions are associated with different orientations of the particles: aligned, misaligned, misaligned and tilted, respectively. Figure 1 shows examples of porous configurations in 2D.
Type Type  Table V: Rhombs composing the packing systems infiltrated with inert boundaries (see Tab. II). For every quantity, multiple values corresponding to different orientations of fillers: aligned, misaligned, misaligned and tilted (see Fig. 1).  Table VI: Fibers added to the packing systems infiltrated without reactive boundaries (see Tab. III). For every quantity, multiple values refer to different orientations of fillers: aligned, misaligned, misaligned and tilted (see Fig. 1). Sergi, D'Angelo, Scocchi, and Ortona (2012). The method consists in decomposing the particles into small cubes. The lengths are measured in units of the small cubes composing the particles according to the rule = (n+1)r. n is the number of cubes and 2r their side length. Here, we consider spheres, rhombs and sticks with a certain degree of randomness for their size. For a given value of n, the associated lengths are chosen with a Gauss distribution of average and standard deviation σ = /4. When the rotation of the particles is allowed, the rotation angles around the three axes are at most θ x = θ y = π/2 and θ z = π/4. The domain is a cube of side 1 divided into N = 256 bins. The porous medium for the LB simulations is obtained by suitably replicating the packing systems. The minimal distance between two small cubes is d min = 2/N. An external pressure is applied toward the center of the domain in order to obtain more compact structures when necessary, as it is the case for higher degrees of disorder, for example. The other parameters for the pressure are the same as in Sergi et al. (2012).
The packing systems aim at reproducing basic characteristics resulting from the use of common commercial powders for preform preparation. In order to make the distinction between round-shaped and sharp-edged morphologies for the    Tabs. I and II. particles, we consider spheres and rhombs. For example, the latter geometry can approximate the elongated shape of SiC particles. Fibers are also used as an attempt to understand the role of long walls. Concerning the size, the distribution is bimodal, multimodal or monomodal (coarse and fine particles). The properties of the packing systems are summarized in Tabs. I-III. Figure 1 shows the section for some systems.

SETTINGS FOR LB SIMULATIONS
As said before, the LB method is based on the discretization of the velocity space and a statistical treatment of the particle motion. In the BGK approximation for the collisions (Bhatnagar, Gross, & Krook, 1954), LB simulations remain versatile, efficient and accurate. Proper hydrodynamic behavior can be recovered in the incompressible limit at low Mach numbers (Chen & Doolen, 1998). The domain for simula-    tions is N x = 1500 lu long and N y = 255 lu wide. The porous structure has a length of L = N x /2 and the solid-liquid interface starts at x = 2N x /5. The total number of timesteps is given by N t = 2 · 10 6 ts. In the region occupied by the porous structure the boundaries of the simulation domain are solid and subject to the common bounce-back rule; elsewhere they are periodic. This is done in order to reduce the phenomenon of pinning near narrow-to-wide parts (   The process of surface growth from the reaction is introduced by relying on a purely heuristic treatment. A description based on solute transport as in Sergi et al. (2014Sergi et al. ( , 2015 exhibits limitations related to the interface for solute transport. Namely, inside the porous structures the interface can form, leading to the diffusion of solute in the whole system. Therefore, in this work the reaction is implemented by making grow the surface of the solid phase behind the contact line at regular intervals of time. Precisely, the first layer around the solid boundaries is converted to solid phase after every N t /10 timesteps. The reaction-rate constant is thus given by k = 5·10 −6 lu/ts. This framework has the advantage to be suitable for a critical volume of simulations allowing to observe effects of pore structure and surface growth on fluid flow. Unless stated otherwise, these settings are maintained in the sequel.

Packing systems
In Tab. IV we report the results of simulations for systems obtained from the packing of spheres in the absence of surface reaction. Here z indicates the penetration attained on average at the end of the process. ∆z is an estimate of the deviations from the average. This quantity is determined from the maximum and the minimum of the front displacement. It is to be noted that the error ∆z is comparable to the infiltrated distance. This means that the fluid advances into the porous medium through selected pathways (see Fig. 2). It can be expected that these fingers then contribute to fill the pores left behind. It is worth noticing that the average error is 82 lu while the size of the largest particles is on average of 51 lu. In Fig. 3 is shown the kinetics of invasion. It turns out that the infiltration is slower for a monomodal distribution of the particle size. By fitting the data with Eq. 3, we derive the infiltration velocity v inf,2 and the effective radius r eff (see Tab. IV). It can be seen that the effective radius is one to two orders of magnitude smaller than the average particle size, as stressed in previous works (Einset, 1996;Patro, Bhattacharyya, & Jayram, 2007). The permeability K is determined under the conditions usually employed to reproduce Poiseuille flow, i.e. with a single fluid, periodic boundaries and an external acceleration (Sukop & Thorne, 2010). In Tab. IV it can be seen that the permeability for spheres is higher with bimodal size distributions. In order to assess the role of interface dynamics, the tortuosity is computed in two different ways. λ 1 is determined in the presence of an external acceleration as for the permeability. λ 2 is instead obtained from the process of infiltration driven by capillary forces. The results of Tab. IV show that λ 2 is two orders of magnitude higher. As a consequence, the interface dynamics induced by the pore structure plays a prominent role. v inf,1 is obtained by an average over the velocities inside the porous media in the direction of infiltration, i.e. x. The difference with v inf,2 is more marked when the process is slower. For a proper basis of comparison with systems in-cluding the reaction, the characteristic dimensionless numbers are calculated by means of v inf,1 (see Tab. IV). It is also interesting to consider the Mach number defined as Ma = v/c s , where c s = 1/3 lu/ts is the speed of sound. From Tab. IV it can be seen that the order of magnitude of Ma is 10 −4 . It follows that in this study the numerical error associated with the LB scheme remains in the tolerance limits (He & Luo, 1997;Sukop & Thorne, 2010).
The results of simulations without reactivity for packings of particles with the shape of rhombs are reported in Tab. V. Figure 3 shows the kinetics in the case of aligned rhombs. Comparison with the outcome for packings of spheres indicates that the structure arising from rhombs is more subject to pinning. It turns out that, as the degree of disorder increases, the infiltration becomes faster. The fastest infiltrations are obtained with tilt and misalignment, even though the tortuosity λ 2 attains a larger value on average. The dynamics with spheres remain still faster or at least comparable. A possible explanation for this observation is that the arrangement of rhombs can give rise to channels that are filled easily. The shortcoming is that the phenomenon of pinning becomes more frequent. In 3D in a real preform, we expect that the advantages should prevail since the connectivity of the pore chambers is higher and the number of pathways increases. For example, the particle walls above or below a 2D section and almost parallel to it contribute to reduce pinning. This point will become clearer from the discussion on the advantages of bimodal size distributions. As expected, the permeability K is higher with increasing orientation disorder (see Tab. V). For aligned systems it also turns out that K attains higher values in the presence of larger particles. Figure 4 shows the range of variation of Ca. The same behavior holds for Re: the same relation between Re and Ca is of course an identity. More marked variations are observed for Q 4 . The increase can amount to one order of magnitude. It is interesting to note that the deviations from the average values are stronger for rhombs since the filling of distinct channels and depinning introduce more accelerations. From the order of magnitude of the average values of Ca, Q 4 and Re, we conclude that the effects of inertia are weak and capillary forces remain the dominant ones. Furthermore, in Fig. 5 it can be seen that for a given capillary number, the tortuosity seems to be slightly lower for aligned rhombs in a comparison with spheres. This tendency is recognized also in the presence of orientation disorder and for fluctuations of Ca. For clarity, this suggests that, when the fluid accelerates, the tortuosity is not higher and thus the resistance of inertia has a weak effect. So, in this case the averages of dimensionless numbers still provide an exhaustive description. As expected, it follows that the interface dynamics induced by the pore structure seems to play a more important role for the pronounced difference between λ 1 and λ 2 . Spurious currents at the interface should average out since we could verify that, when the fits allowing to obtain v inf,2 are more precise, the discrepancies with v inf,1 are smaller. Again, for aligned rhombs, it arises that liquid penetration tends to be better with a bimodal distribution. This trend is less clear with misalignment and tilt. It should be said that, in this case, 2D sections loose to a large extent the specific characteristics of the particles composing the packings. For this reason, it might be presumed that, in 3D, bimodal distributions should remain more advantageous. Examination of the data associated with the maxima and minima for infiltration reveals that pinning is the discriminating factor coming into play. The basic mechanism at the origin of pinning is that the meniscus advances slowly while loosing its curvature, resulting in a drop of capillary pressure. In the presence of channel walls, the front can travel a longer distance since the pinning sites are not vertically aligned. As a consequence, the meniscus can encounter the surface of new particles with higher probability before halting (see Fig. 2). In the packings with small particles the pinning sites are instead generally separated by small distances. The probability that the invading front stops moving is thus higher. It is also interesting to note that rhombs with bimodal distributions pack more efficiently. On average the filling fraction of their 2D sections is 14% higher than in the case of monomodal distributions. However, their infiltrations remain faster on average of 3%. It can be seen that the results are more sensitive to specific porous configurations than to the porosity. For more clarity, test simulations are performed for different porosities for the three bimodal systems composed of rhombs with particles aligned and misaligned. For filling fractions of 36% − 50% the penetration depth varies within a narrow range on average. This is partly in line with our discussion for the conditions helping the advancement of fluid. The porosity seems to start becoming detrimental for a filling fraction approaching 60%.
Packing structures with sticks are also considered in order to make clearer the role of channels arising from the arrangement of particles introducing long walls. Table VI summarizes the results of simulations. It is manifest that the systems with alignment lead to the fastest infiltrations. Furthermore, within this group the dynamics tend to be faster in the presence of longer sticks. These results show that fibers favor the formation of preferential infiltration paths. This comes out also from the high values of the error ∆z on the average penetration depth. With orientation disorder, results of rhombs and fibers are still comparable (cf. Tab. V). Sticks introduce more channels that are filled easily, inducing stronger accelerations. This appears also from the fluctuations of Ca and Re that turn out to be more marked, displaying also higher peaks (results not shown for brevity). Nevertheless, also for this set of data, capillary forces remain the dominant ones. It is interesting to note that the simulations with misalignment and tilt can lead to quite fast infiltrations, even though the number of channels is lower. A possible explanation is that orientation disorder for fibers can be regarded as equivalent to the presence of small particles when 2D cross-sections are considered (see Fig. 2). As a consequence, when the meniscus slides along the walls introduced by fibers, it is likely to encounter the surface of other particles before the interface becomes flat and pinning occurs. So, the results for fibers confirm the advantages identified for pore configurations associated with bimodal size distributions. The permeabilities for packings with fibers are more informative on the performace of the porous systems for capillary infiltration than in the presence of rhombs. It is also interesting to consider the influence of the surface area of the solid phase (Bear, 1972), i.e. the length of its boundary for 2D structures. In principle, for a given volume, spherical shapes tend to minimize the surface area while thin oblate spheroids mazimize it. It is found that for systems obtained from spheres and rhombs, the boundary length of the solid phase is higher for the systems composed of small particles, in particular for spheres. The boundary of the solid phase is longer in the case of aligned fibers. With orientation disorder, the length of the solid boundary is still longer than that of systems with rhombs. By considering for the solid phase its surface area divided by its volume, known as specific surface (Bear, 1972), it is found that the systems with fibers exhibits higher values than bimodal systems on average. The specific surface is in use for the characterization of powders. For our systems, the specific surface is lower for rhombs because it is easier to realize compact packings with coarse particles. The highest values are obtained with fibers that have however the advantage to favor the formation of channels.
Simulations in the presence of surface growth are performed for packing structures composed of rhombs (see Tabs. II and V). A typical configuration is shown in Fig. 2. Figure 6 shows examples of infiltration kinetics. At the interface between the solid and the wetting fluid the front retreats as the surface grows. This effect is more important when agglomerates of particles merge. In the analysis of the data, we thus consider an average over the highest values for the position of the front. Fits to Eq. 2 allow to obtain the initial effective radius r 0 . The effective radius used in order to determine the characteristic numbers is set to the average value: r eff = r 0 /2. The question arises whether infiltrations stop because of pore closure or pinning. Examination of the fastest dynamics indicates that the cause due to pore closure is as frequent as that for the barrier introduced by pinning. This implies that surface growth enhances the phenomenon of pinning. Indeed, it is found that in the period considered for the fits, the infiltrations without surface reaction are on average 5.7% faster. The resistance of surface growth is more marked for the fastest dynamics, i.e. with higher orientation disorder. At the point where infiltration stops according to the theoretical model of Eq. 2, the results with surface growth are comparable to that with inert boundaries. But it is to be noted that r 0 turns out to be underestimated in particular for the longer infiltrations where the front can overcome pinning barriers. More details for the results are summarized in Tab. VII. It can be seen that orientation disorder is associated with faster infiltrations. Bimodal size distributions lead to better results for aligned and misaligned particles. This is no more true with the addition of tilt. In any case, it can be expected that, in 3D porous structures, bimodal size distributions still present more advantages. In Tab. VII, the infiltration velocities are higher than in Tab. V. With surface reaction, the dynamics is shorter and the initial fast infiltration has more weight. It should be kept in mind that now the dominant effect is the infiltration interruption directly caused or induced by surface growth. With surface growth allowed, Ca, Q 4 and Re increase appreciably (see Tab. V). Furthermore, their fluctuations are by a factor of 2 smaller than for inert boundaries. Q 4 gains on average one order of magnitude, but the effects of capillary and vis-cous forces remain more significant. This arises also from the characteristic numbers involving the properties of surface reaction Q 1 , Q 2 and Q 3 . These quantities describe better the systems because infiltration stops under the action of surface growth. Q 1 is more accurate because the effects of capillary forces are second in importance even though their relevance weakens with respect to viscosity and inertia in the presence of surface growth. Only the fluctuations of Q 3 in a few cases can deviate from the average value more than one order of magnitude.

Comparison with experiments
We now want to carry out simulations for typical experimental results with pure Si (Eustathopoulos, 2015;Israel et al., 2010). The initial porosity of the preform can be assumed to be 0.3 with initial average pore diameter of d = 10 µm. The pore size distribution is assumed to be quite narrow, with pores mainly in the range 7 − 13 µm (diameters). The preform is composed of graphite particles with average size of 70 µm. The reaction-rate constant is k = 4 ·10 −8 m/s (Messner & Chiang, 1990). The other parameters of Si are ρ = 2.53 · 10 3 kg/m 3 for the density, µ = 0.94 · 10 −3 N·s/m 2 for the viscosity and γ = 0.86 N/m for the surface tension (Einset, 1996). For the infiltration velocity we use v = 5 µm/s (Eustathopoulos, 2015). Thus, the holding time for the experiment can be estimated to be t pc = 125 sec. From the average pore radius it follows that at this time, indicated by t pc , pore closure occurs and the infiltration stops. The maximal infiltration depth results to be z max = 0.625 mm. The effective radius is estimated using Dullien model as in the article by Einset (1996). By taking into account the thickening of the surface of 2.5 µm leading to the average pore diameter d/2, it is found that an estimate for the effective radius is r eff = 0.25 µm. Indicatively, in industrial practice targeted by us (e.g., manufacturing of burners and radiation plates), the C preforms are obtained starting from powders with bimodal size distributions. The size of the largest particles can vary in the range of 30 − 70 µm. The size of the smallest particles can be assumed to be 5 µm. Given the above experimental conditions, for the simulations we consider systems N x = 1500 lu long and N y = 511 lu wide. Since the width is twice longer than before, roughly speaking the infiltration velocity may double, as suggested by Eqs. 2 and 3. The length of the samples is again L = N x /2, placed in the simulation domain as described in Sec. 4. The total number of timesteps is now set to N t = 4 · 10 6 ts. The surface grows after every N t /10 timesteps. In simulations, the reaction-rate constant is thus given by k = 2.5 · 10 −6 lu/ts. All the other parameters for the LB simulations remain unchanged. The porous structures are obtained from the packing of rhombs with average side and height corresponding to 36 lu and 31 lu. Both misalignment and tilt are allowed. The resolution of the systems does not allow to include small particles as in industrial applications. The porosity of the resulting sys-tem is 0.64. In this case, the sides of the unit cube defining the porous structure are subdivided into 512 bins. Simulations are performed for 30 sections in 2D. The correspondence with experimental results is established by imposing the experimental values for z max and t pc . Figure 7 shows the kinetics of infiltration as obtained from simulations. As done before, an average is taken over the highest values for the positions of the fronts. It can be seen that the model predicts a more abrupt behavior for the interruption of penetration, as observed also in experiments (Israel et al., 2010). It turns out that after an initial fast dynamics, the infiltration velocity decreases and in the last stage the maximal depth is attained slowly. For the sake of clarity, we consider suitable averages in the neighborhood of the point where pore closure occurs and in the ensuing interval of time where a plateau is expected to form. It is found that for 63% of the systems the infiltrated distance increases more than 10% after the time t pc where pore closure is predicted by the model for the data of Fig. 7. On average, the further advancement amounts to 45%. For the remaining 37% of the systems after the time t pc the fronts travel on average an additional distance corresponding to 0.6%. It is verified that in this case the flow interruption can be ascribed to pinning. As a result, the value of t pc determined by the model is practically associated with the first pinning barrier. When the plateau start forming the effect of surface growth seems to remain mainly indirect. It is to be noted that the simulations reveal the tendency that the interruption of invasion is affected by the combination of the effects of surface growth, tortuosity and interface dynamics. In principle, it is reasonable to presume that the infiltration may stop before the time required for the obstruction of the average radius since the process of infiltration is described by an effective radius. However, we find that in the simulations the effect of pinning is more critical than expected in experiments. In order to elucidate the role of the effective radius more experiments are also needed for given well-known reaction-rate constants (Israel et al., 2010).
In Tab. VIII the characteristic numbers of simulations and experiments are compared. The objective is to assess the quality of the equivalence between experimental and simulation conditions. It appears that the effects of capillary forces are too weak with respect to the effects of surface growth and viscous forces. For the capillary number it is well known that high values satisfy Ca > 10 −2 . For this reason, it could be questioned if the difference for the orders of magnitude of Q 1 is of particular relevance. For sure, our approach allows to establish an equivalence for the dynamics of the average penetration depth. Instead, the details of the invasion process can not be taken for granted. For example, for estimates based on the average radius of the porous structure, in experiments, when the infiltration stops because of pore closure, the ratio of the infiltrated distance to the thickness of the reaction-formed SiC is 125. In simulations, it is obtained 75. As discussed above, the relation between the average and effective radii deserve attention but a major hurdle in order to achieve better results is represented by the phenomenon of pinning. These simulations also prove that a higher resolution is not sufficient in order to reduce pinning. Again, we see that the structure is more important than the average radius, of course increasing with a higher resolution. Under the conditions of these simulations, the average size of the particles is 150 µm while the target was 70 µm for an equivalence with a stronger direct correspondence for the porous structure. It is not easy to improve the comparison because with packings composed of smaller particles, pinning would become even more important. Furthermore, if we use the data associated with the first point in Fig. 7 exceeding the average of z max over the last 10 frames taken into account for the fit, the improvement is relative. The reason is that the velocity decreases but the initial effective radius r 0 increases. The result is that Q 1 and Ca remain of the same order. Thus, it is not straightforward in the simulations to reduce the velocity, to keep r 0 small and to weaken the phenomenon of pinning at the same time.

CONCLUSIONS
In this work, LB simulations are proposed for the study of capillary infiltration into porous structures. Optimization principles for the manufacturing of ceramic materials through reactive melt infiltration are at the center of our attention. In this aim, first are considered porous systems realized with the packing of particles with different properties (shape, size distribution, orientation disorder). Without surface growth, our analysis shows that the invading front tend to select preferential pathways. This can result in the formation of fingers. The data also indicate that bimodal size distributions display better infiltrations than monomodal distributions. The improvement is found to be 4.1% on average (see Tabs. IV and V). Pinning of the contact line appears as the primary factor affecting liquid penetration. Structural disorder can reduce its effect. Interestingly, bimodal size distributions present advantages in order to limit pinning. Larger faceted particles introduce walls giving rise to channels. As a result, the meniscus travels longer distances before the loss of its curvature. The probability to encounter the surface of other particles is even higher with the addition of smaller particles. The results for packings with fibers support this view. It is also verified that inertia effects induced by transient accelerations remain of secondary importance. Surface growth results in a clear interruption of the flow. This can occur because pinning turns out to be enhanced or for the pore-closing phenomenon. The retardation effects are more marked for the fast dynamics with a higher degree of orientation disorder. Again, the results point out that bimodal size distributions appear to be able of accommodating optimal infiltration. Simulations with higher resolution are also performed for a comparison with experiments. The relative strength of capillary forces turns out not to be strong enough. As a result, the infiltration dynamics is not reproduced in full details. Furthermore, the effects of surface growth for the interruption of infiltration are basically indirect in this case. We identify two aspects significant for further advances. First, the relation for the interruption of infiltration with the physical (e.g., average and minimum) and effective radii needs to be elucidated. Second, longer infiltrated lengths should be realized in simulations. A major obstacle to im-provements resides in the over-occurrence of pinning. This situation could be improved by considering other combinations of surface tension, viscosity and, especially, contact angle (Chibbaro, Costa, et al., 2009). This could also be an extension of this work to Si alloys. Finally, our study develops a simulation work relevant for increasing the quality of ceramics obtained from LSI. Despite a poor statistics and noisy results, operational guidelines are highlighted. Further research is necessary in order to sharpen the conditions determining the optimal configuration for the porosity. Experimental tests would also allow to understand the robustness of our findings based on simplified systems.