Stokes–Brinkman formulation for prediction of void formation in dual-scale fibrous reinforcements: a BEM/DR-BEM simulation

A numerical study of voids formation in dual-scale fibrous reinforcements is presented. Flow fields in channels (Stokes) and tows (Brinkman) are solved via direct Boundary Element Method and Dual Reciprocity Boundary Element Method, respectively. The present approach uses only boundary discretization and Dual Reciprocity domain interpolation, which is advantageous in this type of moving boundary problems and leads to an accurate representation of the moving interfaces. A problem admitting analytical solution, previously solved by domain-meshing techniques, is used to assess the accuracy of the present approach, obtaining satisfactory results. Fillings of Representative Unitary Cells at constant pressure are considered to analyze the influence of capillary ratio, jump stress coefficient and two formulations (Stokes–Brinkman and Stokes–Darcy) on the filling process, void formation and void characterization. Filling times, fluid front shapes, void size and shape, time and space evolution of the saturation, are influenced by these parameters, but voids location is not.


Introduction
In dual-scale fibrous reinforcements used in composites processing, the permeability of the tows is several orders of magnitudes lower than the permeability of the inter-tow channels. This causes imbalances between the flow inside the tows and the flow in the channels, which can lead to the formation of voids by mechanical entrapment of air. Once these voids are formed, several phenomena can take place, among others: bubble compression, displacement, migration and splitting.
Under the assumptions that the solid phase is stationary, the fluid is Newtonian and incompressible and neglecting the inertial effects, the Navier-Stokes system of equations governing the flow inside the pores of the porous medium (microscale) reduces to the Stokes equation. By an assembled average including both phases of the porous medium (liquid and solid), the Stokes system of equations reduces to Darcy's law, where the velocity (Darcy velocity) is linearly proportional to the pressure gradient, with proportionality constant given by the permeability, Darcy's law was established empirically by Darcy (1856) and later derived formally by different authors using several approaches (among others see Whitaker [66]).
Considering the shearing effect between the fluid and the walls of the porous medium, Brinkman (1947) added the second-order derivatives of the velocity to the Darcy equation, resulting in the Brinkman equation. The averages implicit in the Brinkman equation can be viewed as averages over an ensemble of different scales of the porous medium that interpolate between the Stokes and Darcy equations. On small length scales, the pressure gradient balances the Laplacian of the velocity and the flow is essentially viscous, Stokes flow, while over larger length scales, the velocity is slowly varying and the pressure gradient balances the average velocity as it does in Darcy's law (for more details see Durlofsky and Brady [9]). According to Ochoa-Tapia and Whitaker [40], Brinkman approximation can be applied when the following three length scale constrains are simultaneously satisfied: r 0 2 / l ε l p 1 , r 0 2 / (l ε l v ) 1, l f r 0 , where l ε , l p , l v and l f are the characteristic lengths associated to the porosity, pressure, velocity gradients and fluid phase, respectively, and r 0 is the characteristic length scale of the Representative Elementary Volume (REV). Krotkiewski et al. [28] reported a direct numerical simulation of the flow field in homogeneous two dimensional porous media with characteristic length L and permeability K , concluding that the Stokes solution is dominant for K /L 2 ≥ 10, Darcy law is representative of the flow field if K /L 2 ≤ 10 −4 , while for 10 −4 ≤ K /L 2 ≤ 10 the Brinkman approximation should be used to account for the transition between both flow regimes, Stokes and Darcy.
The Brinkman equation includes an effective viscosity term, μ e f f , to consider the viscous diffusion not deemed in the Darcy law. In order to explain the meaning of this term, it is necessary to consider that the volume averaging method allows describing the porous medium flow in terms of averages of the local quantities by means of the Darcy-Brinkman equation, which is given by (3) considering μ e f f = μ, where μ is the real fluid viscosity. This last equality is valid provided that a non-slip condition on the interfaces between the fluid and solid phases of the porous medium is considered, as done in the traditional volume-averaging method. However, the non-coincidence between μ e f f and μ has been demonstrated in several experimental, numerical and theoretical works, suggesting that the non-slip condition is not necessarily valid in all cases. For instance, Givler and Altobelli [17] experimentally found that μ e f f is 7.5 times the value of μ for high-porous open cell foams and moderate Reynolds numbers, whereas Starov and Zhdanov [70] studied the dependency of μ e f f on the porosity and particle size in porous media composed of equally sized spherical particles, finding that μ e f f can be lower or larger than μ. On the other hand, the numerical flow simulations conducted by [4] in 3D regular arrays of cubes showed that μ e f f < μ, whereas those ones executed by [34] indicated that μ e f f > μ in order to match the Brinkman equation with the numerical solutions in the boundary layer developed in the porous medium domain when it is in contact with a free-fluid domain. A very illustrative theoretical work was published by [63], where it was demonstrated that the effective viscosity, μ e f f , is different to the fluid viscosity, μ, when a slip condition at the fluidsolid interface of the porous medium is considered. Using an up-scaling procedure, a boundary-value problem to compute μ e f f was obtained in such work, achieving important conclusions. In general, the effective viscosity is different to the fluid viscosity, i.e, μ e f f = μ, because it is a term that 'absorbs' the microscale variations of the velocity gradient when a variationless velocity gradient model is used at macroscopic scale. On the other hand, as the prescribed slip coefficient at the fluid-solid interface increases, the boundary layer thickness decreases and a non-slip condition tends to be reached, in such a way that in the limit when the slip coefficient tends to infinity, μ e f f = μ.
In the Stokes-Brinkman approach the matching conditions at the interface channels-tows are defined in terms of the corresponding surface velocity and traction vectors. Continuity of the velocity field is always considered and two types of conditions could be employed for the tractions: continuous [39] and jump stress [40,41]. The continuity of stress was initially implemented in [39,48], however a jump stress matching condition was proposed later in [40] based on the non-local form of the volume-averaged momentum equation for the analysis of the region between the channel and the porous medium. By means of experiments of unidirectional flow in parallel domains, Ochoa-Tapia and Whitaker [41] found that the jump stress tensor across the interface for isotropic porous media, σ i jn j , is linearly proportional to the interface surface velocity u j , i.e. σ i j n j = μT i j u j / √ K , where μ stands for the fluid viscosity, while the second-order tensor T i j is given by βδ i j , with β as the jump stress coefficient ranging between −1 and 1.47; for more details see [23]. However, Angot [2], demonstrated mathematically that the well-posedness of the Stokes-Brinkman problem is only possible when β ≥ 0. The theoretical estimation of the value of β was developed in [64,65], where it was established that β depends on the porosity of the porous medium in the inter-region channel-tow, ε ct .
The estimation of the effective viscosity term, μ e f f , in the Brinkman equation is known to depend on the geometry of the porous medium and the flow itself. As suggested in the literature [23,56,57], an acceptable approximation for μ e f f in dual-scale fibrous reinforcements is to consider that μ e f f = μ when the continuous stress condition is used, whereas μ e f f = μ/ε t is suitable when the jump stress condition is employed, with ε t as the tow porosity. In a recent work [23], focusing in the prediction of the effective saturated permeability of dual-scale reinforcements for unidirectional flow (where a value of β = 0.7 was used), it was found that the type of stress matching condition (continuity or jump) has a strong influence on the boundary layer thickness in the porous medium, but not on the effective saturated permeability, which is dependent on the pore geometry. In this point, it is relevant to mention that the study of the pore geometrical dependency of the effective permeability could be a complicated and computational expensive task, but some recent efforts have been addressed to reduce the computational cost without compromising the accuracy. For instance, a multiscale framework relating the effective permeability with some microstructural attributes extracted by X-ray tomography was proposed in [55], where several computational techniques were efficiently combined (Level Set, Graph Theory and Lattice Boltzman/Finite Element) to determine the tortuosity, porosity and homogenized effective permeability at specimen scale. The particular problem solved in [23] by FEM is also considered here for validation of the numerical technique implemented.
In the present work, simulations of fillings of RUC's (Representative Unitary Cells) of dual-scale fibrous reinforcements are carried out using a Boundary Element/Dual Reciprocity numerical approach (BEM/DR-BEM). The flow in the channel, Stokes flow, is given by its boundary-only integral formulation, Green's formula, and numerically treated by direct BEM, while the flow inside the tows, Brinkman flow, is given by a boundary-domain integral formulation in terms of the Stokes fundamental solutions, where the domain integral is transformed into a boundary integral using the Dual Reciprocity Boundary Element Method (DR-BEM). The present numerical approach is validated by comparison with a benchmark analytical solution used previously in the literature to assess the robustness of FEM-based numerical solutions of flow problems in dual-scale porous media [23,58].
The developed BEM/DR-BEM numerical scheme is then used to simulate the simultaneous filling of channels and tows inside the RUC at constant pressure regime with the purpose of analyzing the influence of the stress matching conditions on the void formation at several values of the capillary ratio, C cap = P cap,max /P in , where P cap,max and P in stand for the maximum capillary pressure and the inlet pressure, respectively. This type of filling problem has been previously considered in the literature using different formulations and numerical techniques. For instance, Gourichon et al. [19] studied the influence of the RUC porosity in the formation of voids using a Darcy-Darcy formulation and the FEM/CV conforming method. On the other hand, Schell et al. [49] studied the influence of the tow porosity, ε t , on the final void content using the same formulation and numerical technique as in [19], while the problem of unidirectional filling of circular tows in cylindrical coordinates was analyzed in [68] using a Stokes-Brinkman formulation and the Finite Volume method, where it was studied the effect of the filling velocity, resin viscosity, inter-tow dimension and intra-tow dimension on the shape of the fluid front.
In comparison with domain-meshing numerical techniques previously used in the literature for the simulation of the filling process in dual scale fibrous reinforcements based on the Stokes-Brinkman formulation [11,23,50,51,58], two main advantages can be identified by the use of the BEM/DR-BEM approach employed in this work. Firstly, the use of BEM-based techniques does not imply any mesh discretization of the problem domain, which is not a trivial task in a domain that is continuously changing in size and shape as in the present case. Secondly, the tracking of the fluid front is directly carried out by an Euler integration of the kinematic condition at the moving interface, which assures a higher order accuracy on the prediction of the fluid front shape without the need of implementing any additional reconstruction algorithm of the moving boundaries, contrary to the case of the Volume of Fluid (VOF) [14,15] or the Level Set Method [12,13,53].
Another important contribution of this work is the study of the influence of two types of stress matching conditions (continuous and jump) on the void formation. Other works have analyzed the influence of these conditions on the effective saturated permeability [23,58], but not on the size, shape and location of the voids formed by mechanical entrapment of air. Additionally, the processes of the void formation using the Stokes-Brinkman and Stokes-Darcy formulations are compared each other, which has not been reported before to the best of our knowledge. Finally, the consideration of a flow direction-dependent capillary pressure in the tows without experimental factors and of the curvature-dependent surface traction effects for the fluid fronts in the channels are ones of the main differences between the present work and previous publications in the same field [1,7,49,59].

