Finite Strain Modelling for Multiphase Flow in Dual Scale Porous Media During Resin Infusion Process

Resin infusion is a pressure-gradient-driven composite manufacturing process in which the liquid resin is driven to flow through and fill in the void space of a porous composite preform prior to the heat treatment for resin solidification. It usually is a great challenge to design both the infusion system and the infusion process meeting the manufacturing requirements, especially for large-scale components of aircraft and wind turbine blades. Aiming at addressing the key concerns about flow fronts and air bubble entrapment, the present study proposes a modelling framework of the multiphase flow of resin and air in a dual scale porous medium, i.e. a composite preform. A finite strain formulation is discussed for the fluid–solid interaction during an infusion process. The present study bridges the gap between the microscopic observation and the macroscopic modelling by using the averaging method and first principle method, which sheds new light on the high-fidelity finite element modelling.


Introduction
In recent years vacuum assisted resin infusion (RI or VARI) becomes an important method for manufacturing the large-scale composite components of structures like airplane wing covers and wind turbine blades [1,2]. A resin infusion process mainly includes two steps, (1) liquid resin infusion through a dry composite preform and (2) curing (i.e. heat treatment) to transfer the liquid resin into a solid matrix for holding fibres together. From the point of view of mathematical modelling, the former may be considered as a fluid transport process inside a porous medium consisting of a fibre network, the latter is a chemo-thermo-mechanically coupled process with phase changes and the creation of residual stress. The present study is restricted to the step 1, liquid resin infusion process. Therefore, chemical and thermal effects are not discussed explicitly. The schematic illustration of a resin infusion system and a representative element of the composite preform accommodating fluids are shown in Fig. 1. The preform shows a dual scale microstructure, i.e. fibres and bundles of fibres, which may experience finite strain deformation during a manufacturing process. For non-crimp unidirectional (UD) fabrics, bundles could be stitched together to form plies and a composite preform. This preform is placed such that one side is sitting on the mould and the other side being covered and sealed by a vacuum bag. Inlets under the atmosphere pressure and outlets under the vacuum pressure are connected to either the mould or the vacuum bag for the liquid resin entering or leaving the porous preform driven by the pressure gradient. The aim of liquid resin infusion is to fully fill in the void space of the preform by the resin. In other words, an optimal infusion process should maximise the wetting of a preform with a minimal amount of resin [2]. Thus the design of inlet and outlet system (e.g. locations and the number of outlets/inlets), the assisting flow media to divert flow paths, the time control of opening and closing of inlets and outlets, and the lay-up of composite plies are the key factors to be considered for a successful manufacturing. The challenges are (i) the component scale dry domains formed due to flows taking the least resistance pathes and (ii) the bundle and fibre scale air bubbles trapped behind flow fronts. The latter, if cannot be eliminated, will end up with defects in a cured component. This is a potential source of material damage like delamination. Therefore, the key questions are to understand the mechanisms of the formation of dry domains and air bubbles and to address the challenges in the design of infusion system and infusion process. Bearing in mind that this is a fluid-solid interaction process, a fluid's motion induces changes of the microstructure of a solid skeleton and in turn affects the flow itself.
Darcy's law has been widely discussed for modelling resin flows in composite manufacturing. The calculation or computation of permeability is one of the key tasks of modelling. Due to the dual scale microstructure of a preform, the liquid resin flow in the void space of a preform is a dual scale flow containing the intra-and inter-bundle flows. Typical length scales of such two types of flow are, < 10 µm and > 100 µm, respectively [3]. According to the theory of permeability of porous media [4][5][6][7], it is clear that the intra-bundle void space has the permeability very different from that of the inter-bundle void space in both magnitude and anisotropy.
Generally, the analytical formulae of Kozeny-Carman [8], Gebart [9] or Berdichevsky [10] may be used to predict the permeability as a function of the volumetric fraction. Based on those formulae, a typical analytical method to estimate the permeability of a dual scale composite preform is to integrate the estimations at different scales together. For example, Gebart [9] proposed a theoretical estimation of the permeability of an idealised unidirectional reinforcement consisting of the regularly ordered, parallel fibres derived from the first principles (i.e. Navier-Stokes equations) both for the flow along and for the flow perpendicular to the fibres. Using Gebart's intrabundle result and the additional estimation of permeability of the inter-bundle channels, Lundström [3] suggested the formulation of permeability of dual scale preforms by integrating the intra-and inter-bundle solutions together. Such an analytical estimation usually considers the fluid as a single phase. Thus air bubbles are omitted. Even for a single phase flow, the dual scale microstructure is still highly complicated so that numerical methods have been discussed extensively for predicting the permeability [11][12][13] and experimental methods have been developed for parameter measurement and model validation [14][15][16][17].
The multiphase flow adds more complication to a dual scale system. Considering its relatively lower density and volumetric fraction by comparison to the resin, the air phase can be modelled in various ways subjected to the assumptions applied. The existing models may roughly be classified into two categories: (1) the void space model without explicitly modelling the air phase motion and (2) the air phase motion model. The void space model naturally refers to the concepts of the saturated and unsaturated permeabilities [18,19,21]. A complicated multiphase flow can be characterised by a phenomenological model of the relative permeabilities with volumetric fractions as the state variables. Such a model generally comes within the scope of Darcy's law which in principle is limited to a homogeneous and random porous medium [4].
By contrast, the explicit modelling of the air phase motion may provide the micromechanical insight into the mechanisms of void formation and motion. Park et al. [22] and Gangloff et al. [23] reported the detailed hypotheses about the void formation by air entrapment and the air bubble motion aiming at the manufacturing strategy to minimise the void phase. The direct computational fluid dynamics (CFD) modelling may provide the enhanced visualisation for improving the understanding of air-resin interaction [24]. In the circumstance that the air phase is taken into account, mathematical modelling has shown explicitly that Darcy's law may not cover some key features observed in the dual scale multiphase resin infusion system. Pillai and coworkers [25][26][27] used the microscopic model and averaging method to demonstrate that the flow flux equation is in a form of Brinkman's equation, which is consistent with the model proposed by Celle et al. [28].
As the resin and air are two immiscible fluids, it is interesting to understand if capillary forces play important role in the resin infusion. In the well-known Buckley-Leverett theory (see e.g. [29]), the capillary force is ignored so that two immiscible fluids have the same pressure. However, Amico and Lekakou [30,31] conducted the experimental study and came to the conclusion that the effects of the capillary forces cannot be ignored when the resin pressure is sufficiently low. More precisely, a dimensionless number, either capillary number, Ca * , or bubble mobility, may be used to characterise the competition between viscous flow and capillary wicking, which may be useful to depict the local flows of resin and air and be important for the removal of voids [22,23].
Bearing in mind that the resin transports through a solid skeleton composed of fibres and fibre bundles [28,32], without being glued by a matrix, i.e. the cured resin, fibres and fibre bundles may experience the elastic finite strain deformation during an infusion process. Larsson et al. [34,35] proposed a model taking the fibre bundle (i.e. the tow including the fibres and intra-tow voids) as a hyperelastic material. Similar to Pillai's model [25], a bundle is considered as a 'sink' to the mass balance of resin in the inter-bundle void space. It is noted that the finite strain deformation is one of the features in the resin infusion process which may not be considered in other dual scale media, e.g. fractured media [33].
In addition to the aforementioned material characterisation and modelling of dual scale porous media, finite element analysis (FEA) has been reported for the component scale modelling of resin flow processes [20,21,36]. Commercial software for both material modelling (e.g. Digimat ® ) and process modelling (e.g. ESI-Composites ® ) have achieved success to a certain extent. However, to the author's best knowledge, especially that of the modelling of industry manufacturing processes of wing covers and wind turbine blades, it is still a great challenge to create a model feasible to cope with the simulation of a manufacturing process of large-scale components with the reasonably useful detail of air bubbles. This is a typical challenge of the time-dependent cross-scale modelling where a component could span up to more than ten metres while the size of individual air bubble may be down to the same scale of fibre diameter. The transient nature of the problem and the pressure-dependency of permeability coupled to the finite strain deformation of the preform during an infusion process make the modelling even more challenging.
In order to address the challenge, the present study proposes a framework of finite strain dual scale fibre-reinforced model explicitly taking into account the multiphase flow of the resin and air. The finite strain modelling of the fibre reinforced soft materials has achieved substantial progress in the last decade [37]. And the image processing and data analysis technologies provided by the advanced imaging software, e.g. Thermo Fisher Avizo ® or Materialise ® , offer the support to construct the more sophisticated microstructure-based constitutive model. Moreover, the increasing data from lab experiments and industry manufacturing practice have not been fully exploited by the mathematical and finite element modelling, let alone predicted by the models. Therefore, it is interesting and important (1) to exploit the advanced finite strain fibre-reinforced model and imaging-based microstructure model for the more sophisticated modelling of resin infusion processes and (2) to construct a high-fidelity modelling framework which is able to accommodate the high-precision data of a dual scale flow and in turn predict the high-accuracy detail of interest. For this purpose, the present study adapts the dual porosity concept of fractured media [5] and extended it to a finite strain model of dual scale fibre-reinforced materials. It is assumed that the same fluid (either the resin or the air) in the intra-and inter-bundle void spaces may be considered as two different phases. Such two phases can interchange mass and momentum on their interface. Therefore, there are five phases, i.e. the resin in the intra-bundle void space, the resin in the inter-bundle void space, the air in the intra-bundle void space, the air in the inter-bundle void space, and the solid phase of a preform. Hence, the proposed method may capture some key features of interest, e.g. local flow fronts and defects induced by air bubbles at the different scales. It is also worth emphasising that the proposed framework of finite strain modelling can be extended consistently to cover the curing process by adding the chemo-thermal equations.
The present work is arranged as follows: First, the finite strain kinematical theory is presented to define the dual scale multiphase flow system. Secondly the fundamental balance laws of the individual phase at the microscopic scale are discussed for understanding the motion of each phase and the constitutive law. Thirdly, the macroscopic governing equations of mass balance and momentum balance are obtained by averaging the corresponding microscopic counterparts. The constitutive laws of macroscopic stresses taking into account the capillary forces are also obtained by using the averaging rules. Finally, instead of deducing the macroscopic constitutive laws of fluid fluxes from the macroscopically thermodynamic constraint, the present work discusses the macroscopic fluxes using the averaged momentum balance equations, which would be helpful for clarifying the multiphase flow in dual scale media and shedding new light on the finite element modelling of the manufacturing process of large-scale composite components.

