2.1 Governing equations, boundary and matching conditions

When the Reynolds number of channels and the tows permeability are small, the coupling between the free-fluid and porous sub-domains depicted in Figure 1 can be defined by a Stokes-Darcy formulation as given by equations (1)-(4).

For the channel (Stokes domain):

For the porous medium (Darcy flow in the principal directions of permeability):

On the other hand, the architecture of the porous media (warps and weft) is considered as a bank of aligned micro-cylinders, and the main permeabilities can be computed using the model proposed by Gebart [²²], equations (5),(6).

Where parameters c, c ₁ and V _f,max depend on the type of array. For a hexagonal array, these parameters are given in [²²].

At the channel-tows interface, the matching conditions are given by equations (7)-(10).

Continuity of normal velocities:

Jump of tangential velocities: this condition has been widely discussed in the literature [²³,²⁴], but the most common is the slip condition of Beavers-Joseph that reads [²⁵,²⁶]:

When the ratio between the height of the porous medium and the square root of permeability is very high, a good approximation for the slip coefficient, γ, is [²⁷]:

A deeper discussion about the Beavers-Joseph condition for this application and its numerical treatment can be found in [²¹].

Continuity of surface tractions:

For the problem represented in Figure 1, the boundary conditions can be classified into three types:

For the inlet and outlet of the channel domain, uniform horizontal tractions are considered, equations (11), (12):

No penetration condition at the Darcy domains (tows) is given by equation (13):

Free-surface conditions at the liquid-air interfaces are stated in equations (14)-(16):

At the liquid-air interface corresponding to the bubble front when void migration occurs, capillary pressure is computed as p _cap = −σκ = −σ(∇ · ˆn), where the curvature at the liquid-air interface, 𝜅, is numerically computed as shown in [²¹,²⁸]. On the other hand, in the porous medium, the capillary pressure for the warps is calculated with equation (17) [²⁹]:

Where the equivalent capillary radius, R _ec , is given by equation (18):

For 0 ≤ φ < π/2, equations for A_{int ,s} and C_int,s are (19) and (20), respectively [²¹,³⁰]:

If the impregnation occurs perpendicular to the fibers (φ = π/2), equations (21) and (22) apply:

For the wefts, capillary pressure is estimated as posed in equations (23)-(27) [³¹]:

2.2 Integral equation formulations and numerical techniques

The numerical treatment of the governing equations, boundary conditions, and matching conditions of the Stokes-Darcy approach was adopted from [²¹]; the present section summarizes the main points of this treatment. The Boundary Element Method (BEM) is used here to solve the governing equations. The boundary integral formulations for the Stokes [³²] and Darcy anisotropic [³³-³⁵] equations are given by (28) and (29), respectively:

where: c _ij = (α/2π)δ _ij , with α as the solid angle at the source point, whose value is α = π for points located over a smooth contour. For points located inside the domain, c _ij = δ _ij . In a similar fashion, c (ξe) is equal to 1/2 for points over a smooth surface and to 1, for points inside the domain. On the other hand, the superscript “e” represents the Equivalent Isotropic System (EIS). The integral kernels in (28) and (29) are given in trms of the corresponding fundamental solutions equations [(30) to (36):

where ξ (ξ^e in the EIS system) and y (ye in the EIS system) are the source and field points, respectively, with the integrals containing the fundamental solutions U_j ⁱ (ξ, y) (Stokes) and p^∗ (ξ^e, y^e) (Darcy) corresponding to the Single Layer Potential (SLP), and those ones containing Kij (ξ, y) (Stokes) and p^∗ (ξ^e, y^e) (Darcy) to the Double Layer Potential (DLP). The factor f ^e comes from the transformation of the normal vectors, as shown in [²¹].

The following characteristics of the numerical setup are deemed in the present work: second-order isoparametric discretization of the boundary and physical variables, non-continuous shape functions for the corner elements, numerical calculation of integrals using Gaussian quadrature, rigid body motion principle [³⁴,³⁶] and Telles transformation [³⁷] for the numerical treatment of strong and weak singularities, respectively. Details about the numerical treatment of the Beavers-Joseph condition, equation (8), are given in [²¹]. By imposing the boundary and matching conditions, a well-posed system of equations is obtained, which is solved using a Singular Value Decomposition (SVD) algorithm [³⁸,³⁹].

Tracking of the free surfaces is done by implementing first-order Euler integration of equation (14). Time step selection according to physical, numerical and geometrical constraints, the implementation of remeshing and smoothing algorithms depending on the pressure and capillary forces, and the numerical calculation of surface curvatures, are critical issues for such tracking technique. This is thoroughly discussed in Appendix of [⁴⁰].

The fluid velocity inside the channel is computed using (28) with c _ij = δ _ij . The integral representation for the pressure is used to compute the pressure field inside the channel, as shown in equation (37) [³⁴]:

Where is the fundamental solution for the pressure equation, and is the gradient of The pressure field inside the tows is obtained from (29) with c (ξ) = 1, while the pressure gradient is found by taking the directional derivative of (29), as shown in equation (38):

The velocities in the interior points of the porous media are given by the Darcy’s law, in terms of the pressure gradient, equation (38).

3.1 Validation of the BEM method for Stokes-Darcy problems

To verify the developed BEM code, the coupled Stokes-Darcy problem of Figure 2 that admits an analytical solution is considered. In this problem, botgh the free-fluid (Stokes) and porous medium (Darcy) sub-domains are totally saturated. For the porous medium, porosity is set to ε _t = 0.6, leading to a slip coefficient of γ = 1.29 according to equation (9). For this particular problem, an analytical solution for the pressure and velocity field in terms of power series was developed in [²¹], which is reproduced in the equations (39)-(52).

For the anisotropic porous medium (0 ≤ y ≤ h _d ):

For the channel domain (−h _s ≤ y ≤ 0):

The Fourier coefficients of (39) and (46) are given by:

The L² relative error norms between the analytical and the BEM solutions for the pressure, p, and velocities, u _x and u _y , are shown in Table 1. Three factors are considered for the analysis of these errors: mesh size (h), anisotropic ratio (K₁ K₂) and slip coefficient (γ). The mesh size is reported as h = e/H, where H = h _d + h _s , is the total height of the domain (See Figure 2] and e represents the element size.

According to Table 1, the maximum error for a determined mesh size is always obtained for the horizontal velocity, u _x , when K₁K₂ = 0.1 and γ = 0.01, but this error decreases as the mesh is finer. In general, the refinement of the mesh leads to the reduction of the L² error in all situations, which means that the BEM code converges to the analytical solution; however, the order of convergence is not necessarily the same in all cases. The vertical velocities, u _y , are not considered in the present analysis since they are very close to zero; however, as observed in Table 1, the accuracy corresponding to u _y is acceptable enough in all situations.

The slip coefficients of coupled Stokes-Darcy simulations of the following sections are calculated using the equation (9) and have an order of magnitude of O (0). Considering the results achieved for γ = 1.29 in the Table 1, it is expected to obtain an acceptable accuracy for those simulations, even for the coarser mesh evaluated here, i.e., h = 1 × 10−1. In order to define the suitable mesh size for simulations of the following sections, both the L² relative error norms and the run times are deemed. The CPU times obtained with an Intel Pentium 2.1 GHz and 2GB of memory are shown in Figure 3. Let us consider the change of mesh size from h = 1 × 10−1 to h = 5 × 10⁻², where average reductions of the L² relative error norm of 33.7% and 32.3% are obtained for pressure and horizontal velocity, respectively, leading to an increment in the CPU time of 6.49 times; in such case, mesh refinement could be considered appropriate. On the other hand, when mesh-size changes from h = 1×10−1 to h = 3.3×10⁻², the increment of the CPU time is 23.40 times, while error reduces by 48.2% for the pressure and 54.3% for the horizontal velocity, which can be considered impractical in this case. Taking this in mind, the mesh-size is taken as h = 5 × 10⁻² in the following simulations.

The slip coefficients of coupled Stokes-Darcy simulations of the following sections are calculated using the equation (9) and have an order of magnitude of 𝒪(0). Considering the results achieved for 𝛾=1.29 in the Table 1, it is expected to obtain an acceptable accuracy for those simulations, even for the coarser mesh evaluated here, i.e., h = 1 × 10⁻¹. In order to define the suitable mesh size for simulations of the following sections, both the L² relative error norms and the run times are deemed. The CPU times obtained with an Intel Pentium 2.1 GHz and 2GB of memory are shown in Figure 3. Let us consider the change of mesh size from h = 1 × 10⁻¹ to h = 5 × 10⁻², where average reductions of the L² relative error norm of 33.7% and 32.3% are obtained for pressure and horizontal velocity, respectively, leading to an increment in the CPU time of 6.49 times; in such case, mesh refinement could be considered appropriate. On the other hand, when mesh-size changes from h = 1 × 10⁻¹ to h = 3.3 × 10⁻², the increment of the CPU time is 23.40 times, while error reduces by 48.2 % for the pressure and 54.3 % for the horizontal velocity, which can be considered impractical in this case. Taking this in mind, the mesh-size is taken as h = 5 × 10⁻² in the following simulations.

3.2 Comparison between the present Stokes-Darcy approach and the classical approach

The main differences between the classical approach, where the channel flow is not modeled and a uniform pressure is supposed in the channel, and the present approach, where the channel flow is modeled by the Stokes equation and it is prescribed a pressure gradient, can be noticed when comparing simulations of Figures 4a-c and 5a-k. For each filling instant, the time normalized with the total time of simulation (t/t _sim ), warps saturation (S _warps ), weft saturation (S _weft ) and total tows saturation (S _t ) are shown; the x and y coordinates are non-dimensionalized as x ^∗ = x/L _RUC and y ^∗ = y/L _RUC . The geometric and material inputs of both simulations are shown in Table 2, where it is also observed that for the simulation with the classical approach [Figures 4a-c) a uniform channel pressure of ⟨P _g ⟩ ^g = 122kPa is prescribed, while for the simulation with the present approach [Figures 5a-k), inlet and outlet pressures of p_in = 125.5kPa and p_out = 118.5kPa are considered, which originates a pressure gradient along the RUC of △P/Δx = 5.83 × 10³kPa/m, with an average pressure of ⟨P _g ⟩ ^g = 122kPa, agreeing with the average pressure of the other approach. For simulation of Figures 4a-c, constant atmospheric pressure is considered for the air given the full air dissolution assumption; on the other hand, in the current approach, Figures 5a-k, a vacuum pressure of P_vac = −75KPa is demeed, varying in the weft following the ideal gas law.

For the classical approach, as observed in Figures 4a-c, the fluid flow motion in warps and weft is uniform, namely, towards the edges and the center of the RUC; this occurs because the pressure is the same at the channel-tows interface. Total saturation is possible with this approach [Figure 4c), with the complete filling of warps happening first [Figure 4b). Thus, according to this approach, full air dissolution is possible. However, this is ideal case does not necessarily correspond to a real application.

On the other hand, for the present Stokes-Darcy approach, considering full air compression in the weft, three stages can be clearly identified for the void in Figure 5a-k: compression, mobilization and migration. In the compression stage, the fluid fronts in the warps are not parallel to inferior and superior edges of the RUC, and the fluid front in weft is barely decentered (see datatips in Figure 5a), which is caused by the imposition of a pressure gradient along the flow direction in the channel; this phenomenon is more perceptible for greater pressure gradients. The filling of the warps occurs at 2.7% of the total simulation time [Figure 5b). When the equilibrium condition is achieved for the trapped bubble, air compression ceases, and void mobilization starts (green line in Figure 5c).

From this time instant until the onset of void migration, weft saturation, S_weft, is essentially constant, reaching an equilibrium saturation where void shifts towards the right extreme of the weft at constant volume [Figures 5d to 5g). Small changes of S_weft during this process of void displacement can be present by numerical errors during the fluid front tracking, but they can be neglected.

The compression phenomenon lasts 13.3% of the total simulation time, whereas the void mobilization takes almost the remaining time, namely, 86.7%. When the void reaches the right extreme of weft, Figure 5f, normalized time is t/_tsim = 1 − 2.80 × 10⁻⁵, which is very close to one because the time elapsed from this instant until the end of the simulation (the time of partial migration of the bubble), is much shorter than the times corresponding to void compression and void mobilization. The process of void migration can be seen in Figs. 5h to 5k. The balance between capillary forces, which depend on the surface tension (λ) and on the curvature of the bubble surface (κ), and pressure forces of the compressed air, which changes with the bubble expansion or shrinking, determines the dynamic evolution of the intra-tow void.

3.3 Velocity field and streamlines

In Figures 6a-c, the velocity fields and streamlines when the bubble is compressing inside the weft and the warps are totally saturated are shown; that is, there is no mass transfer from the channels into the warps; some figures are off-scale for visualization purposes. Several features can be identified in Figures 6a-c for the channel domain. Firstly, the velocity vectors at the interface with the weft are almost tangential to such interface [Figures 6b and 6c), which means that the magnitude of the normal penetration velocity into the weft is very small in comparison with the tangential velocity.

It is also remarkable the increment in the magnitude of the velocity with the reduction of the inter-tow space [Figures 6a-c) because this indicates that the effect of this reduction in the velocity field is more relevant than the effect of the liquid absorption into the tows that reduces the velocity magnitude along the RUC. Regarding the streamlines, the BEM code predicts that the inferior ones (nearby to the lower edge of Figures 6a) tend to concentrate at the symmetric boundary as the inter-tow space reduces. In general, when a streamline traverses the whole RUC, it ends at a lower vertical position with respect to its starting vertical position and this is more notorious as the streamlines are closer to the symmetric boundary.

On the other hand, for the weft domain [Figure 6d-h), it can be appreciated that the velocity at the channel-weft interface is essentially normal to such interface; this means that the tangential Darcy velocity appearing in equation (8), is negligible. It is also worth noting, by comparing the magnitudes of the velocity vectors along the channel-weft interface [Figure 6d), that the mass transfer is greater as the inter-tow space is lower, in such a way that the normal velocities at the channel-weft interface are larger in points close to the center of the weft, whereas they are smaller in the neighborhoods of the right and left edges of the weft. The mass transfer in the left half of the weft is barely greater than the mass transfer in the right half, but this difference is enough to generate a decentered bubble, as it was shown in Figure 5a. As it can be appreciated in Figure 6h, most of the streamlines with positions of starting points beyond the horizontal limits defined by the fluid front, converge in the right and left extremes of such a fluid front, and the starting points of these streamlines coincide with the points of low mass transfer, as it can be noticed by comparing Figures 6d and 6h.

The velocity fields and streamlines when the bubble migration occurs are represented in Figures 7a-g. For the channel domain [Figures 7a-b), it can be appreciated that the highest velocities are obtained in the region of bubble migration, which means that the average air migration velocity, , is greater than the average liquid velocity in the channel, , for this particular case. As observed in Figure 7a, the streamlines starting in the inlet boundary encounter with those ones starting in the contour of the partially escaped bubble and this causes the accumulation of the former ones in the right inferior zone of the RUC. In Figures 7a, likewise to the formerly commented case of void compression in the weft [Figures 6a), a high density of streamlines can be observed in the symmetric boundary of the RUC when the inter-tow distance is the smallest, namely, in the half of the RUC. For the weft domain [Figures 7c-f), the velocity magnitudes are larger in points whose horizontal position is between the horizontal limits defined by the fluid front, as in the last case analyzed [Figures 6d-g).

In this case, most of the streamlines starting in the zone of low mass transfer, where the normal interfacial velocities are smaller, converge in the left extreme of the fluid front (See Figures 7c). According to Figure 7g, in the right extreme of the weft, some streamlines are very short since the interface channel-weft, where the streamlines start, is adjacent to the fluid front, where the streamlines finish.