Holistic modeling of composites manufacturing using poromechanics

researchers, cf. the review article by Hubert and Poursartip 13 , and it is well established that a true processing operation will result in a two-scale flow, cf. 14 Then many studies were carried out considering the coupling effect between the flow motion and the fiber bed deformation, cf. 4,9–11,14–18 In this contribution, the aim is to develop a formulation, which is physically beyond the formulation of the two models presented in Refs. 19,20, where they are developed based on the ideas described in Refs. 6,11.


Introduction
Fiber-reinforced composite materials are used extensively, in particular in the aerospace industry and the use of composites is assumed to grow significantly in automotive industry in near future.There are also societal issues that have been raised for CO 2 emissions reduction, where use of lightweight structures made by fiber-reinforced composites could help to achieve this goal.In these perspectives, the development of new manufacturing technologies has been the focus of the designers and manufacturers in recent years.Because of that, modeling and simulation of manufacturing of advanced composite materials has attracted a lot of attention in order to maintain control of and to reduce cycle times as well as the number of trial and error cycles.Specifically, in automotive industry, high volume production places higher demands to predict the cycle times.To this end, development of simulation tools to model and predict the fundamental manufacturing phenomena happening during a process is necessary.Besides this, manufacturing simulation is a key to predict the final product properties.
Process modeling of composite materials, which during the manufacturing may be considered as a fluid-filled porous material, can be described, on the macro-scale as well as micro-scale, within the framework of the well-founded theory of porous media (TPM), as described in the review by de Boer 1 .][8][9][10] Some further advancements of TPM, cf. 11,12which deals with formulation of the mass balance in terms of logarithmic compaction strain measures, which is used to specify the compressibility of the different phases via the entropy inequality, has been used in this contribution.In addition, the analysis of the resin flow through a porous medium has been investigated by many researchers, cf. the review article by Hubert and Poursartip 13 , and it is well established that a true processing operation will result in a two-scale flow, cf. 1414][15][16][17][18] In this contribution, the aim is to develop a formulation, which is physically beyond the formulation of the two models presented in Refs.19,20, where they are developed based on the ideas described in Refs.6,11.The intention is primarily to develop this simulation tool for research purposes since there is, to the author's knowledge, no commercial software capable of handling such a complex model.The idea of this simulation tool is to give the user the option to choose the boundary conditions, initial conditions, and material data depending on the process he has in mind to simulate, whereas the core of the method is to solve the time-dependent dual-scale compressible resin infiltration problem through the fiber preform assumed as a compressible volumetrically deformable porous material.Different deformation-dependent permeability models are also in place to choose from, depending on the type of the process and also the type of the preform (isotropic or anisotropic).
To describe the proposed objectives, the paper is outlined as follows: in section (2), the details of the homogenized TPM are given.We introduce the assumptions and the micro-constituents of the wet-out problem, where governing equations are formulated with respect to them.Furthermore, constitutive relations are derived from a thermodynamical representation of the problem.In addition, the Darcy law is presented where the deformation-dependent permeability models are introduced.In section (3), finite-element discretization is presented, the weak form and FE equations are derived and the solution procedure is explained.Then, in section (4) a number of different manufacturing methods are explained, simulated and results are presented.The results indicate that the method we have implemented numerically is capable of predicting the physical process behavior.The challenges that were addressed earlier are also represented during the simulation.Our experience with the method is that the formulation is very sensitive to the chosen material parameters.In addition, the time stepping and the geometry discretization affect the constitutive relations at different scales and in different ways.In section (5), the overall formulation and its features are discussed and finally some concluding remarks are presented in section (6).In this fashion, our formulation for simulating the composite manufacturing processes, which is based on the general continuum mechanical framework, cover a variety of processes such as liquid resin infusion (LRI), out of autoclave, and wet-press-forming prepregs.In fact, the different processes are identified in the continuum mechanical framework from the set of boundary conditions defining one and each of the processing methods.