Kinematics theory
We consider an REV in Fig. 1 composed of a deformable solid skeleton and fluids flowing through the void space. Mathematically, the kinematical relations of the different configurations and particles during a motion are shown in Fig. 2.
At a certain initial time, t 0 , the solid skeleton is represented as a continuum by its reference configuration, B 0 , with material points labelled by X s ∈ B 0 . At any current time, t ∈ [0, +∞), a mapping χ s : B 0 → B t is a deformation which maps the reference configuration, B 0 , onto the current configuration, B t embedded into the Euclidean space R 3 . The current position of a material point X s is written as x s = χ s (X s , t) ∈ B t ⊂ R 3 . The deformation gradient is defined as F s (X s , t) = ∂χ s /∂ X s with its Jacobian denoted as J s = det F s > 0.
The solid skeleton may be characterised by a dual scale porosity fibre-reinforced microstructure, i.e. the fibre bundle package in which each bundle is consisting of one family of fibres (see Fig. 1). Let a unit vector, a 0 (X s ), denote the initial fibre direction of the family of fibres attached to a material point, X s , of a solid particle on B 0 .  The current spatial fibre direction at the same material point is where the deformed fibre and its stretch are described by a = F s (X s , t)a 0 (X s ) and λ a = C s : A 0 = tr A, (2) in which ':' denotes the double contraction, 'tr' is the trace operator, the right Cauchy-Green tensor is defined as and A 0 = a 0 ⊗ a 0 and A = a ⊗ a are two structural tensors. For studying the fluid-solid interaction, it is convenient to consider a spatial REV defined by a representative domain of the solid skeleton on B t and the associated current void space. Porosity, i.e. the volumetric fraction in the REV, as the first order correlation function [38] is used to characterise the void space of a dual scale microstructure. Accordingly, the REV can be characterised by two independent porosities, e.g. the porosity of the total void space, φ, and that of the intra-bundle void space, φ f . The porosity of inter-bundle void space is denoted The subscripts b and f stand for the inter-bundle void space surrounding the bundles and the intra-bundle void space surrounding the fibres, respectively.
The void space is initially fully occupied by the air close to the vacuum pressure. Driven by the pressure gradient created by the pressure differences between inlets and outlets ( Fig. 1), the liquid resin may gradually enter and occupy part or all of the void space as illustrated in Fig. 3. Therefore, the resin and air are considered as two immiscible fluids and denoted as wand n-fluids, respectively. In a conventional way, the w-fluid stands for the wetting fluid while n-fluid for the non-wetting one. Together, the wand n-fluids completely occupy the void space in an REV. For effectively modelling the multiphase flow in the dual scale void space, we refer to a fluid in the intra-bundle void space and to the same fluid occupying the inter-bundle void space as two 'apparent phases' [5].  Thus, there are five phases in total, denoted s, wb, w f , nb, and n f , referring to (i) the solid phase, (ii) the w-fluid between the bundles in the inter-bundle void space, (iii) the w-fluid between the fibres in the intra-bundle void space, (iv) the n-fluid between the bundles in the inter-bundle void space, and (v) the n-fluid between the fibres in the intra-bundle void space, respectively (see Fig. 3 for the geometric relations and Fig. 4 for the topological relations). Hence, we have the relations where θ α represents the current volumetric fraction of the α-phase (α ∈ {s, wb, w f, nb, n f }) in an REV. The total volumetric fractions of the wand n-fluids read respectively. Based on the above characterisation of the void space, we can define the kinematics of the fluids in an REV by pulling an α-fluid particle located at x α ∈ R 3 in the current void space at time t back to its position, X α ∈ R 3 , at the initial time t 0 (Fig. 2). Therefore, the kinematics of the solid and fluids may be defined in a single framework. Let X α ∈ R 3 (α ∈ {s, wb, w f, nb, n f }) be a material particle of the α-phase as a continuum. At the current time t, the location of a X α -particle at a point in an REV may be expressed by a spatial coordinate vector, x α = χ α (X α , t) ∈ R 3 via a mapping χ α : R 3 → R 3 . Noting that χ α (X α , t) = χ β (X β , t) (α = β) holds only if two particles, X α and X β , overlap on the current interface between the αand the β-phases.
The material velocity of the α-phase is defined as where˙≡ D/Dt is the material derivative. The spatial velocity gradient of the α-phase is where ∇ is the spatial gradient operator. It can be proven that L s =Ḟ s F −1 s for the solid phase.