Governing equations, boundary and matching conditions
In the Stokes-Brinkman formulation, the governing equations at each medium are defined as: Mass conservation (For all domains): Momentum for the Stokes domain (Channels flow): Momentum for the Brinkman domain (Porous media or tows  flow): where in the right hand side of (3) the double index notation is not considered for i = 1, 2. Here u i , p, μ, μ e f f and K i represent the velocity vector, pressure, liquid viscosity, effective viscosity and main permeabilities, respectively. Let us define the following non-dimensional variables for the length, velocity, time and pressure [46,47]: For all domains: For the Stokes domain (Channels flow): For the Brinkman domain (Porous media or tows flow): where L * and U max are the characteristic length and the maximum velocity of the problem, respectively. In terms of these characteristic values, the non-dimensional form of the governing Eqs. (1), (2), (3) can be written as: where μ e f f is the inverse of the Darcian number, Da * i , in the principal direction i. The non-dimensional matching conditions for the Stokes-Brinkman problem are as follows [2]: • Continuity of velocities: • Normal and tangential component of the jump stress: where s and b represent the Stokes and Brinkman domain, respectively, andt / μ e f f U max σ i j n j are defined as the dimensionless traction vectors, with σ i j as the Cauchy stress tensor. Here K n and K τ are the effective permeabilities in the normal and tangential direction, respectively, with the effective permeability in any orientation, φ i , given by [33]: where φ i is the angle between a given direction and the major axis of permeability. In the above jump stress conditions, β n and β τ are the normal and tangential jump coefficients, which are considered equal in the present work as done in [23,58]. Some authors consider that the minus sign of the right hand side of (13) and (14) is implicit in the stress jump coefficients, β n and β τ [23,58]. When a continuous surface traction at the interface is considered, the corresponding matching condition can be obtained from (13) and (14) by setting the values of β n and β τ to zero.
In the present paper, the RUC filling is carried out at constant inlet pressure for all cases. Therefore, the nondimensional inlet boundary conditions are as follows: At the Stokes domain: At the Brinkman domain: wherep represents the prescribed inlet pressure. No-flux condition,û ini = 0, and zero traction in the tangential direction, t iτi = 0, are used in boundaries where symmetry is specified. At the fluid fronts, kinematic and dynamic boundary conditions are defined. The former condition establishes that the fluid front advances along its normal direction, while the dynamic condition accounts for the discontinuities of normal stress due to the capillary pressure, p cap . Both kind of conditions are represented in (18) and (19), (20), respectively: • Kinematic condition (For all domains): dxi /dt =û n n i = û j . n j n i (18) • Dynamic condition (At the Stokes domain): • Dynamic condition (At the Brinkman domain): where n i , p a andû n are the outwardly oriented normal vector, air pressure and dimensionless normal velocity at the fluid front, respectively. In (18),û n is given by the Stokes velocity in the channel domain and by the pore velocity in the porous medium, which, in turn, is defined as the Brinkman velocity divided by the tow porosity, ε t [32,67]. In the case of the channel flow, the capillary pressure is p cap = −σ κ = −σ ∂n i /∂ x i , with σ as the surface tension and κ as the curvature of the moving boundary. For the flow inside the porous medium, the capillary pressure can be computed using the following equation [35]: where A int,s is the cross sectional area of the solid particles in the flow direction, C int,s is the wetted perimeter of the solid particles, θ is the contact angle and ε t is the tow porosity. In this work, the geometry of the porous medium is considered as a bank of aligned circular fibers and the capillary pressure in any point of the longitudinal tow (warp) depends on the angle between the normal of the fluid front and the principal axis of the fibers, ϕ. If the flow infiltrates the warp in parallel, ϕ = 0, or oblique direction, 0 < ϕ < π/2, A int,s and C int,s are approximated as the cross sectional area and perimeter of a truncated circular cylinder using Gauss-Kummer series for the perimeter of the resulting ellipse [5], as given by: If the impregnation occurs perpendicular to the fibers, ϕ = π/2, the following equations are used: On the other hand, for the transverse tow (weft), it is accepted that the fluid front moves inward perpendicular to the fibers and it is used the mean capillary pressure given in [38] considering an hexagonal array of fibers: where R and d are the fiber radius and half-distance between fibers, respectively. The model of Gebart [16] for a hexagonal array is used to calculate the permeabilities in the principal directions: where V f = 1 − ε t is the fiber volume fraction, while the parameters c, c 1 and V f,max are given by the following expressions:

Integral equation formulations and numerical techniques
The Stokes-Brinkman problem of this work is solved using direct BEM and DR-BEM for the channel and tows domains, respectively, with the corresponding matching and boundary conditions presented in Section 2. The integral formulation for the Stokes equation is the following: where c i j = (α/2π ) δ i j , with α as the solid angle in the source point, whose value is α = π for points located over a smooth contour. For points located inside the domain, c i j = δ i j . In (36), u j and t j are the velocities and tractions in the field points, respectively. The fundamental solutions of the integral kernels of (36) are given by: where ξ and y are the source and field points, respectively. On the other hand, for the anisotropic Brinkman equation, Khor et al [26] deduced fundamental solutions in the transformed Fourier space, which reduce to the isotropic fundamental solutions in the real space when χ 1 = χ 2 [27,62]. However, the transformation of these functions to the real space in the anisotropic case is not a trivial problem and to avoid these difficulties, a boundary-domain integral formulation in terms of the Stokes fundamental solutions is considered for the Brinkman equation and the resulting domain integral is transformed into a boundary integral using DR-BEM [42]. The integral formulation is as follows: In DR-BEM, the non-homogeneous term, g j (y) = χ 2 j u j (y), is approximated using Radial Basis Function (RBF) interpolation given by Augmented Thin Plate Splines (ATPS). The augmented part of a generalized thin plate spline of order n is a polynomial of order n−1 that is added to obtain an invertible interpolation matrix [18]. In the present work, n = 2 and the form of the ATPS is as follows: where N B is the number of boundary points, N L is the number of interior points and r (y, z m ) = |y − z m | is the distance between the field points, y, and the trial points, z m . Accordingly, the non-homogeneous term can be expanded as follows: where α m l represent the approximation coefficients in the direction l. The ATPS represented in (41) to (44) requires the addition of orthogonality conditions, as shown in the following equation: After substituting (45) into (40), the integral representation takes the following form: The transformation of the domain integral into a boundary integral is accomplished by defining the following auxiliary Stokes field: with the particular solutions for the ATPS given in [10]. The substitution of the auxiliary field defined in (48) and (49) into (47) and the application of the Green's identities in the domain integral, lead to the following boundary-only integral representation: where the coefficients α m l are given in terms of the inverse of the interpolating matrix obtained by collocation of (45) at N B boundary nodes and N L internal nodes.
The singularities arising in the kernels K i j and U j i of (36) and (50) are dealt with the rigid body motion principle [45] and the Telles transformation [61], respectively. In these equations, the boundary and the physical variables are discretized using quadratic isoparametric interpolation. Additionally, discontinuous shape functions with a collocation factor of α dis = 2/3 are employed at the corners of the problem domain [69]. After the discretization of the contour and variables and the assembly of the resultant matrices, by considering the corresponding boundary and matching conditions, a linear system of equations is obtained, which is solved using singular value decomposition (SVD) to consider cases where the condition number of the global matrix is very large due to large differences on the magnitude of the problem parameters. In the DR-BEM formulation the coordinate systems of the saturated porous domains are continuously updated as the fluid front advances, in such a way that each coordinate system is located in the corresponding centroid of each saturated domain to avoid the increment of the condition number of the final system as the filling takes place.
First order Euler integration of the kinematic condition (18) is used to track the moving fronts. As an explicit time stepping-algorithm is used, the time step needs to be restricted to small values. The Courant-Friedrich-Levy (CFL) condition is used to guarantee the stability of the solution as time progresses, with the CFL constant changing as a function of the capillary ratio, C cap . As C cap increases, smaller values of this constant are required, leading to shorter time intervals. Additional geometrical restrictions shall be considered to compute the time interval in order to avoid the crossing of the evolution points as the fluid front advances and to avoid the advancement of the points beyond the limits between the two regions considered (channels and bundles) or beyond the boundaries of the global domain. At a point x i on the moving fluid front, both the unit normal vector,n i , and the curvature, κ i , are computed numerically using the following fourth-order lagrangian polynomial [52]: Once the meshes of the moving fronts have been reconstructed at the current time step and the normal and curvatures have been computed, the BEM/DR-BEM algorithm is used to calculate the velocity of the moving fronts and the cycle is repeated again in a quasi-static approach given the low Reynolds number of the problem. A more detailed description of the tracking technique of the fluid front can be found in a work recently published [43].