A homogenized TPM General
In this section, we will reiterate on the biphasic TPM which at macro-level contains a mixture of solid and fluid phases.Based on this macro-scale continuum theory, we are describing the behavior at meso-scale consisting of the resin flow motion, resin flow front, preform deformation, and micro-infiltration through the fiber plies.The fiber bed during any manufacturing process is being considered as a (solid) porous fiber network.The pores are partially filled with resin and contain resin filled and unfilled void space.Figure 1 shows a schematic view of a two-phase continuum body that undergoes finite deformation.The resin flow through the preform is the single most important event happening during a process and the preform deformation is considered as coupled to the resin flow motion.This flow is of dual-scale nature, where we have to consider the micro-infiltration through the plies in addition to the macroscopic resin flow.During a process, two coupled flows may be present: inter-and intra-ply (bundle) flows.The inter-ply flow is the flow through the wide channels between the fiber tows, cf. Figure 2, whereas the flow between the fibers inside a fiber tow is the intra-ply flow, cf. 21,22Figure 2 shows the schematic flow motion in both scales.

Micro-constituents
The micro-constituents occupying the domain, with volume V, are considered homogenized.Below we have the assumptions and the respective notation for each media.
• Compressible solid particles, p, with the volume fraction ϕ p = V p /V, representing the fibers.• Incompressible liquid constituent, l, with the volume fraction ϕ l = V l /V, representing the resin.• Voids, v, with volume fraction ϕ v = V v /V, embedded in the fiber plies.• Compressible gas constituent, g, with volume fraction ϕ g = V g /V, occupying the voids.To ensure that each representative elementary volume is occupied by fluid and solid phases, we have the saturation constraint as where n s and n f are the macroscopic volume fractions for the solid and fluid phase, respectively.In turn, the solid phase is considered subdivided into particle and void constituents. (1) Likewise, fluid phase is also subdivided into liquid and gas constituents.

The solid phase and its subconstituents of fiber particles and micro-voids
We can formulate the fiber content as a function of micro-saturation degree ξ s , defined as the degree of wet out within a representative fiber ply, ξ s = (ϕ p − ϕ pd )/ϕ p , and densities ρ p and ρ s leading to the formulation of the compaction strains.In that sense, we can divide the partially saturated fibers, namely particles, into a wet portion and a dry portion ϕ pd .Figure 3 represents fiber content volume fraction ϕ induced by fluid pressure at a certain level of saturation.
Assuming a typical vacuum-assisted process leads to n s ρ s = ϕ p ρ p , where we come to It may be noted that having the initial condition at ξ s → 0 leading to 0 = s 0 ∕ p will also lead to p → s 0 ∕ 0 , which in combination with equation (3) yields the compaction of the solid phase as a function of irreversible wetting factor ξ s and the reversible fiber content volume fraction ϕ as Taking the logarithm of equation ( 4) yields the additive decomposition of the total solid compaction strain ɛ s as where the reversible compaction strain ɛ se and the irreversible wetting compaction ɛ sp are presented, both as a function of the fiber content volume fraction ϕ and the micro-saturation ξ s .In particular, the stationarity of the solid phase leads to , whereby the solid volume fraction, in view of equation (5) where the solid compaction strain ɛ s is defined, is governed by the relation Now by introducing the logarithmic compressibility strain εs = − ρs ∕ s and εf = − ρf ∕ f , and combining equations (8)   and (10) we get the local format of total mass balance of our binary mixture as It should be noted that equation (1) and its rate are used to work out ̇ns and ̇nf .
In order to assess the liquid saturation degree ξ f , we would refer you to Larsson et al. 19 for extensive details.However, we formulate the liquid content as M l = Jϕ l ρ l , and keeping in mind that ρl = 0. Assuming that the liquid velocity is the fluid velocity meaning that l = f we represent the liquid mass balance as where we also express this relation as the balance of liquid volume fraction By combining this equation with equations ( 1) and ( 7) and the fact that ϕ l = n f ξ f then we obtain the governing equation for the liquid saturation evolution as where it can be noted that the sources governing the rate of the liquid saturation degree ξf are (i) the volumetric solid deformation, (ii) the liquid divergence term ∇ ⋅ d , (iii) the fluid volume fraction n f and also (iv) the solid compaction strain ɛ s .In Section Flow front tracking problem, we use this relation as a governing equation to establish flow front tracking feature of this contribution.