Microscopic balance equations
Let's consider a current moving material domain of the α-phase, α (X α , t), mapped from a reference material domain, 0 α (X α )(≡ α (X α , t 0 )) by χ α : R 3 → R 3 . The surface of α (X α , t), denoted α (X α , t), is a material surface with respect to the α-phase moving in a velocity V α . There is no mass flux of the α-phase passing through the α -surface. The intrinsic density of the α-phase (per unit volume of α ) is denoted ρ α which generally is a function of pressure and temperature.

Mass balance of a phase
In the case that there is no mass source (growth) or sink (absorption) in α (X α , t), a mass balance simply reads D Dt The Reynolds transport theorem [5,39] states that D Dt in which ν α is the outer normal vector of the α -surface. Thus, by using Gaussian theorem, Eq. (8) becomes D Dt is an arbitrary domain, the local format of the mass balance of an α-phase is obtained from Eq. (10) as

Linear momentum balance of a phase
The linear momentum balance of the moving material domain of an α-phase, α (X α , t), bounded by the surface, where σ α and f α are the Cauchy stress tensor and body force vector applied on the α-phase, respectively. Following the similar process to obtain Eq. (11), the local format of the balance of linear momentum is obtained as where ρV α ⊗ V α and σ α are interpreted as the rates of momentum gained by the advection and diffusive momentum transfers, respectively, f α represents the rate of supply of momentum by a body force, all terms are per unit volume of the α-phase. By using the mass balance (11), Eq. (13) can be re-expressed in terms of the velocity V α , known as equation of motion, as or in an equivalent form