Assessment of accuracy and convergence
To validate the proposed BEM/DR-BEM scheme for the solution of coupled Stokes-Brinkman problems a fully developed flow in two horizontal layers, porous medium and channel, is considered as shown in Fig. 1. This problem admits analytical solution and has been used by other authors to validate FEM codes [11,23,58]. In this problem, the inlet and outlet boundaries are subjected to pressure boundary conditions and the upper and lower ones, to symmetry boundary conditions. Both continuous and jump stress matching conditions at the interface between the two layers are considered.
Equations (55) and (56) give the analytical solution for the velocity profiles (with φ = μ/μ e f f ), which is only valid provided that the boundary layer thickness of the Brinkman flow is smaller than the height of the porous medium, in such a way that the solution tends to a Darcy flow in the lower part of the porous medium domain.
This problem was solved by [23] using a modified Brinkman approach for the whole domain that reduces to the Stokes flow in the channel domain. In this approach, the stress jump condition is incorporated by a level-set formulation , μ e f f = μ/ε t with ε t = 0.5, and K = 10 −4 cm 2 , where K is the permeability of the porous medium in the flow direction. Considering a characteristic length of L * = H for this problem, the inverse Darcian number is χ 2 = 5.00×10 3 . Four mesh-sizes were evaluated, h = 0.1 cm, h = 0.05 cm, h = 0.01 cm and h = 0.005 cm, using a uniform distribution of a regular mesh of square elements over the entire domain, with h as the length of the side of one square element, obtaining a total number of elements in each case equal to 0.1/h 2 . In a similar fashion, in the present case, h represents the size of one quadratic element of the contour mesh, which leads to  Table 1, as well as the convergence rate, n c . The L 2 relative error norm is defined as: where u x anal and u x num are the analytical and numerical solutions of the dimensionless horizontal velocity, respectively. The convergence rate, n c , corresponds to the exponent of the power curve that fits to the data of L 2 relative error norm versus Element Size, i.e., L 2 = a.h n c , where a is a constant and n c is the slope in a log-log plot of the L 2 relative error versus Element Size (see Fig. 2 with the corresponding coefficient of determination R 2 = 0.988). According to [3], the application of the Finite Element Method (FEM) to fluid flow problems implies the use of mixed formulations where multiple field variables shall be considered, like the velocity and pressure in the case of an incompressible fluid flow. In such a cases, the discretization scheme of the domain should fulfill three conditions with the purpose to assure the solvability, stability and optimality of the FEM solution, namely, consistency, ellipticity and infsup condition, being the last one the most difficult to satisfy due to the choice of numerical constants to be introduced in the formulation or to the modification of the original FEM scheme in order to satisfy implicitly such a condition. The statement of the Inf-Sup condition depends on the problem being analyzed; in the case of Stokes-Darcy problems, this condition is detailed in [36]. For both Stokes and Brinkman flows, some distorted FEM meshes could not satisfy the Inf-  Sup condition, resulting in spurious pressure modes (for more details see [3]). On the other hand, in the BEM/DR-BEM formulation for both flow fields, spurious pressure modes are not possible due to the unique relationship between the velocity and pressure fundamental solutions in the integral formulation of the problems. However, significant numerical errors and inaccuracy can be found by using very distorted BEM meshes, due to ill conditioning of the global matrix system and/or near singularities in the numerical integrations. Therefore, in our scheme, it is convenient to avoid a distorted mesh, i.e., a mesh having adjacent elements with very dissimilar sizes, to preserve the accuracy of the solution. As mentioned before, since the points at the fluid front are not uniformly spaced just after the fluid front advancement, we have implemented a remeshing algorithm with the purpose to obtain a balanced mesh in every time instant and avoid in this way the loss of accuracy in the solution.
The comparison of the present BEM/DR-BEM scheme with the FEM scheme of [23] is shown in Fig. 3a, b, where a semi-log plot is adopted in order to distinguish the velocities in the porous domain, which can be several orders of magnitude lower than the velocities in the channel. In the Brinkman domain, the positions of the interior trial points of the RBF interpolation are determined by an extension into the domain of the boundary mesh points, with exception of the points close to the corners. In both schemes, BEM/DR-BEM and FEM, the numerical solution converges to the analytical one, reaching a high order of accuracy when using the finer mesh, which corresponds to a value of h = 0.005 cm, i.e., 4.000 square finite elements uniformly distributed for the FEM scheme, and, in the case of the BEM/DR-BEM scheme, 460 boundary elements distributed over the external boundary and interface between the two domains (Stokes and Brinkman), with 484 boundary and 7761 interior interpolation points in the porous domain (Brinkman), as shown in Table 1.
Two main differences between the two results can be observed. Firstly, the BEM/DR-BEM shows high accuracy in the Stokes velocity profile, channel flow, for all mesh sizes, while the Brinkman velocity profile is over-predicted for the coarser meshes (Fig. 3a). On the other hand, the FEM scheme always predicts an accurate velocity profile in the porous medium (Brinkman), but the velocity profile in the  (Fig. 3b). As commented in [23] the observed behavior in the FEM solution is due to the interpolation scheme used in the level-set formulation. To be more specific, it is necessary to mention that a single equivalent momentum equation for the coupled Stokes-Brinkman domain was considered in [23], which is essentially a Stokes equation modified with permeability and jump stress terms to account for the fluid flow in the porous medium and the stress matching condition in the interface, respectively. Additionally, the domain geometry was defined by a level set function, Φ, where Φ = 0 at the interface, Φ > 0 at the free-fluid domain and Φ < 0 at the porous medium. When |Φ| < ε int , with ε int as the half thickness of the diffuse interfacial region, an interpolation function is defined to express the permeability and jump stress terms as a function of Φ and ε int , and the errors of such interpolation led to greater numerical errors in the Stokes velocity profile, as explained in [23]. On the other hand, the observed behavior of the numerical solution with the mesh size in the BEM/DR-BEM approach is due to the approximation error on the evaluation of the volume integral in the integral representation formula (40) used to represent the Brinkman flow. As mentioned before, in the case of the Stokes flow an exact only-boundary integral formulation is known, equation (36), requiring only the discretization of the boundary integral densities, u i and t i , without any additional approximation. The discretization of the densities, u i and t i , is also required in (40), as well as the corresponding approximation of the volume integral. The transformation of the domain integral appearing in (40) into the boundary integrals arising in (50) by DR-BEM approximation involves the interpolation of the non-homogeneous term, g j (y), using Augmented Thin Plate Splines (ATPS) as shown in (45). In this case, the approximation error of such interpolation, which is greater as the permeability is lower and/or the mesh is coarser, is the principal error source of the Brinkman velocity profiles. It is worth-mentioning that in the BEM /DR-BEM results no oscillations are present in points close to the interface Stokes-Brinkman for any mesh-size, contrary to the FEM scheme used in [23] where oscillations can be noticed for the coarser meshes (see details of Fig. 3a, b).
To analyze the influence of the non-dimensional parameters χ 2 and β on the accuracy of the solution, several simulations were performed using the data summarized in Table 2. As observed, simulations for values of χ 2 = 1 × 10 4 , χ 2 = 1 × 10 3 , χ 2 = 1 × 10 2 and χ 2 = 18 were carried out, considering in each case three different values of the jump stress coefficient, namely, β = 0, β = 0.5, and β = 1.0. In all cases, the same number of boundary elements and interpolation trial points were considered, namely, 49 boundary elements, 56 boundary and 95 interior interpolation points in the porous domain (Brinkman). In Table 2 the obtained L 2 relative error norm for each case is also reported. As observed, the accuracy of the solution improves as the value of β increases for χ 2 = 1 × 10 4 , χ 2 = 1 × 10 3 and χ 2 = 1×10 2 , however, for the case of χ 2 = 18, this behavior is reversed. Besides, for a given value of β, a higher accuracy is found as χ 2 decreases, which is reasonable because in the BEM/DR-BEM scheme the approximation error of the nonhomogeneous term of the Brinkman equation has a stronger influence on the results as χ 2 is larger.
The velocity profiles corresponding to the simulations of Table 2 and their comparison with the analytical solutions given by (55) and (56) are presented in Fig. 4a-d. The magnitude of the Stokes velocity increases as the value of β reduces, being this effect more significant as χ 2 is lower (see Fig. 4d). As expected, for a constant β, the thickness of the boundary layer increases with the reduction of χ 2 . On the other hand, for a constant value of χ 2 , the boundary layer thickness increases as the value of β reduces, being this change more notorious as χ 2 is smaller. Similar results to those ones shown in Fig. 4a-d were reported in [58] using a FEM approach, which means that the present scheme is consistent with previously reported results.