Momentum balance
Considering the localized format of the momentum balance in the spatial format, we obtain the balance relations for quasi-static behavior of the mixture where = s + f is the total Cauchy stress.In turn, is related to the effective stress and the fluid pressure p via the Terzaghi effective stress principle as = − p1.

Balance of energy and entropy
Similar to the binary mass and momentum balance relationships, a binary balance of energy yields the change in the internal energy of the mixture material.Likewise, the entropy inequality for the mixture material is obtained for the isothermal behavior which combined with the energy equation lead us to (11) where J = det( ) and = ∕ is the deformation gradi- ent, cf. Figure 1.

The fluid phase considered as a gas-filled liquid
To represent the fluid content expressed in equation ( 1), the fluid mixture is formulated by defining ϕ l = n f ξ f leading to the evolution of the liquid volume fraction as where ξ f is defined here as the liquid degree of saturation on the macro-level.Further assessment of the liquid saturation evolution ξ f leads us to the initial condition that ξ f = 0 will correspond to completely gas-filled pores as we have

Governing equations
In order to simplify our model here, we further assume that the void motion is affine with the particle (solid) motion; and also the fluid content, which is considered generally compressible, and it can be reduced to one phase using the relation ϕ l = n f ξ f .So the micro-level problem of four actual constituents in the volume of material is reduced to a two-phase continuum problem.

Balance of mass applied to solid and fluid phases
Due to the mass balance of the solid phase in the continuum, we have Ṁ s = 0 where, M s = Jn s ρ s is the solid content.In par- ticular, the stationarity of the M s using equations ( 5) and ( 6) leads to with s being the solid-phase velocity.
The respective mass balance of the fluid phase, where we have M f = Jn f ρ f , is non-stationary with respect to the solid phase due to the relative motion between the phases.It can also be noticed that the compressible fluid density ρ f is written as using a linear mixture relationship and the definition of the liquid degree of saturation, whereby the fluid density is compressible due to gas density g p and the evolution of the liquid saturation degree.
The mass balance of the fluid (liquid and gas) phase has the form where d is the Darcian velocity defining the transport of liquid phase within the partially saturated continuum.We also obtain Ṁ f based on the compressibility of the ρ f and combine it with equation ( 9) to get ( 6) Neo-Hookean like material using isochoric-volumetric split written as where ̂ = J −2∕3 t ⋅ .The parameters involved in this expres- sion are the packing exponent β and the k s E factor.We remark that, by having equation ( 21), the elastic small shear deformation properties of the fiber bed (represented by the elastic shear modulus G) are generalized into the finite deformation regime.As to the response due to volumetric deformation, it is assumed that this response is represented by the same packing law as is used for the fiber plies.In this case, the argument of packing is the fiber volume fraction p = p 0 ∕J where p 0 is the initial volume fraction of particles.It should also be noted that the assumption of isotropy of the fiber bed response means that the fiber plies possess no preferred orientation and may be justified within certain isotropy planes of loading of the fiber bed.

Solid compaction and micro-infiltration
It is concluded that the solid-phase compaction consists of an irreversible wetting process corresponding to the exclusion of voids in the micro-constituents and reversible component, basically due to elastic packing also induced by fluid pressure as shown in Figure 3.As to elastic fiber packing, we note that uniaxial compression tests of fiber beds have been studied quite extensively in the literature, e.g.Kim et al. 23 In the present paper, the semi-empirical elastic fiber packing law proposed in Toll 24 is directly generalized to the compressive response of the non-saturated region of a fiber bed consisting of voids and dry particles as where p is the excess fluid pressure, relative to p 0 , acting on a representative fiber ply, as shown in Figure 3.Moreover, p 0 = kE 0 is the configurational fluid pressure in the ini- tially non-saturated representative fiber ply.The homogenized free energy is obtained from the energy equivalence n s ρ s ψ s = (ϕ pd + ϕ v ) of the representative volume leading to where the first equality is specified in material arguments ψ s (ϕ, ξ), whereas in the last equality a change in arguments has been made using the continuum compaction strain ɛ s and the elastic wetting strain ɛ sp , which is at hand for the numerical implementation.Moreover, we also introduced the ratio γ = ϕ/ϕ 0 formulated in arguments γ(ɛ s , ɛ sp ) as Straightforward application of equation ( 20) 2 yields the fluid pressures p (relative to the configurational pressure p 0 ) as p = p 0 (γ α − 1).( 21) where the different terms D s , D nvf , and D i represent dissipation contributions due to solid deformation, fluid compressibility and Darcian interaction, respectively.The details regarding how equation ( 16) is formulated in terms of free energy ψ can be found in Refs.19,20.