Energy balance of a phase
The density of the total energy of the α-phase, denoted e α , per unit volume of α (X α , t), is the sum of its thermodynamic internal energy part and kinetic energy part, where I α is the specific internal energy (internal energy per unit mass), * = √ * · * is the Euclidean norm of a vector.
The energy balance of the α-phase states that D Dt where the various terms on the right hand side (r.h.s) of Eq. (17) contributing to the energy balance are: (i) the advective energy supplied through the α -surface in the first integral on the r.h.s, (ii) the work of forces acting on the phase by the surface force, ν α · σ α , in the second integral, (iii) the work of forces acting on the phase by the body force, f α , per unit mass in the third integral, (iv) the diffusive energy supplied through the α -surface by a heat flux, j H α , in the fourth integral, and (v) the rate of production of energy within α by a heat source, H α , in the last integral.
The local format of Eq. (17) can be worked out as By using the mass balance (11) and linear momentum balance (14), Eq. (18) can be re-expressed in a concise form For obtaining the constitutive model feasible to FEA, the Helmholtz free energy of solid phase, ψ s , and Gibbs free energy of an α-fluid, ϕ α , are introduced respectively as follows: where T is temperature, s α (α ∈ {s, wb, w f, nb, n f }) is the specific entropy (per unit mass) of the α-phase. Given the α-fluid's pressure, p α , defined on the current fluid particle at x α , is the pull-back of p α from x α to X α in the way as transferring the Cauchy stress to the second Piola-Kirchhoff stress. By substituting Eq. (20) 1 into Eq. (19), the energy balance of the s-phase becomes Similarly, by substituting Eq. (20) 2 into Eq. (19) and using the relationships, the energy balance of an α-fluid (α ∈ {wb, w f, nb, n f }) is obtained as Whereby we use the stress decomposition of an α-fluid, which decomposes the Cauchy stress σ α into the viscous stress τ α and pressure p α .

Entropy balance equation
Equations (22) and (24) are the expressions of the changes of specific entropies derived from the energy balances.
On the other hand, the general form of an entropy balance can be expressed as where j S α and S α are the diffusive entropy flux and rate of production of entropy per unit mass of an α-phase, respectively. The second law of thermodynamics states that the dissipation of an α-phase, denoted D α , must satisfy the constraint, By comparison of Eq. (26) with Eq. (22), it can be identified that for the solid phase the expressions of j S α and S s are and respectively. Similarly, for an α-fluid we may compare Eq. (26) with Eq. (24) and identify that and

Constitutive model of a solid
In order to further explore the thermodynamics constraints upon the constitutive laws, the free energy functions ψ s and ϕ α are specified in the forms as By substituting Eq. (32) 1 into Eq. (29) and applying the inequality (27) (α = s), the constitutive constraint upon the solid phase is obtained as where D s is the rate of deformation of the solid skeleton defined by The inequality (33) indicates the constitutive laws and constraint, The above inequality in (35)  Bearing in mind that the solid skeleton has a dual scale microstructure, i.e. fibres and bundles, it may be more convenient to introduce the constitutive law of the fibre bundle as one family fibres at the microstructure scale and to estimate the macroscopic property over an REV by the homogenisation of a package of the bundles. Thus, the free energy functionψ s is proposed to have the form as [40] whereC s = J −2/3 s C s , and the invariants are defined as As an example, a particular neo-Hookean free energy function [40] is adopted herê where κ and μ are interpreted as the bulk modulus and shear modulus, respectively, of the fibre bundle and κ 1 and κ 2 are the parameters of the elastic response of the fibres. The explicit expressions of the stresses and elastic tangent modulus derived from the free energy (38) may be referred to [41].

Constitutive model of fluids
In the similar way, by substituting Eq. (32) 2 into Eq. (31) and applying the inequality (27) Each term on the left-hand side (l.h.s) of inequality (39) is discussed individually. First, for the arbitraryṖ α we obtain Secondly, let's consider the Newtonian fluid with a constitutive law, where C denotes the fluid's viscosity coefficient and D α = 1 2 (L T α + L α ) is the rate of deformation of the α-fluid. In this case, the constraint (39) requires that which means the positive-semi-definite condition of C . Hence, assuming that the molecular structure of the fluids (resin and air) are statistically isotropic, the coefficient C is proposed in a form of an isotropic tensor, where a and a are two scalars, I is the fourth-order unit tensor with components I i jkl = δ ik δ jl + δ il δ jk , in which δ i j is Kronecker delta. By substituting Eq. (43) into (41), the constitutive law of a compressible isotropic Newtonian fluid (say, air) reads and where μ α = a is the α-fluid's dynamic viscosity or shear coefficient viscosity while λ α + 2 3 μ α is called the dilatational viscosity or coefficient of bulk viscosity in which λ α = a [5]. Substituting either (43) or (44) into the constraint (42) yields 2 ) and (tr D α ) 2 , are independent from each other and both non-negative.
For an incompressible isotropic Newtonian fluid (e.g. resin), Eq. (44) reduces to and the mean normal stress, is computed by using the incompressibility condition, ln J α = 0 or tr D α = 0, while p α itself may be treated as a Lagrangian multiplier [42].
Thirdly, similar to (35) of the solid counterpart, inequality (39) indicates which also leads to the Fourier's law of heat flux,

Macroscopic model of the multiphase flow in a dual porosity medium
As indicated by the specification of free energy functions in (32), a macroscopic model may take the averaged values of the solid's motion χ s and fluids' pressures p α (∀α ∈ {wb, w f, nb, n f }) over an REV as the set of macroscopic degree of freedom (DoF). Accordingly, from the point of view of the variational principle and finite element method, the governing equations for computing a macroscopic system may include: (i) the macroscopic mass balance of each fluid phase, (ii) the macroscopic linear momentum balance of the solid skeleton, (iii) the macroscopic constitutive law of the effective stress acting upon the solid skeleton, and (iv) the macroscopic constitutive laws of the massor volume-weighted fluxes of the fluids. One may observe that the macroscopic mass balance, macroscopic linear momentum balance and macroscopic constitutive laws of stresses can be obtained by averaging their microscopic counterparts over an REV, which is shown in this section. However, the macroscopic constitutive law of a fluid flux, e.g. Darcy's law, does not have an apparent microscopic counterpart. Instead, as an application of the first principle method, it can be derived from the averaged linear momentum balance equation of the fluid [5,25].