Statement of problem of void formation, simulation data and void characterization
In the processing of composites materials, different sorts of RUC architectures can be examined for dual-scale fibrous reinforcements. In this section, it is considered a 2D geometry emulating a longitudinal plane of a cross ply fabric (Fig. 5).
The unidirectional filling of 3D RUC's of fabrics has been recently considered by superposition of 2D simulations at several longitudinal planes of the RUC, showing good agreement with experimental results [7]. Therefore, we can infer that the parametric study for the 2D geometry represented in the Fig. 5 can be useful to understand the influence of some factors on the void formation in cross ply fabrics. Simulations considered in the following sections are classified into four different series. Series 1 to 3 provide different data to be used in the simulations performed with the Stokes-Brinkman formulation, while Serie 4 defines the data of the simulations run with the Stokes-Darcy formulation. The parameters of the Serie 1, which is taken as the reference case, are presented in Table 3, where the fixed parameters are prescribed and the calculated parameters are determined from the former ones. The jump stress coefficient used in the simulations of Serie 1, β = 1.24, was calculated using the model of Valdés-Parada et al. [64,65] together with the Larson-Hidgon coefficient in the channel-tow interface, γ * [29] : where the porosity of the channel-tow interface, ε ct , is approximated as the porosity of the bundle, ε t , and μ e f f = μ/ε t . In Serie 2, the jump stress coefficient is changed to β = 0.7, corresponding to the approximation considered in [23], and μ e f f = μ/ε t as before. For Serie 3, the continuous stress condition is considered, i.e. β = 0, and μ e f f = μ.  For the jump stress cases (Serie 1 and 2), considering a characteristic length of L * = H t , where H t is the tow height (see Fig. 5), the inverse Darcian numbers in the principal directions are χ 2 1 = 1.40 × 10 3 and χ 2 2 = 6.80 × 10 3 , whereas for the continuous stress case (Serie 3), the corresponding values are χ 2 1 = 2.26 × 10 3 and χ 2 2 = 1.10 × 10 4 . For simulations of Serie 4, a Stokes-Darcy formulation is considered and therefore the stress jump condition is not applicable. For this type of formulation, uniqueness of the solution requires continuity of surface tractions and a jump condition for the interfacial tangential velocities given by the Beavers-Joseph slip condition, for more details see (71) and (72). The slip coefficient, γ , of the slip condition, see (71), is usually found by mean of experiments, but when the ratio between the height of the porous medium and the square root of the permeability is very high, a good approximation for this coefficient is γ = 1/ ε t 1/2 [22]. The formulation is completed by requiring continuity of the normal component of the interfacial velocities, see (70). In Section 4.3, the results obtained using the Stokes-Darcy formulation are compared with those ones obtained with the Stokes-Brinkman formulation.
As a constant inlet pressure condition is considered in all simulations, the comparison between the pressure and capillary forces is described by the capillary ratio, C cap = P cap,max /P in , as defined in [30], where P cap,max and P in stand for the maximum capillary pressure and inlet pressure, respectively. For Series 1, 2 and 4, three different capillary ratios are contemplated: C cap = 1×10 −2 , C cap = 1×10 −1 and C cap = 5 × 10 −1 , whereas for Serie 3 only the lower value of C cap is taking into account, i.e, C cap = 1 × 10 −2 . For each simulation the following parameters are achieved: • Size, shape and location of voids, expressed as (see Fig. 5): where A void is the void area, A RU C is the RUC area, a and b are the semi-major and semi-minor axes of the ellipse circumscribing the void, a 1 is the semi-major axis of the weft and l void is the distance between the front edge of the weft and rear edge of the void. • Dimensionless time,t, defined by (6), and normalized time that is defined as t * = t/t f illing , in which t f illing is the total filling time of the RUC. In this work, two adjacent RUC's are considered and the filling stops when the partial equilibrium of the bubble formed in the second RUC has been reached. • RUC's saturation versus Normalized time • Compression of the voids, defined by the following compressibility ratio: where

Comparison between Stokes-Darcy and
Stokes-Brinkman approaches

Description of the Stokes-Darcy approach
Two coupled domains shall be considered in the mesoscopic modeling of dual-scale fibrous reinforcements, namely, channels and tows. The modeling of the flow in the porous media (tows) can be done by the Brinkman equation, which reduces to the Darcy equation when the permeability is very low.
Despite that the majority of cases of composite processing involve low permeability tows and, therefore, are prone to be modeled by the Darcy approximation, some authors prefer to use the Brinkman equation in its original form [11,23,58]. This is due to the possibility of imposing explicitly the matching conditions at the channel-tow interface in terms of the interfacial velocity and surface traction components, given that the Brinkman and Stokes partial differential equations are of the same order. As expected, several authors have used the Darcy equation to model the flow in the tows [24,43,60], but, considering that this is a first order partial differential equation, a jump matching condition on the tangential velocities at the channel-tow interface is required, which is linear proportional to the Stokes surface traction in the channel and is function of a slip coefficient, γ (Beavers-Joseph slip condition). The value of the slip coefficient, γ , appearing in this condition is still matter of controversy. Considering that both approaches, Stokes-Brinkman and Stokes-Darcy, have been used in the literature and are consistent with the problem dealt here, the present section is devoted to compare the results of void formation obtained by both of them, namely, Serie 1 for Stokes-Brinkman and Serie 4 for Stokes-Darcy. It is important to mention that a coupled solution system is considered for both approaches, i.e., the equations of the freefluid and porous medium domains, as well as the matching conditions, are directly included in a single solution system, which could be ill-conditioned as mentioned before. In other works, decoupling strategies have been used for the solu-tion of Stokes-Darcy problems, such as: iterative subdomain methods [8], Lagrange multipliers [31], two-grid method [37], among others. For instance, Mu and Zhu [36], who used the Saffman matching condition for the tangential velocities (which comes from disregarding the tangential Darcy velocity in the Beavers-Joseph condition), proposed and assessed a decoupling FEM methodology based on interface approximations via temporal extrapolation for non stationary cases. This methodology allows solving two decoupled subproblems independently by invoking conventional Stokes and Darcy solvers. After analyzing the behavior of the convergence rate and approximation errors with the time step and element size regarding a coupled strategy, it was concluded that the proposed methodology is computationally effective for this kind of problems.
In dimensionless form, the Darcy equation is as follows: are the non-dimensional coordinate and non-dimensional velocity in direction i, respectively. The integral representation formula for the pressure field when the Darcy velocity (64) satisfies mass conservation (1) is given as (for more details see [6,13,44]): (65) where: At the channel-tow interfaces, the flow fields defined by (36) and the corresponding pressure (Stokes flow), and (64) and (65) (Darcy flow), have to satisfy the following matching conditions: • Continuity of normal velocities: • Beavers-Joseph slip condition for tangential velocities: • Continuity of surface tractions: where s stands for Stokes and d stands for Darcy. Heren i is the normal vector outwardly oriented from the Stokes domain, γ is the slip coefficient, τ i is the tangential vector, K 1 and K 2 are the main permeabilities,σ (s) i j is the dimensionless Cauchy stress tensor in the Stokes domain andp (d) is the dimensionless pressure in the Darcy domain. In the Beavers-Joseph condition (71), the tangential Darcy velocity cannot be disregarded a priori and is approximated in terms of the interface Darcy pressure using Lagrange interpolation functions for quadratic elements, L j (ζ ), as follows (Sum on i and j) : At the inlet of the Darcy domain, the non-dimensional boundary condition is: while, at the fluid front, the kinematic condition is given by (18) and the dynamic condition is given by: where the parameters appearing in (74) and (75) were previously defined in Section 2. The Stokes-Darcy formulation used in this work was previously validated in [43].