Constitutive equations
In this section, the constitutive equations with respect to the formulated mass and momentum balance in the content of the twophase mixture will, in the view of equation ( 16), be formulated.

Compressible fluid-filled fiber network response
Using the dissipation inequality in equation ( 16), the proper stress response of the compressible (partly) fluid-filled fiber network may now be formulated.To this end, we formulate the dissipation produced by the solid phase for the entire component as where the last expression is obtained after pull-back to the reference configuration B 0 , as dV is the infinitesimal volume element in B 0 , while dv is representing infinitesimal volume element in the current configuration B. It may be noted that the integrand JD s in equation ( 17) is rewritten in terms of the Kirchhoff stress = J and spatial velocity gradient l or the second Piola-Kirchhoff stress S and the right Cauchy-Green deformation tensor C and also the intrinsic fluid pressure p as where, in particular, the last equality is obtained from mass conservation due to the stationarity of M s = M s 0 = n s s 0 .If we assume, for simplicity, hyper-elasticity of the effective stress response we obtain the dependence ψ s ( , ,  p ) in the free energy for the solid phase.In view of equation (18), the reduced dissipation is obtained as where we assume an additive split of the free energy into macroscopic and mesoscopic parts.In view of equation ( 18), the final conclusion in equation ( 19) corresponds to the state equations where = − J −1 p is the effective second Piola-Kirchhoff stress due to the Terzaghi effective stress principle.In the following, we subdivide the stress response in terms of the effective fiber bed response and the solid-phase compaction associated with micro-infiltration induced by fluid pressure.

Effective fiber bed response
The effective response of the fiber bed is considered hyper-elastic governed by the stored energy function for ( 17) is the permeability factor due to shape and amount of void space at fully saturated conditions.We consider the permeability factor k in terms of the so-called Kozeny-Carman equation; cf.Gebart 25 , which is formulated as where r s is the particle (or fiber bundle) radius, C ≈ 10 is the Kozeny constant and ϕ is the fiber volume fraction.
Anisotropic permeability.A special model is developed by Rouhi et al. 20 to account for (i) the macroscopic flow between the layers and (ii) the infiltration flow into the dry (mesoscopic) fiber plies, cf. Figure 2, as where K mes is the permeability through the fibers, K Ch is the permeability through the channel, ϕ l is liquid volume fraction and = ⊗ is the structure tensor related to the director field t, cf.Figures 2 and 4. The permeability through the channel is approximated considering the resistance to viscous flow within a rectangular channel, cf. Figure 5 where || is the continuum stretch in the direction of t, H 0 is the half thickness of the initial prepreg, h s 0 is the initial thickness of fiber bed and h f 0 is the half of the initial thickness of the channel between the plies.

Compressible liquid-gas response
In order to assess the pressure dependence in the fluid density ρ f = ξ f ρ l + (1 − ξ f )ρ g , it is assumed that the same pressure prevails in the liquid and gas constituents and that the highly compressible gas constituent is pressure dependent in the spirit of the ideal gas law, i.e. ρ g = k g p, where typically the gas-compliance k g is determined by k g = m g /Rθ, where m g is the molecular mass of the gas, R is the universal gas constant, and θ is the absolute temperature.It should be noted that the rate behavior of the fluid density is characterized in terms of the compression modulus of the liquid-gas mixture defined as Indeed, the value of K f increases for increased saturation and decreased gas-compliance k g .For continued saturation toward ξ f → 1, we obtain that K f → ∞ and ξ f = 1 leading to fluid incompressibility, i.e. ρ f → ρ l .