Averaging rules
Let's consider a spatial point, x ∈ B t ⊂ R 3 , and its associated REV, denoted U 0 (x) ⊂ R 3 , which is fixed in the Euclidean space. The subdomain inside the REV occupied by an α-phase is denoted as U 0α (x, t) ⊂ U 0 ⊂ R 3 , which is time-dependent.
In order to differentiate a macroscopic model from a microscopic one, a local coordinate system, x , is assigned to the REV for the microscopic description (see Fig. 1). This x coordinate system is locally isomorphic to the global coordinate system, x, for the macroscopic description. Therefore, the configuration of an α-phase in the REV is Given a microscopic tensorial field, G(x , t; x), its averages over U 0α and over U 0 are defined, respectively, as and where is the characteristic function of the α-phase. The volumetric fraction of the α-phase, θ α , is the ratio of the volume of the α-phase, U 0α (x, t), to the volume of the REV, U 0 (x), Consequently, we have the decomposition, whereG is the deviation of G from G α .
Some useful averaging rules are summarised here (please refer to [5] for the details of derivation). The averaging rules of sum and of product of two microscopic tensorial fields, G 1 and G 2 , state that and The averaging rules of time derivative and of spatial derivative, are expressed respectively as and where ∇ x and ∇ x are the gradient operators in the x and x coordinate systems, respectively, and we adopt the symbol [5], to represent the surface average over an interface, S αβ(α) (x, t), in the REV. The set, symbolically represents the union of all of the neighbouring phases of an α-phase. Hence S αβ(α) (x, t) is the the union of all of the interfaces between U 0α and any other phase, U 0ζ (ζ = α), moving with a local velocity u. S αβ(α) is the area of S αβ(α) (x, t) whileS αβ(α) is the specific area (per unit volume of the REV), defined as

Symbolic expressions of the topological relations between phases
Bearing in mind that β(α) refers to the set of all of the neighbouring phases surrounding an α-phase in the REV, this naturally requires the definition of the topological relations of the different phases as shown in Fig. 4. In order to effectively characterise the topology of a dual scale multiphase system, we introduce the symbolic phase matrices, α and β, to represent the phases and their complete set of neighbours, respectively. We write and where the components, β The interface between the α i j -phase and its complete set of neighbours, β i j , is denoted as S α i j β i j . Hence we have Noting that there is no summation over the subscripts i and j.
It is proposed that β f i j may be expressed as a function of the α's components in such a form where e i jk is the alternating symbol [42] (or the Levi-Civita symbol). Equivalently an explicit matrix form of β f is which apperantly can be decomposed into where It is clear that either ∪ i, j∈1,2 S α i j β i j or ∪ i, j∈1,2 S α i j β i j is a complete set of the fluid-fluid interfaces in the REV but with opposite normal directions. This explains the meaning of the negative sign in β f .
The matrix β can be further decomposed into where Hence the set, ∪ i, j∈1,2 S α i j β c i j , is the complete set of the interfaces between two immiscible fluids (w and n). On the other hand, S α i j β h i j represents an interface between two phases which are the same fluid (either w or n) occupying the intra-and inter-bundle void spaces, separately. In other words, S α i j β c i j is an interface between two materials while S α i j β h i j an interface between two scales.
The other form of decomposition of β f is where are the miscible and immiscible (unmixable) neighbours of α, respectively. This complete the discussion about β f in Eq. (62). The remaining term on the r.h.s of (62), β f s , is written as which does not contribute to mass balance since any fluid-solid surface is a material surface with respect to fluid mass given that there is neither absorption nor dissolution. It is useful to introduce two matrices, The former represent the set of α i j -phase's identical fluid (either w or n) occupying the intra-and inter-bundle void spaces separately while the latter the set of α i j -phase's neighbouring immiscible fluid phase and solid phase. Therefore, S α i jαi j is not a material surface while S α i jβi j is such a one. In order to establish the connection between two index systems, i.e. α and α, through out the present work the identity, holds wherever referring to either α or α i j . By comparison to β(α), β and its various decompositions provide a more convenient and clear description about the topological relations.

Macroscopic mass balance
By using the averaging rules, the averaging of the microscopic mass balance equation (11) yields the macroscopic mass balance equation of an α-phase, Using the proposed phase matrices, Eq. (73) may be re-expressed in a physically equivalent form as where all material surfaces satisfying ν α · (V α − u) = 0 are omitted from the surface integral term. The remaining , represents the mass exchange of an identical fluid (resin or air) between the intra-and inter-bundle void spaces. Thus it may be interpreted as a 'sink' on one side of the S α i j β m i j -interface and a 'source' on the opposite side.

Macroscopic linear momentum balance
Averaging the microscopic linear momentum balance equation (13) with the help of the averaging rules, we obtain To combine with the mass balance equation (73), Eq. (75) becomes Wherebyρ αV α α and(ρ α V α ) ⊗V α α are the mass and momentum dispersive fluxes, respectively. As one may be aware that(ρ α V α ) ⊗V α α is in a form similar to the so-called Reynolds stress in the Reynolds-averaged Navier-Stokes (RANS) equations of turbulence theory, the dispersive flux need the additional constitutive law and other equations for the closure of model. As shown by Bear and Bachmat [5], a non-advective flux, including the diffusive and dispersive fluxes, may be grouped with an advective flux into a single Darcy's law in a phenomenological model. Therefore, the present study focuses on advective fluxes and leaves dispersive fluxes to be discussed elsewhere.