Comparison of the RUC filling process for
In Fig. 6, different instants of the filling process predicted by the Stokes-Darcy (S-D) and Stokes-Brinkman (S-B) approaches are compared; in each instant, the fluid front position along the channel is the same for both approaches, S-D and S-B, but different evolution times are predicted. The fluid front shapes at the first instant look very similar for both approaches (Fig. 6a vs. b), with S-D predicting a slightly greater arrival time than S-B. As expected, due to the low value of the capillary ratio considered, C cap = 1 × 10 −2 , the fluid fronts in the channel surpass the fluid front in the tow, and, after a while, the liquid completely surrounds the first weft, see Fig. 6e, f. At the instant corresponding to Fig. 6c, d, when the liquid at the channel is approximately half the way along the first weft, 6.57% of the total filling time has elapsed according to S-D, while S-B predicts a value of 5.78%. Additionally, the minimum position of the fluid front in both the warps and the weft is barely larger for the S-D simulation. When the channel fluid fronts totally surround the first weft and merge one another, the air is trapped and the void compression starts. According to Fig. 6e, f, the fluid fronts in the warps and the weft for S-D are ahead with respect to S-B, resulting in a smaller initial void (bubble) in the case of S-D. Moreover, the real arrival time to this position is longer for S-B, but the normalized time is shorter instead. This last feature is common for all the cases analyzed in Fig. 6a-h, where the normalized times of the fluid front evolution are always shorter for the S-B simulation. Figure 6g, h corresponds to the case when the fluid front in the channel is approximately at 80% of the total length of the domain; at this point of the simulation, the time in the S-D approach is 57.2% of the total evolution time, while the time in S-B is 54.7%. This is a manifestation of the reduction of the saturation rate as the fluid front progresses with respect to the initial rate, as it is confirmed later. When the channel fluid front totally surrounds the second weft another bubble is formed, undergoing a compression process until the partial equilibrium is attained (see Fig. 6i, j), as in the first bubble. In these last two figures, it can be observed that the voids predicted by the S-D simulation are smaller than those ones predicted by S-B, with the bubble of the first weft smaller than the one of the second weft for both approaches. In general, the influence of the type of approach, S-D or S-B, on the void location is not as significant as the influence on the void size and shape, for more details see Table 4 4

.3.3 Effect of the capillary ratio on the void characterization
Let us now consider the effect of the capillary ratio on the characterization of the voids (size, shape and location) formed in the first and second weft, and the differences among these voids when they are predicted by the two formulations considered here, S-D and S-B. In our analysis three different capillary ratios are considered: C cap = 1 × 10 −2 , C cap = 1×10 −1 and C cap = 5×10 −1 ; the solutions and corresponding analysis for the former capillary ratio, C cap = 1 × 10 −2 , were already discussed in the previous sub-section (Fig. 6).
As it can be noticed from (21) to (30), the maximum capillary pressure in the present case is always achieved at the symmetric boundary of the warp, and it is given by: Consequently, in all the cases considered here, where the values of σ, θ, R and ε t are kept constants, the value of the maximum capillary pressure, P cap,max , is always the same and therefore different values of the capillary ratio, C cap , corresponds to different values of the inlet pressure given by: P in = P cap,max /C cap .
As it was observed before in Fig. 6i (S-B) and j (S-D), for C cap = 1 × 10 −2 the bubble of the first weft is smaller than the one of the second weft. In such case, this happens because the first bubble supports a higher pressure considering that it is closer to the inlet and consequently experiences a higher compression. The values of the compressibility ratios obtained with both formulations, S-D and S-B, are: for the bubble of the first weft, ψ = 0.512 (S-B) and ψ = 0.587 (S-D), whereas, for the second one, ψ = 0.304 (S-B) and ψ = 0.332 (S-D). Figure 7a-d show the results for the other two values of the capillary ratio (C cap = 1 × 10 −1 and C cap = 5×10 −1 ), where the flow in the channel has advanced enough to allow the formation of the bubbles in both wefts. For these two cases of C cap , the bubble at the first weft is larger than the one at the second weft for both formulations, S-B and S-D, which is opposite to the case described before for C cap = 1 × 10 −2 , where the second bubble was larger than the first one. This apparent unexpected behavior is due to the larger values of the capillary ratio considered in these last two cases, C cap = 1 × 10 −1 and C cap = 5 × 10 −1 , in comparison with C cap = 1 × 10 −2 . From the previous analysis of the maximum capillary pressure it follows that when P cap,max is constant, as in the present cases, larger values of the capillary ratio, C cap , correspond to lower values of the inlet pressure, P in . At lower inlet pressure, the fluid fronts at the channel move more slowly; besides, the capillary forces always promote the weft impregnation. These two simultaneous effects result in smaller differences between the location of the fluid fronts at the channel and the location of the fluid front at the wefts before the channel fluid fronts enclose the wefts, and, consequently, in smaller initial bubbles, as C cap is higher. This consequence is in turn more pronounced at the second weft due to the slower channel flow velocities obtained there regarding the ones obtained in the first weft, given the increase of the flow resistance as the fluid front progresses; accordingly, the initial void size at the second weft will be smaller than the one at the first weft. On the other hand, at a lower inlet pressure, the bubbles are subjected to lower compression, and consequently, at larger values of the capillary ratio (C cap = 1 × 10 −1 and C cap = 5 × 10 −1 ), the initial void size is very similar to the final void size for both wefts; as the initial void is smaller for the second weft, the final void is too, as can be observed in Fig. 7 for both approaches, S-B and S-D.
The results of the void characterization obtained with both formulations, S-D and S-B, are summarized in Table 4, for all values of C cap considered. As can be observed from that table, the largest voids in the first weft are obtained for C cap = 1 × 10 −1 in both formulations and the smallest ones, for C cap = 5 × 10 −1 , with the void size corresponding to C cap = 1×10 −2 in between them. On the other hand, for the second weft, the void is always smaller as C cap is higher, with no void formation, Si ze * = 0, for C cap = 5 × 10 −1 when using the S-B approach (see Fig. 7c), while S-D predicts the formation of a very small bubble (see Fig. 7d). Regarding the void shape, the increase of C cap generates bubbles with a higher aspect ratio for both wefts using both approaches, S-D and S-B.