Flow front tracking problem
In those applications of the model, where we are dealing with a moving surface, one may generally distinguish one In view of equation ( 22) we obtain, in compliance form, the compaction as a function of the wetting and the fluid pressure given as Now, considering the rate form of equation (25) we formally obtain where it was used that the wetting is a diffusive process represented by the viscous evolution εsp = −p∕.The parameter μ represents the viscous resistance for penetration of liquid into the bulk fibers, which is defined as where ν is the fluid viscosity and ζ is the wetting length.The mesoscopic permeability K mes for fiber plies, which is the permeability through the fiber bed, is calculated based on the Gebart equation 25 for hexagonal fiber packing where r is the fiber radius.Here, we also made the assumption that the fiber content ϕ is constant and in turn ϕ 0 was used in Gebart equation.

Darcian solid-fluid interactions
The last contribution D i ≥ 0 of equation ( 16) represents dissipation induced by drag interaction between the phases.To accommodate this dissipation, the effective drag force f e (or hydraulic gradient with negative sign) is chosen to ensure positive dissipation via Darcy's law where mac is the permeability tensor.In order to derive the permeability for the considered preform, we have different options which are tested and used in this routine.
Isotropic permeability model.Based on the development in 26 , we formulate macroscopic permeability as a constant isotropic positive definite tensor mac = k where we have where h is the distance between layers and p 0 is the initial particle volume fraction.