Macroscopic constitutive model of stresses in a dual porosity medium
Equation (77) indicates the momentum exchanges between the different phases, which is one of the features of a mixture system. Therefore, rather than simply averaging the constitutive laws of each and every phase, the interactions between phases need to be considered for the derivation of a macroscopic constitutive model. This is demonstrated in the derivation of the effective stress of solid skeleton and that of the fluids' macroscopic fluxes, respectively.

Total stress
Summing Eq. (77) over all α-phases yields (bearing in mind that α ≡ α i j ) where the symbol, [ * ] α i j ,β i j , denotes the jump of ' * ' across S α i j β i j from the α i j -phase to the β i j one, the total stress is defined as a volumetric averaged stress as and the total body force per unit volume of a porous medium is Whereby we use (i) the fact that S α i j β c i j ∪ S α i j β f s i j is a material surface, (ii) the continuity of density and velocity on , which requests to differentiate the effective stress from the total stress; (ii) the interface between the immiscible fluids, S α i j β c i j , which leads to the concept of capillary pressure; and (iii) the interface between the intra-and inter-bundle void spaces, S α i j β h i j , which contributes to the coupling effects between two scales. One may observe that the last term on the r.h.s of Eq. (78) plays a similar role as that of a body force in the total momentum balance. It may be interpreted as the rate of supply of momentum by the difference in the macroscopic velocity, V , between the (same) fluids in the intra-bundle and in the inter-bundle void spaces. This is consistent to some other researchers' postulation of a 'sink' term (see e.g. [25]) in the momentum balance of an inter-bundle fluid where the intra-bundle void space plays the role of a sink through the absorption of resin from the inter-bundle void space. However, in the present study there is no need to identify if it is a source or a sink as the formulation allows either direction of momentum exchange. And it is worth noting that such momentum exchange is governed by the mass balance so this term does not request an additional constitutive model.

Capillary pressure
The momentum exchange on the S α i j β c i j -interface between the immiscible fluids is shown in the first surface integral on the r.h.s of Eq. (78). The continuity of momentum transfer [5] on such a material surface reads where γ αβ c is the surface tension between the α and β c fluids, r * defined by 1/r * = (1/2)((1/r ) + (1/r )) is the mean radius of curvature of the S α i j β c i j -interface given two principal radii of curvature, r and r , and t = ∇ x γ αβ c /|∇ x γ αβ c | is a tangential vector on the interface. Equation (81) represents the force balance on the S α i j β c i j -interface which may be decomposed into the force balances in the normal direction ν α i j and tangent direction t, respectively, as and If the viscous stress τ and force balance in tangential direction are negligible, Eq. (81) (and (82) as well) reduces to the so-called Laplace formula which defines the microscopic capillary pressure, Hence, the difference in radius between the air bubble in the (intra-bundle) f -void space and that in the (inter-bundle) b-void space indicates the difference in capillary pressure. Accordingly, the macroscopic capillary pressure is defined as which can contribute to the formulation of the effective stress of solid skeleton as shown below. To be consistent with the microscopic model (84), the phenomenological model of the macroscopic capillary pressure generally is a function of the fluids' volumetric fractions and surface tension [4,5], which can be measured by experimental methods.

Effective stress of a solid skeleton
By using the stress decomposition (25) for a fluid, the total stress may be expressed as Assuming that the viscous stress τ α i j is negligible, Eq. (87) reduces to where the effective stress [5,43,44] is and the averaged pressures of all fluids as a whole over the REV and over the whole void space in the REV are defined, respectively, as The subscript 'v' refers to the whole void space in the REV. The macroscopic capillary forces defined in (85) may be explicitly introduced into the expression of p v v as where the macroscopic capillary pressures in the b-and f -void spaces, respectively, are and the averaged pressure of the w-fluid (i.e. resin) and averaged capillary pressure of the whole void space are obtained as Given the constitutive models of the macroscopic capillary pressures, p c b and p c f (or p c v ), as the functions of the fluids' volumetric fractions [4,5,45], Eq. (91) indicates that one may choose the macroscopic pressure of the w-fluid (resin), p w w , as the DoF in modelling. This is suitable for FEA and experimental tests as the resin's pressure at the locations of interest can be measured by using pressure gauges during an infusion process.

Macroscopic stresses of individual phases and homogenisation methods
With the capillary pressure, effective stresses, and cross-scale momentum source/sink terms defined, the momentum exchanges in a dual scale multiphase system are clarified. The remaining work to complete the macroscopic constitutive model of stresses is the homogenisation of each phase's stresses. The macroscopic stresses of each phase may be obtained from averaging the microscopic constitutive laws of the solid, (35), and that of a fluid, (41), For a compressible fluid, averaging Eq. (40) yields There are different homogenisation methods in literature to obtain the averaged constitutive laws. Roughly those methods were classified into three categories, the Voigt method (assuming an uniform strain in an REV), the Reuss method (assuming an uniform stress in an REV) and the energy equivalency method, e.g. the self-consistent method [38]. In FEA, one may only need to consider the constitutive laws at the discrete Gaussian points of each element [39]. In other words, the averaging equation (49) is computed at the discrete Gaussian point, say x = x Gaussian . If χ s and p α (α ∈ {wb, w f, nb, n f }) are taken as the DoF of the finite element system, the Voigt method supposes that the gradients of the DoF are uniform in the REV of each (macroscopic) Gaussian point, x Gaussian . This means that strain measure F s is assumed uniform in the solid phase while p α varies monotonously within its domain in the void space of the REV, i.e.
Such assumed fields for homogenisation are consistent with the finite element's convergence criterion, e.g. the patch test [39].