Comparison of saturation curves
The time evolution of the saturation curves predicted by the S-D and S-B approaches for the three capillary ratios are compared in Fig. 8a, where it can be observed a similar general behavior for all curves. Figure 8b shows the saturation curve obtained with the S-D simulation for C cap = 1×10 −1 in order to describe the general behavior of all cases reported in Fig. 8a. In Fig. 8b, different filling instants are highlighted using marker points, which are numerated from 0 to 5. As can be observed, the highest saturation rate, i.e., the slope of the curve, is reached at the beginning of the injection (from 0 to 1). As soon as the fluid front arrives at the first weft, point (1), the saturation rate decreases until the merging of the channel fluid fronts (from 1 to 2) due to the flow resistance exerted by the first weft. When the channel fluid fronts have surrounded the first weft and encountered one another, point (2), an increase of the saturation rate can be noticed. This rate is kept almost constant until the arrival of the fluid front to the second weft (from 2 to 3), moment from which the saturation rate progressively decreases again until the second merging of the fluid fronts (from 3 to 4) due to the resistance exerted by the second weft. In point (4), saturation rate barely increases, and from (4) to the end of simulation (5), saturation rate is essentially constant. The nearly constant saturation rate between 2 and 3 is expected since no porous obstacles are present in the channel between these points, in contrast with what happens between the points 1 to 2 and 3 to 4, where the saturation rate decreases due to the presence of the first and second weft, respectively.
According to Fig. 8a, for C cap = 1 × 10 −2 and C cap = 1 × 10 −1 , the saturation curves predicted by both approaches (S-D and S-B) are almost identical, but the final saturation is slightly higher for the S-D result (see detail in Fig. 8a). The highest differences between the saturation curves of both approaches (S-D and S-B) are observed for C cap = 5×10 −1 , where the capillary effects are the most relevant; in this case, the final RUC saturation is barely higher for the S-B scheme (see detail in Fig. 8a).

Influence of the jump stress coefficient
Let us first consider the influence of the jump stress coefficient, β, on the void formation by comparing the results obtained from Series 1 (β = 1.24) and 2 (β = 0.70), for the  three different values of C cap considered here. The total fillings of the RUC for Serie 1 and the different values of C cap were already analyzed in the previous section and reported in Fig. 6i (C cap = 1 × 10 −2 ), 7a (C cap = 1 × 10 −1 ), and 7c (C cap = 5 × 10 −1 ). Equivalent results for Serie 2 are shown in Fig. 9a-c, where, in the case of C cap = 1 × 10 −2 (Fig. 9a), different fluid fronts with positions in the channel corresponding to the same ones reported in the first column of Fig. 6 for the Serie 1 are highlighted using arrows. As can be seen from these figures (first column of Figs. 6 and 9a), the real, dimensionless and normalized arrival times corresponding to these positions are shorter for β = 0.70 (Fig. 9a). In Table 5, the effect of C cap on the void size at the first and second wefts for both values of β (β = 1.24 and β = 0.70) is presented, with the smaller voids found for β = 0.70 in both wefts; the void size difference between both values of β is less significant for C cap = 5 × 10 −1 , i.e. the larger capillary ratio examined here. Although the size of the first bubble is not a monotonic function of C cap for both values of β, with the smaller voids corresponding to C cap = 5 × 10 −1 and the larger ones, to C cap = 1 × 10 −1 , the size of the second bubble is a decreasing monotonic function of C cap , and the resulting total void size, which is the sum of the size of both bubbles, also decreases with C cap . This last behavior is in agreement with results previously reported in the literature [30], where unidirectional macroscopic simulations were carried out using the software LIMS and it was found that the saturated length corresponding to the arrival of the  fluid front to the end of a cavity is larger as C cap increases, which, from a mesoscopic viewpoint, can be interpreted as the reduction of the total void size with the increment of C cap for a given position along the cavity. In general, the void shapes obtained in this work are physically consistent with some experimental works. For instance, Hamidi et al. [20,21] studied the void morphology in the Resin Transfer Molding process (RTM), observing different void shapes inside the bundles: • Cigar-shaped bubbles with a very high aspect ratio, as those ones reported in Figs. 7c-first weft, 7d-second weft, and 9b-second weft. • Elliptical bubbles with a lower aspect ratio than the previous ones, see for example Figs. 6i, 7a, b, 9a, b-first weft. • Almost circular bubbles, which are less frequent, having an aspect ratio close to one, see Fig. 10 According to the results in Table 6, the void aspect ratio at the first and second wefts increases with C cap for both jump stress coefficients, β = 0.70 and β = 1.24. In the case of C cap = 1 × 10 −2 , the void aspect ratio does not change considerably with β for both wefts, while for C cap = 1 × 10 −1 , a significant reduction of the aspect ratio with the increase of β is found. On the other hand, when C cap = 5 × 10 −1 , only the first bubble is formed for both values of β and the aspect ratio increases with β.
The obtained numerical results do not reveal a relevant influence of β and C cap on the location of the void, which is always formed at the rear edge of the weft. This is also in agreement with the experimental results of Hamidi et al. [21], where radial injections in a circular mold were conducted, founding that the most of voids were formed at the rear edge of the tows with respect to the flow direction, both inside the tow and in the transition tow-channel. Other numerical works, [24,25,43,54], have also predicted the void formation at similar locations.
In summary, it can be concluded that the process of void formation by mechanical entrapment of air is defined by a dynamic balance between the inlet and capillary pressures, which is characterized by the capillary ratio, C cap . The void size and shape are also influenced by the magnitude of C cap , as well as by the value of the jump stress coefficient, β. On the other hand, the void location appears to be independent on both C cap and β, since the voids are always located at the rear edge of the wefts, which is in agreement with previous numerical and experimental works.

Comparison of void size between continuous stress and jump stress simulations
The influence of the type of interface condition on the void formation for C cap = 1 × 10 −2 is studied by comparing results of Serie 2 (stress jump condition, β = 0.7) and Serie 3 (continuous stress condition, β = 0). The total RUC filling for the continuous stress condition is presented in Fig. 10, where the fluid fronts having the same positions in the channel to those ones depicted in Fig. 9a are highlighted by the corresponding arrows. The comparison between both simulations (Fig. 9a vs. Fig. 10) indicates that for the same fluid front position in the channel, the positions of the fluid fronts in the warps and wefts are more advanced in the continuous-stress case (β = 0), and consequently, the initial size of the first and second bubbles, corresponding to the instant when the chan-nel fluid fronts completely enclose the first and second weft, respectively, is smaller in the continuous stress conditions. The initial and final void size at the first weft in the case of the jump-stress condition, Fig. 9a, is Si ze * = 1.52 × 10 −2 and Si ze * = 6.57 × 10 −3 , respectively; on the other hand, these values for the continuous-stress condition, Fig. 10, are Si ze * = 6.30 × 10 −3 and Si ze * = 2.69 × 10 −3 . In the case of the jump-stress condition the compressibility ratio is ψ = 0.567, whereas for the continuous stress this value is ψ = 0.573. Accordingly, the void compression does not have a relevant influence in the difference between the final void sizes for these two cases (continuous-stress and jumpstress) since the compressibility ratios are almost the same, and the main reason of this difference is the initial void size, which is smaller in the case of β = 0, leading to a smaller final void in such a case.
A similar behavior can be appreciated for the void size at the second weft, namely, the initial void is smaller in the continuous-stress case and the compressibility ratio is almost the same in both cases (continuous-stress and jump-stress); the difference between the final void sizes at the second weft for these two cases is even higher than at the first weft, with Si ze * = 8.80×10 −3 in the case of the jump-stress condition, and Si ze * = 2.02×10 −4 for the continuous-stress condition.