Deformation-dependent
Kozeny-Carman model.Once more assuming mac as positive definite, we restrict to isotropic permeability condition, mac = k where k ( where d is represented by equation ( 28).Whereby the saturation degree variable is regarded as a local field variable f = f [ , t] or, simply, as an internal variable governed by (34).

Numerical solution procedure
In order to numerically solve the established framework, the finite-element representation of the involved primary fields and governing equations of momentum and mass balances are considered.To be able to formulate the respective weak form of the governing equations we identify the primary variables for the coupled set of equations in (11) and ( 15) where a finite-element subdivision of the region D into elements D e , e = 1,..., NEL is made.It is assumed that each element has the interpolation (34) fluid-saturated portion B 0 → B 0 [t] and one non-saturated one-phase portion C 0 → C 0 [t] separated by the free surface boundary Γ I [t], as shown in Figure 6, in order to formulate the governing equations of the two different continua.However, as exploited in the foregoing subsections, one of the unique ideas of this contribution is to consider the motion of Γ I [t] in terms of the evolution of the fluid saturation field Clearly, at the initiation of a wetting process the initial condition is that f [ , t] = 0 in order to define the one-phase non- wet the fully saturated region is defined by .Let us next identify and formulate the coupled problem between liquid saturation, momentum and mass balance pertinent to our mixture continuum.As to liquid saturation, we recall the relation ( 14) along with the Darcian flow model (28), which is to be satisfied locally according to where p denotes the prescribed pressure on the external boundary.

Application
In order to numerically show the applicability of the formulation for wide range of different processes, we will run a few examples.The simulated processes are representing setups of LRI and compression molding.The choice of the manufacturing methods to study will dictate the type of boundary conditions and material parameters to be used for numerical implementation.

Liquid resin infusion
Here, we have considered the planar infusion problem shown in Figure 7.The applied pressure load due to the vacuum is prescribed along all the external boundaries, except at the impermeable lower boundary.The solid deformation is also controlled via prescribed displacement, in particular, along the lower, impermeable boundary in the figure.The wetting process is driven by the Darcian solid-fluid interaction force as induced by the pressure gradient.This process is manifested by a migrating resin flow front represented by the saturation degree evolution with time in the porous fiber bed model.Typically, the migrating front of fully saturated integration points is developing quite fast in the beginning of the wetting process, whereas later on in the infusion process the front speed reduces due to the diminishing pressure gradient at the front.In Figure 8, the deformation of the preform is shown along with the distribution of saturation degree and pressure.Preform deformation has a direct influence on the permeability, porosity and the Darcian liquid flow advancement.Figure 9 shows the global saturation degree versus time, where it is also compared with the cases that are studied in 19 , where the micro-infiltration and solid-phase compressibility were not considered.As we can see, the saturation of the domain takes longer time to be completed in the current study since the resin infiltration through the fiber plies play a role as a sink for the resin flow.Figure 10 shows an example of a hat-stringer simulated also with vacuum infusion.The flow front and the preform deformation are captured by the routine even at the corners of the mold.It is clear also that the fluid pressure is compensating some part of the deformation as the flow advances.

Compression molding
In a press-forming process, we apply an elevated load, which is ramped up until certain compaction level and then it is relaxed.It affects the process time from being a matter of minutes to seconds.For this type of processes we consider a compressive relaxation test applied to a fluid-filled fiber network, which is related to press forming of vacuum-bag-only (VBO) prepregs.The analyzed planar specimen is shown in Figure 11 along with the applied loading consisting of a prescribed vertical displacement r on the top boundary and fixed displacements on the lower boundary of the plate.As to the fluid-phase boundary conditions, two different types are considered.The first case corresponds to macroscopically undrained conditions, where all outer boundaries are considered impermeable where ̂ e u contains the element shape functions and φe is the vector element nodal placements.We then obtain the spatial velocity gradient as where N I [ ] are the element interpolation functions and φ I are the corresponding element nodal placements.Likewise, the fluid pressure field is approximated, i.e. it is assumed that each element has the interpolation where e.g.̂ e is the vector nodal pressures.We thus obtain the pressure gradient as ∇p = p ̂ p where p is the consequent pressure gradient interpolation matrix.The rate of pressure is relative to the solid reference configuration and is defined by whereby the integrated nodal displacements ̂ and nodal pressures ̂ is assessed as Then, the weak form is given as where it was used that J = ∇ ⋅ .Here in turn P is the function space containing the virtual displacement field [ ], and S is the function space containing the virtual fluid pressure field [ ].We also introduce prescribed nominal traction ̄ represent a system of fiber plies of carbon fiber and epoxy resin and a complete list of material data can be found in Ref. 20.

Discussion
The aim of developing this framework is to show a general poromechanical model, solved by finite-element method, with constitutive equations in different scales, is capable of handling complex, coupled and interrelated phenomena happening simultaneously during composite manufacturing.
Modeling challenges such as highly deformable preform and infiltration of the resin undergoing large elastic and inelastic deformations involving compressible as well as incompressible formulations are developed.The coupling effect between the flow and deformation of the preform is a source of complexity and is considered, which is mostly treated by others as decoupled.
where there is no way that the resin can be drained during the process.The resin will be dispersed through the fiber plies and channels between the fiber tows.The second case is the partly drained condition, corresponding to drainage along the vertical boundaries of the specimen, where the fluid pressure is prescribed to zero in order to simulate the drainage of the excess resin.The prescribed displacement is thus considered as a controlled displacement r on the top boundary, which is ramped based on a chosen time step up to the total compression r/H = 0.2, where H is the height of the specimen.Figure 12 shows an example of the simulation where the initial loading rate is ̇r = −0.5 m/s.
The considered relaxation test for a VBO prepreg considered by setting ̇r = 0 (i.e.r is the rate of displacement) after the specimen has been compressed 20%.The results are shown in Figure 13, where on the basis of the above loading program we monitor the wet out in the specimen depending on the rate of initial loading ̇r and relaxation time for different loading rates for the undrained case.In addition, the permeability model discussed in equation ( 31) is used to model the macroscopic Darcian permeability.The routine is tested using other permeability methods that are elaborated in the Section Darcian solidfluid interactions as well.In a previous contribution, results with respect to constant and deformation-dependent permeability are represented, cf. 19The involved material parameters industrial tool for simulation of manufacturing processes.The technical details are in place, every corner of the routine is known, the physics are understood and implementation is complete.The possibility of extension of the model to 3D is also possible since most of the constitutive models are already implemented and simplified for a 2D problem.One can introduce new 3D elements and apply it in the routine.However, a higher order element is needed and necessary for 3D simulations.In addition, the proposed methodology, in terms of different choice of processes, needs verification and validation against experimental tests.At the same time, calibration of the material data with experimental response is also possible and will help the stability and accuracy of the simulation to increase. 27

Conclusion
In this contribution, we have developed and implemented a generic simulation routine to be used for manufacturing modeling of composite materials.This work is the extension, improvement, optimization, and combination of our previous developments in the same area, cf. 19,20We have understood the essence of such a simulation tool for manufacturing methods and put effort to develop such a routine to be able to cover simulation of all these methods under the same umbrella.
To the author's knowledge, there is no available general poromechanical formulation for composites processing methods.In that sense, we have developed a special porous media formulation in order to model a dual-scale coupled flow-deformation process with constitutive relations concerning four different mechanisms governing all the processes involving (1) constitutive effective stress response of the fiber bed, (2) infiltration of resin into fiber plies, (3) elastic packing response of the plies, and (4) Darcy's law governing the macroscopic interaction between the two phases.In addition, if necessary, the model accounts for the non-saturated behavior typical for the transition region at the flow front between full and non-saturation.
Flow front tracking in combination with a deformable fiber bed during the infusion is in place which saves all the time and efforts from using external methods such as level set, etc.
Defining the fluid saturation field f [x, t] is the key to consider- ing the flow front and distinguish between the fully saturated and unsaturated areas.
Liquid permeability which is the key in the macroscopic flow advancement by Darcy law is modeled by different approaches and is considered to be dependent on the deformation and cover the isotropy and anisotropy of the preform.
Considering the microscopic structure of the preforms, we developed a dual-scale flow formulation, inter and intra-ply (bundle) flows.The inter-ply flow is the flow through the wide channels between the plies, whereas the flow between the fibers inside the plies is the intra-ply flow.
In addition, constitutive relations concerning constitutive effective stress response of the fiber bed and elastic packing response of the plies are also modeled and considered using available idealistic models.
Considering all the challenges present in the model, different applications are considered and different manufacturing methods are tested using the routine.There are different factors that play important role during the simulations.Material data is one.The routine is quite sensitive to the material parameters.The choice of the data should be wise and realistic.The data used in our simulations are from the literatures in the case of vacuum infusion 19 ; and those related to the press-forming process are characterized from experiments. 20The relation between time step and element size is the second.The choice of the time step should be such that saturating an element, considering incompressibility of phases, liquid permeability, and post-flow modeling like porosity can be synchronized.For example, solving an element in one time step can alter the accuracy and results in numerical errors.We should also prevent overshooting or overfilling, where an element is going to be solved (or filled) in less than a time step.
For future development, there is potential for introducing a user interface to use this routine as both academic and

Disclosure statement
No potential conflict of interest was reported by the authors.
Finally, a staggered finite-element-based solution procedure was put forward for the total solution advancement, involving, on the one hand, the saturation degree-dependent porous media formulation, and, on the other hand, the assessment of the saturation degree using a post-processing of the Darcian velocity field.Also by introducing the dual-scale formulation of the flow-deformation and an anisotropic permeability model, using packing response of the plies, we are able to track the saturation in both micro-and macro-scale.

Nomenclature
right Cauchy-Green deformation tensor D i dissipation induced by 'drag'-interaction between phases D s dissipation produced by (homogenized) solidphase material D nvf dissipation developed by fluid compressibility deformation gradient {\mathbf{g}} gravity

Figure 1 Figure 2
Figure 1 Deformation of a two-phase continuum body.B 0 represents the undeformed region of the solid phase and B f 0 represents the fluid reference configuration and at the meeting point x, at each time instant t in B, we have the mixture continuum of solid and fluid phases

Figure 3
Figure 3 Hyperelastic packing of fiber content ϕ in dry region induced by fluid pressure of representative fiber ply in the fiber bed at a fixed value of micro-infiltration Ξ s

f [ , t] = 1 .
We thereby replace the strictly discontinuous free surface problem by a smooth transition of the liquid front in terms of the evolution of the fluid saturation field ξf = ξf [ , t]

Figure 4 Flow channel with orientation t and fiber stacksFigure 5 Figure 6
Figure 4 Flow channel with orientation t and fiber stacks

1
and flow Q on the outer surface Γ 0 .The next step is to establish the set of discretized finite-element equations pertinent to the present choice of interpolation of displacements and fluid pressure.In view of the weak form of momentum and mass balances (40)-(41) we obtain where explicit expressions for the unbalances e − ext e and e ] I = Ne u φe ⇒  =  = ]p I = Ne p pe ⇒  = NODE p ∑ I =1 N I [] I = Ne p ηe ,

Figure 7 LRIFigure 8 Figure 9 19 Figure 10 Figure 11 Figure 12
Figure 7 LRI of fiber bed principle for specimen subjected to vacuum induced pressure load

Figure 13
Figure 13 Relaxation test at globally undrained condition with respect to different initial loading rates