Macroscopic fluxes in dual porosity media
As aforementioned, the macroscopic constitutive law of the mass-or volume-weighted flux of a fluid phase in a porous medium may be derived directly from the thermodynamical constraint, objectiveness, and material symmetries of the specific material at the macroscopic scale [44]. By contrast, taking the averaged momentum balance as the first principle, the micromechanical method directly estimates the macroscopic flux in an REV and consequently provides the constitutive law of the fluid's flux, e.g. Darcy's law or its enhanced formats, and the corresponding material property of the porous medium, e.g. permeability [3,9]. Aiming at understanding the complicated resin infusion process and shedding new lights on a high-fidelity FEA, the present study adapts and extends the micromechanical method proposed in [5] to the finite strain dual scale multiphase system discussed above.
In the view of averaged spatial derivative (58), stress decomposition (25) and the specification of body force as gravity, f = −g∇ x z , Eq. (77) can be re-written as Noting that the l.h.s of above equation is in terms of macroscopic variables defined at x while the r.h.s in terms of microscopic ones defined in an REV. Equation (98) is taken as the first principle for estimating the macroscopic fluxes at x. For this purpose, first, each term on the r.h.s need to be expressed in terms of the macroscopic variables, e.g. V α α and p α α , to transfer (98) into a fully macroscopic equation. Secondly, a reduced format of (98) is used to obtained the first order formulation of the macroscopic fluxes in a closed form which is consistent to Darcy's law. Thirdly, the enhanced formulations, e.g. the quadratic formulation and Brinkman-type formulation, are also presented for further exploring Eq. (98).

Averaging the pressure gradient and gravity terms in Eq. (98)
The pressure gradient term on the r.h.s of (98) reads where is the complete surface of U 0α . By using the homogenisation assumption (97), it can be proven that (see Appendix A) where the tortuosity of the α-fluid is As the resin generally moves slow enough in the preform (except the vicinity of inlet points perhaps) under the proper control in a well-designed infusion system, the resin flow is always at a very small Reynolds number (Re 1) [23,28]. Thus, one may assume that in the vicinity of S α i jβi j , the normal components of both the inertial force and the viscous resistance to the flow are negligible by comparison to the normal component of pressure gradient and gravity, i.e.
as indicated by Eq. (15). Thus, Eq. (15) projected in ν α direction reduces to which means the microscopic pressure gradient is proportional to the body force in this circumstance. Hence, by additionally assumingρ α ρ α α , i.e. ρ α ≈ ρ α α = const in the REV, the surface integral on the r.h.s of Eq. (100) Using the identity, Eq. (104) becomes Substituting (106) into Eq. (100) yields which completes the estimation of the pressure gradient term on the r.h.s of (98). Combining (107) with the gravity term also on the r.h.s of (98), the two terms are expressed in a concise form, which is a function of the macroscopic variable, p α α , and its macroscopic gradient.

Averaging the viscous resistance term in Eq. (98)
The third term on the r.h.s of (98) represents the viscous resistance to the flow which naturally may be expressed as a function of the macroscopic velocity, V α α , by using a constitutive law of the viscous stress τ α .
As an example, the constitutive law (44) of Newtonian fluid is used here. If the gradients (and deviations) of the material properties, μ α and λ α , in the REV can be ignored, Eq. (44) and the averaging rule (55) indicate that Applying the averaging rules (57) and (58) upon the r.h.s of (110), we obtain where q α = θ α V α α is the volume-weighted flux (or specific discharge [5]) of the α-phase. By adopting an additional no-slip condition in the tangential direction on the fluid-solid surface, Under this assumption and using the relation (72) 2 , the first surface integral on the r.h.s of (111) becomes Noting that specifying G = 1 for Eq. (57) yields ∇ x θ α = − ν α αβ(α)S αβ(α) , where β(α) equals toβ so that By using Eq. (114) and an approximation, Consequently, the first two terms on the r.h.s of Eq. (111) are expressed together as where q r α = θ α (V α α − V s s ) is the macroscopic volume-weighted flux of the α-phase relative to the motion of the solid skeleton.
In a similar way, the fourth and fifth terms on the r.h.s of Eq. (111) are estimated together as With Eqs. (116) and (117) obtained, the remaining work to complete the estimation of the viscous resistance is the estimation of two surface integrals in the third and sixth terms respectively on the r.h.s of (111). For estimating the third term on the r.h.s of (111), we may calculate that Noting that ≥ 0 so that x − ν α ∈ U 0α (x, t). By assuming Replacing V s (x ), V α (x − ν α ), and V α (x ) in above equation by V s s , V α α , and V α α i j β u i j , respectively, and introducing a characteristic length, α , to replace , we complete the estimation of the third term on the r.h.s of (111) as where C α and C α are the macroscopic dimensionless shape factors associated with the S α i j β f s i j respectively, on the α side of the surfaces.

Averaged momentum balance equation in terms of macroscopic variables
Governed by the last term in mass balance equation (74), the last term on the r.h.s of Eq. (98) represents the momentum exchange between fluids in the intra-and inter-bundle void spaces. This term may be taken as a nonlinear contribution of V α to (98) and simply estimated as where u ≈ V s s on S α i j β m i j since this interface actually is defined by the solid skeleton (see Fig. 3). Finally, replacing V α α by 1 θ α q α on the l.h.s of Eq. (98) and substituting (109), (125), (126) into Eq. (98), we obtain the averaged momentum balance equation of the α-fluid in terms of the macroscopic variables as

First order formulation of macroscopic fluxes
We consider the special case in which (i) the effect of inertia and (ii) the effect of internal friction in the fluid can be neglected in comparison with the viscous resistance force at the fluid-solid and fluid-fluid interfaces, and (iii) the nonlinear term of V α is negligible as well. In such circumstance, Eq. (127) reduces to Our aim is to estimate the macroscopic fluxes in a closed form using the above reduced format of momentum balance equation. For this purpose, the surface integral of microscopic velocity, V α α i j β u i j , in Eq. (128) need to be expressed as a linear function of the macroscopic fluxes.