Comparison of saturation curves between continuous stress and jump stress simulations
Finally, let us analyze the differences between the predicted saturation curves of both conditions, continuous stress and jump stress. Saturation curves are very important on the quantitative description of the time and space evolution of the saturated volume, and they can be very useful for the comparison of filling processes at different filling instants. The obtained saturation curves for the continuous (C-S) and jump (J-S) stress conditions are compared in Fig. 11. From this figure, it can be observed that the time evolution of the corresponding saturation curves is very similar to the one presented in Fig. 8b, with the C-S condition predicting a higher value of saturation at any normalized time, t * = t/t f illing . In the saturation curves shown in Fig. 11, the representative filling instants reported in Fig. 8b (i = 1, 2 stand for the RUC saturation obtained with the C-S and J-S simulations, respectively. It is very important to take into account that, in this case, ds is not referred to the saturation difference for a determined value of the normalized time, t * , but to the saturation difference corresponding to any of the instants highlighted in Fig. 11, which do not necessarily correspond to the same values of t * in both simulations, C-S and J-S. In the first filling instants (1), when the flow reaches the first weft, the RUC saturation is very similar for both simulations (ds is very small) because it is primarily determined by the saturation of the channel, which is almost the same in both cases. During the impregnation of the first weft (from 1 to 2), ds increases until the merging of the channel fluid fronts at (2), due to the higher weft saturation obtained in the C-S simulation during this period, as commented before. The change of the saturated volume of the first weft is not completely defined when the channel fluid fronts enclose it since the void undergoes a continuous compression during the filling process of the RUC until the partial equilibrium is reached. Therefore, the change of the saturated volume of the first weft from the time of merging of the fluid fronts to the time of arrival of the fluid front to the second weft (from 2 to 3) is due to the compression of the bubble. Given that during this period the filling of the channel and the warps is almost identical for the C-S and J-S simulations, the change of ds is essentially determined by the compression of the bubble. Accordingly, as the compressibility ratio, ψ, of C-S and J-S is very similar for the first bubble, as mentioned before, and considering the definition of ψ in (63), the change of the volume of the bubble is larger in the case of the larger initial void size, which corresponds to the J-S simulation. Therefore, the change of the saturated volume of the first weft is also larger for the J-S simulation between the time instants analyzed (from 2 to 3) and this leads to the reduction of the saturation difference, ds (see Fig. 11).
In the C-S simulation, a larger infiltration at the second weft occurs between the instant of the arrival of the fluid front to this weft and the instant of merging of the channel fluid fronts to form the second bubble (from 3 to 4). During this time interval, both simulations (C-S and J-S) predict almost the same changes on the saturated volume in the channels and the warps, and the change on the saturated volume in the first weft due to the void compression is negligible in comparison with the change of the saturated volume in the second weft. These factors lead to an increment of ds that is mainly determined by the change of the saturation in the second weft, which, in turn, is larger for the C-S condition as mentioned before (see Fig. 11). At the end of the simulation, the RUC saturation is s = 0.928 for C-S and s = 0.904 for J-S.

Conclusions
In this work, BEM techniques have been implemented for the simulation of voids formation in dual-scale fibrous reinforcements using Stokes-Darcy (S-D) and Stokes-Brinkman (S-B) mathematical formulations. The BEM solution for the S-D formulation was previously validated in [43], while a coupled Stokes-Brinkman problem with analytical solution was used here to validate the developed BEM/DR-BEM scheme (BEM for the Stokes domain and DR-BEM for the Brinkman domain), obtaining a very good accuracy. The error of the numerical solution increases with the magnitude of the inverse Darcian number, χ 2 , for all values of the jump stress coefficient, β, considered here, and decreases with the magnitude of β for each value of χ 2 , except for the smallest one, χ 2 = 18, where the error increases with β (see Table 2).
The present numerical solution predicts a reduction of the Stokes velocity in the channel and of the boundary layer thickness in the porous medium with the increase of β and χ 2 , which is in agreement with the results in [58]. Moreover, our solution is also compared with a FEM solution previously reported in [23], finding that both approaches converge satisfactorily to the analytical solution, with larger errors in the Brinkman domain than in the Stokes domain in the present BEM/DR-BEM scheme, contrarily to the FEM scheme of [23] where the errors in the Stokes domain are greater. Conversely to the FEM solution, no numerical oscillations in the boundary layer were obtained with the present BEM/DR-BEM scheme.
The problem of void formation during the filling of two adjacent RUC's, for three different capillary ratios, C cap = 1 × 10 −2 , 1 × 10 −1 , 5 × 10 −1 , was analyzed by using Stokes-Darcy (S-D) and Stokes-Brinkman (S-B) formulations, resulting in a similar behavior of the void size and shape with C cap for both formulations. At the first weft, the change of void size with C cap is not monotonous, with the largest voids obtained for C cap = 1 × 10 −1 , the smallest ones for C cap = 5 × 10 −1 , and the intermediate ones for C cap = 1 × 10 −2 . For the second weft, the void is smaller as C cap increases. The total void size, which is the sum of the size of both bubbles, also decreases with C cap , which is in agreement with the results of [30]. On the other hand, the void aspect ratio in both wefts increases with C cap .
The results of the RUC's fillings also allow concluding that the differences between the S-D and S-B formulations depend on the capillary effects. For instance, in the particular case when the capillary effects are weak, i.e., C cap = 1 × 10 −2 , the S-D formulation predicts smaller voids and lower aspect ratios than S-B, in both wefts. On the other hand, for a same value of C cap , similar saturation curves were obtained with both formulations, S-D and S-B, with exception of the case when the capillary effects are dominant, i.e., C cap = 5×10 −1 , where important differences between the saturation curves of both formulations were appreciated.
According to the present results, the type of matching conditions considered in the Stokes-Brinkman formulation affects the void size and shape. When using the jump stress condition (J-S), the smallest voids at both wefts were obtained for the lowest value of the jump stress coefficient, β = 0.7, for all values of C cap , with the exception of the case of the second weft when C cap = 5 × 10 −1 (when capillary effects are dominant), where no void formation was obtained for the two values of the jump stress coefficient (β = 1.24 and β = 0.7). Regarding the void shape, the aspect ratio always increases with C cap for each value of β in both wefts (see Table 6). According to the BEM/DR-BEM simulations, the void location is not significantly affected by C cap and β, and is formed at the rear edge of the weft.
The general behavior of the void size and aspect ratio with C cap in both wefts, which was described above, did not change with β, as well as the decreasing behavior of the total void content with C cap , which was obtained in all cases of the S-B formulation. On the other hand, for C cap = 1×10 −2 , the lowest void content was obtained for the continuous-stress condition (C-S), which is consistent with the larger saturation predicted by the C-S condition for all filling instants in comparison with the one predicted by J-S. Despite the saturation values differ in all filling instants, both saturations curves, C-S and J-S, showed a similar behavior.