Averaged fluid velocity on the S α i j β u i j -surface
The expression of V α α i j β u i j depends on the force balance on the S α i j β u i j -interface between two immiscible fluids. Noting that the force balance equation (81) on the S α i j β c i j -interface can be equivalently re-formulated to the S α i j β u i jinterface as which shows that the tangential force balance equation on a S α i j β u i j -interface has the same form as (83). Thus, by using the constitutive law (44) and kinematical relation, By assuming that the dominant fluid velocity variation in the close vicinity of the S α i j β u i j -interface is the variation of the tangential velocity along the normal direction ν α , Eq. (130) reduces to Taking average of Eq. (131) over the S α i j β u i j -interface within the REV, we come to an estimation, Proceeding in a way similar to that used in Eqs. (118)-(120), Eq. (132) yields Adopting the additional no-slip condition on the material surface S α i j β u i j leads to Therefore, we can work out V α

Closed form solution of fluxes
Inserting (135) into Eq. (128), we obtain the reduced format of momentum balance equation of the α-fluid as a linear equation of the macroscopic fluxes, where in which By interchanging α with β u i j , we have the other equation of the phase pair, α − β u i j , where and in which Noting that we use a transfer, α I J = |β u i j |, where, given the indexes i, j, it can be proven that there exist I, J ∈ {1, 2} satisfy such a transfer.
Solving the set of linear Eqs. (136) and (140), we obtain the macroscopic fluxes of the phase pair, α − β u i j , where By specifying α = wb, w f in (144), we can calculate two pairs, q r wb q rnb and q r w f q rn f , respectively, which completes the estimation of the macroscopic fluxes of the dual scale multiphase flow in a closed form. It is interesting to observe the driving forces of the macroscopic flux of an α-phase shown in Eq. (144). First, the α-phase's own macroscopic pressure gradient plays a role as described by Darcy's law. Secondly, there is an explicit coupling between two neighbouring immiscible phases. The forces due to pressure gradient and gravity in one fluid cause the motion in the other one, e.g. b β u contributing to q r α , due to the momentum exchange at their interface. Thirdly, the inhomogeneity in the surface tension, ∇ x γ αβ u α i j β u i j , acts as an additional driving force. Finally but implicitly, the interaction between the intra-and inter-bundle flows which is represented in the mass balance equation (74) and tortuosity T * α , is computed on the S α i jαi j -interface.

Quadratic formulation of macroscopic fluxes
The other reduced format of momentum balance derived from Eq. (127) is By comparison to Eq. (128) which leads to the closed form solution (144), Eq. (146) takes into account a nonlinear term of the velocities which explicitly represents the interaction between the intra-and inter-bundle flows. We use two approximations, = −θ α ∇ x p α α + p α α ge 3 · T * α .
By using the expression (136), Eq. (157) becomes It is difficult to obtain a closed form solution of the above Brinkman-type equation for an REV domain. On the other hand, from the point of view of FEA, one may take Eq. (158) as a point-wise constitutive equation to be computed at the Gaussian points of each element. We shall discuss a suitable numerical method for computing (158) elsewhere.

Conclusions
Vacuum assisted Resin Infusion is an important manufacturing method for composite materials. However, it is a great challenge to design the infusion system and infusion process based on an effective prediction of the resin flow in a dual scale fibre-reinforced composite preform. At the component scale, it may need the costly and timeconsuming small scale and full scale trial and error tests to mitigate the dry domains due to the effects of flow paths. At material scale, the air may flow and be trapped at the intra-and inter-bundle void spaces, which results in the defects after curing and may need the expensive nondestructive defect detection to evaluate the risk. Moreover, the fluid-solid interaction, which is a transient process due to the change of the resin pressure during an infusion process, makes the material property, e.g. permeability, process-dependent and thus adds further difficulty to the process modelling.
One of the essential challenges in finite element modelling of RI is the constitutive model capable to capture and feasible to predict the above key features of the resin infusion. Macroscopic constitutive modelling based on thermodynamic constraints and material symmetry, which generally leads to a phenomenological Darcy's law, may lack the high-fidelity to accurately depict the phenomena in dual porosity fibre-reinforced media and accordingly is not able to provide the effective prediction to support the design of infusion system and process, especially for large-scale composite components.
By taking the same fluid in the intra-and inter-bundle void spaces as the different phases, the present study propose a mathematical framework aiming at effectively clarifying and predicting the dual scale multiphase flow phenomena in the resin infusion. First, the fundamental physics laws of each phase are proposed in the finite strain format to fully understand the system at microstructure level. The microscopic mass balance and linear momentum balance equations play the role to understand the motion of phases while the energy balance and entropy balance lead to the constitutive laws of each phase. Secondly, a symbolic system, denoted phase matrices, is proposed to clarify the topological relationship between different phases. Using this symbolic system, the averaged mass balance and linear momentum balance equations are obtained as the governing equations of the macroscopic modelling. And the macroscopic finite strain constitutive laws of stresses can be obtained by the homogenisation and development of the microscopic counterparts. This is important to address the fluid-solid interaction. Finally, the fluid fluxes are directly derived from the averaged momentum balance equations, which is aiming to clarify the mechanisms of the multiphase flow in a dual scale medium and sheds new light on a constitutive model suitable for FEA. Analytical analysis and finite element modelling based on the present study will be discussed in future work for further validation and application.
In brief, the present study serves to clarify the microscopic mechanisms and to bridge the gap between microscopic observation and macroscopic modelling of the resin infusion process. The proposed method potentially may be helpful to the research in other areas involving multiphase flows in dual scale fibre-reinforced porous media.