A 3D process simulation model for wet compression moulding Part

Wet Compression Moulding (WCM) provides large-scale production potential for continuously fibre-reinforced structural components due to simultaneous infiltration and draping during moulding. Due to thickness-dominated infiltration of the laminate, comparatively low cavity pressures are sufficient - a considerable economical advantage. Similar to other Liquid Compression Moulding (LCM) processes, forming and infiltration strongly interact during process. However, the degree of forming is much higher in WCM, which disqualifies a sequential modelling approach. This is demonstrated in this work via experimental characterisation of the interaction between compaction and permeability of a woven fabric and by trials with a transparent double dome geometry, which facilitates an in situ visualization of fluid progression during moulding. In this light, and in contrast to existing form filling approaches, a forming-inspired, three-dimensional process simulation approach is presented containing two fully-coupled macroscopic forming and fluid-submodels. The combined model is successfully benchmarked using experimental double dome trials with transparent tooling.


Introduction
Wet Compression Moulding (WCM) provides high-volume production potential for continuous fibre-reinforced components. Due to simultaneously draping, infiltration and curing, low cycle times can be achieved [1][2][3][4][5]. Due to thickness-dominated flow progression and short infiltration distances within the mould, comparatively low cavity pressures are sufficient, which makes the process attractive compared to conventional liquid composite moulding processes such as (Compression) Resin Transfer Moulding (RTM/CRTM) or Vacuum-assisted Resin Infusion (VARI). Similar to other LCM processes, experimental and numerical investigations regarding WCM processes proved strong mutual dependencies between involved physical mechanisms, especially between resin flow and textile forming [6,3,2,[7][8][9]. As illustrated by the principle process steps in Fig. 1, the individual fabric layers are cut, stacked (1) and resin is applied on the topside (2), commonly using a wide slit nozzle. The resin seeps in while the partly infiltrated stack is transferred to the mould. In the third and most important step (3) infiltration and textile forming take place simultaneously. Thus, injection is not a downstream process step in which the fabric is almost constrained. Moulding of partly infiltrated stack leads a strong, intrinsic interaction between viscous draping and infiltration mechanism as schematically shown in Fig. 1 (green box). However, the individual contribution and relevance of the involved mechanism varies during moulding. In this regard, the tool stroke can be subdivided into three phases (cf. Fig. 1). In the first approximately 85% of moulding (I), the partly pre-infiltrated textile stack is shaped. This phase is dominated by draping mechanisms which are related to the membrane, bending and interface behaviour of the partly infiltrated material. When compaction begins, relevant fluid pressure starts to emerge in some areas of the mould. This coupled phase (II) is critical in terms of Fluid-Structure Interaction (FSI) since the stack is not fully constrained yet, but already locally exposed to increased fluid drag forces. Finally, a fluid phase is entered (< 5 %) which is dominated by infiltration effects (III) such as porous through-flow and curing. Here, draping, except compaction, is almost completed. Processing times typically range from 1.5 to 3 min [6]. In the following, related experimental and numerical studies are outlined with special focus on woven fabrics as this material is used in this study.

Related experimental work
Understanding of the occurring physical mechanisms is crucial for proper process development and final part design -in practice as well as by simulation. Starting with draping, numerous studies concerning the different deformation mechanisms are conducted. Comprehensive reviews in this field are for example published by Boisse et al. [10,11], Gereke et al. [12] and Sachs et al. [13]. Preliminary work with special regard to the implications of pre-infiltration on membrane [7], bending [8] and contact [14] behaviour proved significant influences during moulding.
During the last decades, compaction (consolidation) of composites has been subject of many studies, as it is of general concern in many composite processes. This applies to the field of thermoplastic unidirectional (UD) tapes [15] as well as to liquid moulding processes [16][17][18][19][20][21][22]. It is generally acknowledged, that a characteristic non-linear shape for the reaction forces can be expected when engineering textiles are subjected to transverse compression. Beyond that, even higher reaction forces are obtained when dry textiles are additionally subjected to shear deformation, which originates from increased material density as outlined by Ivanov and Lomov [23]. Moreover, obtained reaction forces increase further when pre-infiltrated specimens are exposed to compaction because of the superimposed fluid pressures below the punch as for example shown by Bickerton and Abdullah [20]. Additional compaction related effects, which are neglected in this study, are stress relaxation, rate-dependency and hysteresis effects [20].
Strongly interconnected with compaction, determination of permeability is of high importance in all LCM processes. Among others, Arbter et al. [24] and Vernet et al. [25] provide comprehensive reviews and recommendations. On macroscopic scale, Darcy's law [26] is usually applied to describe the fluid progression within a homogenous porous medium. Several studies are conducted with scope on the anisotropy arising from fibre reorientation in woven fabrics during shear deformation [27][28][29]23,25,[30][31][32]. Demaria et al. [29,28] showed, that geometrical reorientation of fibres and change of fibre volume content are the two factors responsible for the change of permeability during deformation. Beyond that, a complex interaction of superficial squeeze flow [33] and porous media through-flow arises inside the mould within the last phases (II,III) of the tool stroke (cf. Fig. 1). The same mechanisms are observed during CRTM, VARI and other LCM processes [34,35]. An experimental assessment of this interaction is considerably more demanding and industrial application therefore tries to avoid free surface flows wherever possible [3] for the sake of process robustness. Still, it cannot be prevented completely, especially when dealing with more complex components or unsuitable process parameters. Moreover, it can even be beneficial in terms of fast infiltration of large structures which is for example utilized in CRTM.

Related numerical work
Draping. Usage of Finite Element (FE) based forming simulation enables a detailed analysis of the resulting deformation by means of constitutive modelling of the material behaviour, considering process and boundary conditions via a Lagrangian description. Modelling of macroscopic forming behaviour is based on constitutive modelling of the relevant deformation mechanisms, which are usually separated into intra-ply mechanisms for the single plies (membrane and bending) and interface mechanisms between the stacked single plies [36,37]. Several modelling approaches have been published in recent years. Regarding engineering textiles, studies have been presented among others by Boisse et al. [38], Hamila et al. [39], Dangora et al. [40], Schirmaier et al. [41], and Harrison et al. [42]. Comprehensive reviews regarding FE forming simulation modelling approaches are presented by Boisse et al. [11,43] and Bussetta and Correia [44]. Several of these models are already applied in industry as part of commercial available codes, for example PAMFORM (ESI-GROUP), ANIFORM (ANIFORM ENGINEERING B.V.) or recently SIMUDRAPE (SIMUTENCE). Given the fact that all forming simulations contain large contact areas between the individual plies, explicit time integration schemes are almost always applied for this type of simulations. A variable thickness to account for transverse compaction is still state of research [45][46][47]. Unfortunately, this compaction is a key mechanism during WCM. Thus, an enhanced continuum-shell-based approach is proposed in this work, which is well suited for minor challenging mould geometries.
Mould filling. Similar to the field of forming simulation, numerous models for form filling applications have been proposed in recent years [48][49][50]. Their capabilities strongly depend on their field of application and the amount of allowed deformation. Accordingly, several modellings strategies and solving techniques are proposed, ranging from usually applied Finite Volume (FV) [51,52] techniques or Finite Differences (FD) [53] codes to Finite Elements with additional control volumes (FE/CV) [54][55][56][57][58][59] and coupled approaches using discrete interfaces such as the Eulerian-Lagrangian-Methods (CEL/ALE) [60]. The disadvantage of the latter one is the missing fluid exchange at the stack interface between free surface fluid and porous medium. Usually, Stokes flow is applied to account for the superficial resin, whereas Darcy's law Fig. 1. Graphical abstract | WCM process steps, crutial mechanisms and aim of this study. C.T. Poppe et al. is utilized to model porous media flow [35]. Implicit integration schemes are usually applied in order to achieve a stabilized solution for fluid velocity and pressure field [50]. Similar to the draping models, established modelling approaches are exploited in industrial applications for example as part of PAM-COMPOSITES (ESI-GROUP) for continuous reinforcements which covers most of the LCM processes -except the WCM process.

Motivation and structure of this paper
In this work, a macroscopic, fully coupled, three-dimensional process simulation model for WCM is proposed, driven by a forming-based approach. In contrast to existing forming simulation approaches, a suitable forming submodel, able to account for non-linear compaction behaviour is presented, parametrized and verified by experimental trials using a transparent tooling system. Subsequently, a fluid submodel is presented, able to describe the simultaneous fluid progression during moulding in a homogenized manner, even when the domain is subjected to finite strains. The combined process model is benchmarked with experimental trials on part level. In contrast to the ahead outlined fluid models, an explicit time integration scheme is applied to solve for the pressure field within the fluid. This allows for the fluid submodel to be directly solved along with the forming submodel which significantly enhances numerical performance as all iterative cycling is prevented. In the following, the process related terms "infiltrated" or 'pre-infiltrated' are used consistently throughout this work, to address wetted, impregnated or soaked woven fabric from micro-to macro-scale.

Experimental trials
In the first part of this work, experimental trials are conducted to provide an insight into the key physical mechanisms during WCM and provide characterisation and verification data for the process simulation model presented subsequently. A non-bindered 12 K carbon fibre woven fabric (T700SC-12 K-50C by ZOLTEK) in combination with a silicon oil is utilized. The latter provides a rate-and time-independent viscosity with a standard deviation of 5mPas. Due to previously conducted studies with this material combination, experimental data regarding membrane [7], bending [8] and contact behaviour [14] are already available for dry and pre-infiltrated states. Moreover, first results for the in-plane permeability of the unsheared fabric are reported in [5]. Still, key material information about transversal compaction behaviour and permeability during shear-deformation is missing.

Compaction
Compaction behaviour is analysed using a simple punch-to-plate setup mounted on a universal testing machine as presented in [5]. Unsheared fabrics are investigated along with pre-sheared fabric with 25 • and 50 • . The latter is close to the locking angle of this textile. A layup of [0] 4 is used, yet to demonstrate nesting effects, a second stack with different relative fibre orientation [0/45] 2 is investigated additionally.
A constant punch speed of 5mm/s is applied during all trials. Ratedependency is excluded in this study. In addition to dry stacks, preimpregnated stacks are investigated using a predefined amount of silicon oil with a constant viscosity of 20mPas. To ensure the fluid being pushed through (porous medium flow) and not above (superficial flow) the porous fiber stack, the punch is slowly brought into contact with the top ply before starting a trial. Each measurement is conducted 3 times. The recorded reaction forces during compaction are depicted in Fig. 2. Starting from the left to right, the impact shearing (a), pre-infiltration (b) and nesting (c) on the reaction forces are displayed. An expected, non-linear behaviour is observed during all trials. Compaction stiffness increases along with increasing shear angle. As stack height and material density are increased by the shear deformation, stack stiffness increases much earlier. Similar effects are observed under the impact of preinfiltration. The superimposed fluid pressure below the punch leads to a shift towards higher compaction forces. In contrast to the dry specimens, the pre-infiltrated result curves of different shear angles do not have a significant offset from another. This indicates, that during viscous compaction the presence of incompressible fluid is more important than the increased transversal stiffness of the sheared specimens. Whereas the dry results are used to parametrize downstream presented the compaction model, viscous compaction results will be used to verify fluid-structure-interaction in forthcoming work. Finally, the impact of nesting (c) is outlined in Fig. 2. Nesting reduces the stack stiffness by fibre overlapping at the interface. As outlined in the plots, this effect is more pronounced when fibres are aligned in parallel ([0] 4s ).

Permeability
Permeability, which is later needed to account for fluid progression throughout the porous stack is measured using a radial test setup with circumferential clamping. The setup has already been described and verified in [5]. Similar to comparable tests [25], a one-dimensional formulation of Darcy's law is applied to calculated the permeability based on the pressure difference between inlet and outlet in stationary regime. A constant inlet flow rate is applied in this work. This implies, that the completely saturated and not the transient permeability is characterized [61]. In [5] we verified this procedure via a direct comparison between a test setup measuring the transient permeability and the here applied one which show no significant differences for the permeability values for the in this study applied material. This cannot be assumed for every material [61]. In addition to undeformed results, the setup is used to measure permeability of four ply stacks, priorly subjected to conditioned shear deformation [25]. This provides an alternative method compared to the commonly applied picture frame-like  clamping frame as for example presented in [27,28,23,32]. The measured areal weights (w γ 12 A ) corresponding to in-plane shear angle γ 12 are summarized in Table 1 and are found to be in good agreement with w A,γ = w A,0 /cos(γ 12 ) [28,62].
In this work, a forming-driven approach is pursued, where usually non-orthogonal material frames are utilized [11]. We take advantage of this in the following by expressing the permeability tensor in terms of such a frame. Fig. 3  Usage of an orthogonal the permeability tensor K, which is reasonable in fluid mechanics, leads to several challenges when the fabric is subjected to shear deformation [28,29]. First, an anisotropic permeability in both ellipse directions (K 11 , K 22 ) with non-constant ratio arises. Secondly, the orthogonal ellipse frame rotates relative to the first fibre orientation, often referred to as ellipse rotation angle (cf. β in Fig. 3). Preliminary trials revealed a circular flow front progression for the undeformed balanced fabric K 0 X1 = K 0 X2 , or K 0 f1 = K 0 f2 , respectively. As outlined in Fig. 3 it can geometrically be shown, that the aspect ratio a = const. = K f1 /K f2 between both fibre direction remains constant when the non-orthogonal material frame is applied. This is supported by the findings of Lai and Young [27] for a comparable material. In this regard, a direct calculation of the material frame permeability value K f1 based on the combined considerations of both fibre directions is enabled via Result. The experimental result are plotted in Fig. 4. Every dot illustrates an investigated configuration. On the left side (a), the experimentally determined permeability K X1 with respect to the fix orthogonal frame { X i } for different shear angles is plotted in grey lines. The permeability values in the orthogonal frame K X1 strong decreases with increasing fibre volume content (fvc) and slightly increases with shearing. The same permeability values are also plotted with respect to the current fibre orientations K f1 in the non-orthogonal frame { f i } using green lines. An overlap of the resulting green curves with the undeformed values (black curve) indicates, that usage of a nonorthogonal frame can bypass the above outlined challenges resulting from shear deformation in combination with orthogonal frames, because it intrinsically captures the underlying mechanism outlined by Demaria et al. [29]. Meaning, the material frame accounts for the fiber reorientation of the permeability values. The change of fvc is introduced via the above outlined measured area weights at different shear angles. This is later covered by the volume integration within the numerical model. This reduces experimental effort significantly. A verification of the resulting anisotropic ellipse-shaped flow front progression during shearing is given in Section 3.3.3. Transverse permeability. Furthermore, permeability in thickness direction is needed to account for gravity-and compaction-induced flow progression (cf. Fig. 4(b)). For this, measurements are conducted at the Fraunhofer IGCV in Augsburg (Germany). A test bench presented by Graupner and Drechsler [63] is applied using a 12 layer stack. The trials are conducted using a multiply stack compacted in thickness direction between two spacer plates. Conventional distribution structures (perforated plates) distribute the inlet flow. Similar to in-plane trials, a one-dimensional Darcy flow in stationary regime is assumed for evaluation. The permeability is evaluated at different compaction states via the resulting pressure difference within inlet and outlet camber. The large number of plies is necessary due to the gaps of the 12 K material. The values for transverse permeability are significantly lower than the values obtained for the in-plane directions. However, it is assumed that the values would be much higher in case of the four ply stack, which is investigated in this work, because of the rovings gaps which are more likely to create flow channels if only a few are stacked. The findings are approximated via model functions using a regression approach for the fibre-frame permeabilities K f1,f2 (ϕ) = 0.029e − 13.755ϕ and K f3 (ϕ) = 0.001293e − 15.950ϕ for subsequent simulation.

Transparent double dome trials
Experimental trials on a double dome shaped transparent geometry are conducted to achieve a better understanding of interaction between stack deformation and infiltration in WCM processes. To enable an optical insight during moulding and simultaneous fluid progression, a transparent double dome geometry as shown in Fig. 5 (a) is used. The tool was manufactured by PELZ TECHNIK GMBH in Heilbronn (Germany). An electronically controllable cylinder with integrated position control (EMMS-AS-140-L-HS-RSB motor in combination with ESBF-BS-100-400-5P cylinder by FESTO) is used to close the mould precisely. A circumferential blank holder is included into the upper mould. Fluid progression is observed by two cameras, one focused on the top, the other one filming from the bottom side by means of an 45 deg mirror. All trials are conducted using a two-parted four second closing amplitude to reduce the tool cavity from 62 mm to 2 mm. During the first two seconds a constant closing speed of 25mm/s is applied which is reduced to 5mm/ s for another two seconds when a remaining cavity gap of 12 mm is reached. The amplitude is inspired by commonly applied WCM process amplitudes, where closing speed is strongly reduced towards the end of moulding. Additional sealing are placed inside the upper and lower tool to prevent leakage.
To conduct a trial, an undeformed four ply stack [+45/ − 45] 4 is placed on the lower tool. A distinct amount of coloured silicon oil (120ml) with a constant viscosity of either 20 or 200 mPas is centrally applied on the top layer. Viscosity values are chosen to account for reasonable process limits. To ensure reproducibility, a constant soaking time of 15 s or 30 s in case of the 200 mPas trials is waited for before  Results. Obtained results using the low viscous oil are presented in the top row for discrete, remaining cavities indicated by the numbers (1)- (6). All pictures are divided into top and bottom view to allow direct assessment of the effects in thicknesses direction. During the trials with the low-viscous fluid (20 mPas), though-flow in thickness direction leads to an accumulation of fluid in the lower tool as the top left pictures (1,2) show. Although this is only a bit compared to the total amount, it unfortunately prevents an optical evaluation from the bottom side for these two picture. The gaps within the four ply stack facilitate gravity-induced infiltration of the fluid thickness direction. Following the tool stroke from left to right, this implies that the stack is homogeneously infiltrated from both side. This can for example be illustrated by the top row pictures (3) and (4), corresponding with cavity heights of 4 to 3 mm. The fluid which has accumulated in the lower mould is forced back into the stack and only at last, fluid is pushed over the flanges.
In contrast to the trials with the low viscosity fluid, the high viscosity of 200 mPas prevents leakage in thickness direction prior to moulding. No fluid, apart from first drops can be detected at the bottom view until a cavity height of 4 mm is reached, as shown in the bottom line picture (4). Consequently, the stack is infiltrated mostly from the top side as the pictures (3) to (6) indicate, especially when compared with the lowviscosity results. The fluid progression appears to be more superficial (200 mPas), because the fluid is pushed over the flanges much earlier at the top side. In the same manner, the drops at the bottom side lead to some infiltration regions in the dome area, which is later overflown by a more superficial based fluid progression from the centre (cf. (3)(4)(5) in the bottom row). Although the final infiltration state does not differ significantly between both viscosities, the infiltration path is significantly different and the infiltration remains more top-ply-based for the 200 mPas trials. This is reasonable, because the thought-flow resistance in thickness direction increases with viscosity, which favours the surface flow at the top ply.
In addition to the fluid progression within the mould, resulting reaction forces are of high interest for the process development in various liquid moulding processes [48,20,21,64]. Recorded reaction forces for the three configurations, dry, 20 and 200 mPas are displayed in Fig. 5 (b). In agreement with existing literature and the outlined experimental compaction results, higher forces are obtained with simultaneous fluid progression and increasing viscosity.

Process simulation model
The experimental results demonstrate the complex interaction between stack deformation and simultaneous fluid flow during moulding. As outlined in the introduction Section 1, current modelling strategies for moulding processes often choose a Eulerian or partly Eulerian approach [50] extended by a suitable forming or compaction behaviour. In this work, an opposite approach is being pursued, because for the WCM process, forming-induced mechanisms such as membrane,    bending and interface friction are predominant during most part of the tool stroke. Even though fluid progression is of equally high importance towards the end of the tool stroke, transient material conditions such as fvc and compaction directly rely on the local deformation state. The modelling strategy proposed in this work is guided by a Lagrangian view based on a reliable and suitable forming model, which is than enhanced to account for simultaneous fluid progression during moulding. In this regard, a solution based on the explicit solver of ABAQUS is presented. On one hand, this avoids interfaces between different solving algorithms, on other hand, this allows for a numerical efficient, fully parallelized and coupled modelling approach.

Modelling strategy
In this light, the modelling approach displayed in Fig. 7 is pursued. Starting from the left side, a schematic section cut of a viscous forming trial is shown. As common in the field of FE forming simulations, tools are described as rigid surfaces and the stack is represented by individual single plies, which are interconnected using a suitable contact algorithm. Additionally, fluid progression across ply interfaces is required, which is referred to as Solid-Fluid-Contact. Thus, each single ply is represented by two coupled submodels, a forming submodel which accounts for displayed deformation mechanisms, and a fluid submodel which addresses mechanisms such as forced infiltration and fluid contact. Superficial fluid is not directly accounted for by explicit modelling in this work, but implicitly via an increased stack thickness.
Several user subroutines are mandatory and combined during runtime in order to implement this approach within the explicit solver of ABAQUS (cf. Fig. 7 (2)). An VUEXTERNALDB subroutine is called at several locations and acts as a master control layer to gather and share model information, ensure thread security and compute results. The two submodels are implemented within a user material routine VUMAT superimposed with a user element routine VUEL, which enhances the compaction behaviour of the forming submodel and contains the fluid submodel. Beyond that, a nodal field-based subroutine VUFIELD and an integration-point-based subroutine VUSDFLD are used to access and share informations during runtime (3).
The forming submodel itself addresses membrane, bending and compaction behaviour, via a built-in Continuum shell element (SC) with user membrane material and superimposed user-enhanced compaction behaviour. Usage of SC elements provides the advantage of a lockingfree three-dimensional element, which directly makes use of a inplane material subroutine, however at the cost of coupled membrane and bending behaviour and a constant linear thickness modulus. In general, assessment of compaction and its accountancy within a fully decoupled three-dimensional formulation during forming is especially challenging and it is often neglected in current forming simulation approaches because it requires a significant more complex element and material formulation in order to avoid numerical locking effects [45][46][47]. In this work, decoupling between membrane and bending behaviour within each ply is neglected as SC elements are utilized. This does not imply, that bending behaviour as a whole is neglected. Rather, neglected decoupling means that the bending stiffness cannot be formulated independent from the elasticity modulus of the membrane part. This is commonly done in forming simulations to contribute to the transversal shear between the filaments and rovings [37,11]. In Section 3.2.3 we demonstrate, that bending stiffness of each ply can still be approximated quite accurately when the transversal shear stiffness of the SC element is adapted accordingly. This procedure is valid for minor challenging geometries and forming setups, including the here investigated double dome geometry. Moreover, the applied circumferential blank holder favours membrane deformation even further.
Despite the bending behaviour, a constant linear thickness modulus proves a rather poor approximation of real material compaction behaviour as illustrated in Fig. 2. Therefore, the compaction behaviour needs to be enhanced by a mechanical part within the superimposed user element that supplies a constitutive equation in thickness direction.
The fluid submodel is embedded in a VUEL subroutine which is superimposed to the SC element. Addressed key mechanisms of the fluid model are anisotropic permeability, compression-induced infiltration, transient interface flows and implications of gravity. Heat conduction and curing are neglected in this first approach.
The user element is formulated in an isoparametric manner by means of linear shape functions matrix N. The operator matrix in the current configuration x is calculated using the spatial derivatives of these shape functions with respect to the current coordinates. The deformation gradient F within the elements is computed using the Jacobian matrix J 1 with respect to the current configuration x and the Jacobian matrix J 0 referring to the reference configuration X.
For material modelling, the Green- as a measure of deformation in the initial configuration is applied, where I provides the second order unit tensor. The Green-Lagrange strain rate Ė is calculated by means of the current time increment size ▵t inc and previous values. Moreover, the rate-of-deformation tensor is used to calculate strain rates with respect to the current configuration.

Constitutive modelling
In the following, constitutive equations are presented to account for the outlined deformation mechanisms. Test cases and comparison with experimental data is used to parametrize and verify the outlined model.
Membrane behaviour. The implemented hyperviscoelastic formulation is based on invariants of the Green-Lagrange strain tensor E and the rate-of-deformation tensor D presented more detailed in [7]. Two deformation modes, elongation in both fibre directions and shear, are taken into account and thus contribute to the potential of free energy Ψ according to where I elast is calculated by superimposing the elastic and viscous shear part following an Voigt-Kelvin approach. The second Piola-Kirchhoff stress tensor S is transformed (pushed) to the material fixed frame x via the stretch tensor U to provide the corresponding Cauchy stress with the Jacobian determinant J = det User-enhanced compaction behaviour A hyperviscoelastic orthotropic three-dimensional formulation using the Voigt-Kelvin model is implemented within the user element to account for a non-linear compaction behaviour under finite strains. While the elastic part only affects the thickness direction, an isotropic rate-dependent numerical damping is superimposed in all directions for stability reasons. Voigt notation is used in the following and indicated by σ = [σ 11 , σ 22 , σ 33 , σ 23 , σ 13 , σ 12 ] ⊤ .
For the elastic part, Kirchhoff stress τ el is calculated by means of deformation gradient F , Green-Lagrange strain Ẽ and the stiffness matrix C according to: Only the third component of the stiffness matrix , which only contributes within elements that are stronger compressed than the experimentally acquired data range. By this means, collapsing is prevented for temporarily high compressed elements via: The material part E Mat 3 is superimposed to a required base stiffness E SC 3 = 7.156⋅10 − 4 N/mm 2 of the SC element for compaction strains above ε 33,0 .
Parameters are obtained via comparison with the experimental data for shear angles of 0 ∘ , 25 ∘ and 50 ∘ and are summarized in Table 2. Shear dependency is accounted for by linear interpolations between the data sets. The viscous part of the Kirchhoff stress τ v = 2η Mat 0 Dis calculated via the rate-of-deformation tensor D and is used to account for a basic material damping in all direction. A rate-dependency of the pre-infiltrated compaction trials as known from literature [20] could simply be included here, but this is neglected in this study due to a lack of experimental data.
Following the Voigt-Kelvin approach, both elastic and viscous part of the Kirchhoff stress are superimposed according to τ =τ el +τ v , which are subsequently used to determine the nodal reaction force matrix f R with respect to Green-Naghdi's frame used by ABAQUS via Bending behaviour. The bending behaviour is supplied by the built-in continuum shell (SC) element in ABAQUS which includes the membrane behaviour via a user material routine. However, it is only an approximation of the real material bending behaviour, because membrane and bending behaviour cannot be described in a decoupled manner. This leads to an overestimation of the bending stiffness when the conventional shell theory is applied, because of the high prescribed in-plane fibre stiffnesses. To partly compensate this, the transversal shear stiffness is reduced (see Section 3.2.3). Furthermore, usage of the SC element implies a purely linear-elastic, isotropic bending behaviour which neglects the anisotropy caused by the fibre orientations and ratedependencies. This is suitable for geometries where bending is of minor importance compared to membrane or contact behaviour, and where the material itself provides a comparable bending stiffness in both directions. This applies for the double dome geometry with a circumferential blank holder where membrane behaviour dominates.

Solid-interface behaviour
Contact behaviour between the plies is implemented by a built-in general contact algorithm of ABAQUS, using friction coefficients of μ ToolPly = 0.29 and μ PlyPly = 0.223 [14]. The penalty contact normal stiffness E IFace 3 = 20N/mm 2 is set with the intention to reduce penetration between plies and tools as much as numerically possible during the whole moulding process to ensure suitable compaction values within the plies. Usage of the built-in general contact algorithm implies an isotropic in-plane friction behaviour which is a common assumption for woven fabrics.

Parametrisation and verification of the forming submodel
User-enhanced compaction behaviour. The two plots in Fig. 8 provide the combined results of the shear-dependent parametrization of the compaction behaviour. For the sake of clearness, only the mean curves of the experiments are presented. Standard deviations are already displayed in Fig. 2. A good agreement is observed. The benefits of the userenhanced compaction behaviour becomes particularly obvious when comparing the outlined results with an fitted built-in solution with constant thickness modulus (see Fig. 8(b)). In any case, a constant modulus would lead to an overestimation of the initial stiffness as well as an underestimated locking towards higher compaction states. The identified data sets are summarized in Table 2.
Bending. A cantilever setup (cf. Fig. 9 (a) is used to obtain suitable bending parameters for the transversal shear stiffness. Target values for the cantilever bending length l bend. in fibre direction are obtained by ) experiments outlined in [8] for the same material. Using this approach, transverse shear stiffnesses are determined to κ 11 = κ 22 = 0.03 N/mm 2 and κ 12 = 0.03 N/mm 2 . To evaluate the suitability of this procedure, a more complex geometry than the double dome is used to compare a properly decoupled reference model [8] based on superimposed membrane (M3D3) and conventional shell elements (S3R) with the above parametrized SC model (SC6R) with approximated bending properties.
The bending within the reference model is reduced to a purely elastic material behaviour to ensure comparability. In the same manner, an upscaled, compaction modules E SC 3 = 10 N/mm 2 is used to reduce the impact of possible compaction. Fig. 9 (b) provides a direct comparison of the two approaches for a two-layer forming setup exposed to gravity. Apart from minor differences in slag and shape, a comparable displacement field is obtained. Hence, the SC-based approach provides a    suitable approximated bending behaviour in this complex case, even though a decoupling is neglected. Since the experimentally applied double dome geometry is simpler and contains a circumferential blank holder that constrains the draping even further, the outlined approach is well suited to capture the forming behaviour of the applied double dome geometry.

Constitutive Modelling
Given the objective that the fluid submodel should be superimposed to a solid-mechanical FE-Forming mesh in forthcoming analyses, a wedge-shape element as illustrated in Fig. 10  Gaussian integration scheme in conjunction with six integration points is applied to ensure an accurate integration of the system diffusion matrix K diff during large strains. Following the modelling strategy outlined in Fig. 7, a fully coupled thermo-mechanical analysis within ABAQUS is utilized to model forming and fluid progression simultaneously [59,5]. Doing so enables the additional usage of nodal pressure p N i degree of freedom. Overall implementation is mainly concerned with the following three key mechanisms [5]: 1. Solving of the pressure field within the saturated domain 2. Flow front progression and mass conservation 3. Modification of the pressure field during compaction Saturated domain. To solve the pressure field, Richard's law [65], which provides a combination of mass conservation and Darcy's law is applied in a transient explicit formulation By this means, transient single phase fluid progression can be modelled with an implicit accordance of the velocity field by introducing a nodal hydraulic capacity C N hyd which is used to calculate the mass matrix M mass = ρC N hyd I→0 of the system. In context of our study, C hyd provides a purely numerical parameter, which is required to apply an explicit temporal integration scheme to solve the pressure field. Physically, it relates the change of volume (dV) to the change of pressure (dp). Incompressibility of the fluid would lead to infinite pressures in case of the smallest numerical error during temporal integration. The latter cannot be prevented with finite precision. Hence, we introduced a slight compliance (1≫C hyd > 0) into the system, just enough to prevent this.
Within one element, viscosity η, ϕ and fibre-frame permeability K f are constant. In Eq. (6), an 'e' is superscripted as an identifier for a nodal vector representation according to p e = [p (1) , p (2) , …, p (N) ] ⊤ with N = 1, 6. For a suitable representation of a quasi-incompressible fluid, M mass should be as small as numerically feasible [59]. The diffusion matrix K diff is calculated based on fibre frame permeability tensor K f in the current deformation state x i . Gravity is accounted for by the body force vector b grav,e . The operator matrix B provides the spatial derivatives of the shape functions in the current configuration. The Jacobian determinate det ⎛ ⎝ J ⎞ ⎠ accounts for a change of element volume. Finally, a nodal source vector S +,e is introduced to modify the pressure field.
Flow front progression and mass conservation. A physics-based implementation of the flow front progression is needed to ensure mass conservation. Incremental flows V (k) across every element surface k (cf. Fig. 10(a)) are calculated by means of a midpoint Gaussian integration of the surface integrals where n (k) provides the current normal surfaces vector. Thus, the accumulated incremental volume flow across surfaces of the same element is zero and their location (inlet, vent, wall or flow front). The latter is additionally used to account for process boundary conditions. Further, an algorithm based on the pressure gradient and the free area within the filling elements is implemented to determine the partition of the summed propagation flow. By this means, a mesh-independent flow front progression is achieved [59]. Modification of the pressure field. An essential mechanism is the local response of saturated elements to change of volume, for example during compaction. The local overflow V + The underlying idea of this concept is that at any point in time (increment) for which a stable pressure distribution is reached, the summarized incremental new overflow must be equal to the incremental summarized amount of fluid that trans-passes the flow front. In this manner, overflow is only reduced when the fluid has actually crossed the flow front. This further implies that the overall pressure field is forced to adapt to the current overflow and the conduction situation as a whole. This makes the solution for the pressure distribution unique.
In contrast to implicit time integration methods which are commonly applied for fluid models to enable a converged solution of flow velocity field and pressure distribution, the outlined explicit approach uses a forward modification loop of the pressure field, which only considers the fluid velocity field implicitly. This approach provides an elegant and efficient way to consider the fluid velocity field without actually solving for it directly. Yet, the presented approach only provides reliable results if the system is able to adapt its pressure distribution sufficiently fast. The verification cases presented in Section 3.3.3 along with the results presented in [5] demonstrate that this applies, when temporal integrations steps ▵t inc are small, the hydraulic nodal capacity C N hyd , which represents a purely numerical parameter in this context, is sufficiently low and penalty parameters β 0 and β 1 are chosen adequately. During explicit FE-forming simulations, sufficiently small temporal integrations steps (increments) are ensured intrinsically. However, an additional stabilization mechanisms on domain level is implemented and demonstrated in Section 3.3.4. Hydraulic capacity and the penalty parameters are optimized during the verification test cases. For this purpose, the total incremental overflow is minimised during the virtual trials.

Interface flow modelling
To account for fluid progression in thickness direction caused by gravity and/or compaction, interface flows between the single plies must be captured. This implies the detection of the current flow front's position across interfaces as well as pressure conduction between the nodes of different model instances within the saturated domain. The latter can be achieved via an additional gap conduction option α Gap for the already presented built-in contact model in Section 3.2.2. To identify the flow front, each element must be able to identify not only its own status (saturated, filling or empty) but additionally the status of current neighbours in contact. This can be complicated by large slip in combination with large deformation of the individual instances which are eminent to forming simulation. Moreover, information consistency must be ensured between interfaces that may be part of other processor threads during runtime.
An algorithm is implemented which allows for a multi-thread optimized neighbourhood search by means of scanning spheres for the individual nodes and the element midpoints. During initialization, neighbour elements are identified within the plies but also throughout the contact interfaces. A relational network of element, surfaces and node connectivity is created via combined information obtained from subroutines. As scanning of multiple elements can require a considerable amount of computation time, dynamic pre-partitioned lists are utilized, which are constantly updated during runtime. However, some simpler geometries or forming setups, where only minor slip is induced between the plies do not need an advanced transient contact tracking. Thus, a 'static' contact mode is additionally implemented, which assumes a constant element connectivity at the interfaces.

Verification tests for the fluid submodel
Anisotropic fluid progression. The experimental results presented in Section 2.2 indicate, that the implemented fibre-parallel description of the permeability tensor K f could be beneficial because it intrinsically account for permeability changes during shear deformation. To prove this, an anisotropic flow progression within a pre-sheared balanced fabric is considered and compared with the experimental results outlined by Pierce et al. [32]. They characterised the permeability values of differently pre-sheared fabrics by means of a central injection between transparent plates. Since they later applied ANSYS FLUENT for the flow simulation, they extracted the permeability values using the ellipse method [29] as shown in Fig. 11(a). A comparable result can be reproduced with the presented fluid submodel as Fig. 11 (b) illustrates, even though only values of the undeformed permeability are provided. In both cases, an ellipse-shaped flow front progression and ellipseshaped pressure distribution establishes. For clarification, an ellipseshape is added to the simulation result in Fig. 11 (b), which provides a comparable rotation angle of β = 25 ∘ between warp fibre orientation and principal ellipse axis. The comparability of the results supports Demaria's work [29,28], because it numerically proves, that in case of a woven fabric, anisotropic permeability caused by in-plane shear deformation is intrinsically accounted for, when a fibre-parallel description of the permeability tensor is applied (cf. Fig. 3).
Two-dimensional compaction. Compaction-induced infiltration provides an important mechanism during WCM. A verification of the basic model formulation using a linear test case is already reported in [5]. In case of a radial setup with a predefined amount of central applied resin, the radial flow front progression r f (t) = directly follows from mass conservation. Among others, Bickerton et al. [20] provides the following equation for the transient pressure distribution For this verification test, a constant compression rate ḣ = 0.055 mm/s is used. Initial resin distribution is set to a radius of r = 60 mm and a constant fibre orientation of 0 • is used. An initial height of 0.4 mm is applied. While compaction proceeds, resin is forced to infiltrate the dry fabric areas. The isotropic in-plane permeability of the unsheared deformation results in a circular propagation of the flow front (cf. Fig. 12 (a)). A structured mesh is used, thus mesh independence is proved by the circular flow front shape. The flow front progression (a) is in good agreement with the analytical solution. Beyond that, pressure history (cf. Fig. 12 (b)) at two distinct points p 0 (2 mm) and p 1 (70.7 mm) is compared to the analytical solution. To demonstrate the numerical impact of an increased hydraulic capacity of C hyd = 1⋅10 − 2 mm 3 /MPa, corresponding curves are added to both plots. Whereas the flow front is not noticeably affected, pressures tend to lag slightly behind when hydraulic capacity increases.

Domain stabilisation
A strong local source term S +,e in combination with a low hydraulic capacity C N hyd lead to a oscillation-prone system. Moreover, the strong coupling between deformation and source term makes it challenging to effectually estimate a stable time increment in a common manner. As a monolithic system combining both submodels is solved, a stable time increment is required which conserves stability for both.
Therefore, an adaptive scaling approach based on the forming submodel-based time increment estimation is implemented, using an instability criteria. In order to judge stability, the evaluation of slightly negative pressure values at the current flow front has proven to be a very useful indicator. If instability is detected among saturated elements, the current estimated time increment of the forming submodel is downscaled. Under stabilized conditions, the time increment is upscaled again. To demonstrate the suitability of the presented stabilisation approach, the above outlined radial test case is exposed to different stability scenarios via modification of the nodal hydraulic capacity C hyd (cf. Fig. 13).
While the fluid propagation progresses, emerging instability is detected within all scenarios and the time increment is reduced accordingly. The smaller C hyd is set, the stronger the time increment is reduced. Although different hydraulic capacities slightly affect the predicted pressure history as outline in Fig. 12, all three scenarios provide a stable accurate solution. An error is raised if the domain cannot be stabilized.

Fluid interface
A one-dimensional test case in through-thickness direction is used to verify the fluid contact behaviour outlined in Section 3.3.2. For this, an eight stack ply is exposed to a constant inlet flow at the top ply according to Fig. 14(a). A gap conduction of 5000 W/mmK (thermo-mechanical functionality in ABAQUS) is used. In contrast to in-plane flow progression, the flow front must now progress through the interfaces and not only between element surfaces of the same mesh instance. Pressure history at the top nodes is monitored and compared with the analytical solution (cf. Fig. 14) (b,c)). Flow front progression and pressure history provide a good agreement to the analytical solution. Moreover, the pressure history (c) verifies the modification of the pressure field. As mass is constantly added to the top layer elements, a local overflow is generated which modifies the pressure distribution until a stable solution is established. A slight peak emerges during later transitions of the flow front. This is due to a slight oscillation within the one-dimensional pressure distribution. This is observed when one-dimensional flows and few elements are evaluated. This was not found to be an issue for increased element numbers or multi-dimensional flow cases. As the bottom ply is filled beyond 0.29 s and no fluid outlet is assigned, the pressure immediately increases considerably. In summary, the implemented fluid contact enables a suitable out-of-plane flow front progression between the individual plies of the stack.

Pure forming of dry fabrics
The presented three-dimensional forming submodel is applied to a double dome geometry as outlined in Section 2.3 to ensure its suitability. For this purpose, a simulation setup as displayed in Fig. 15(a) is utilized   based on the modelling approach displaced in Fig. 7. The tool setup comprises rigid surfaces to account for lower and upper tool, as well as blank holder and seals. The proposed forming model is implemented within every single ply of the four ply woven fabric stack with initial fibre orientations of ± 45 • . A predefined ply thickness of t = 0.7 mm is used, obtained by measurements on the stack. Boundary conditions are applied according to the experimental trials. A displacement is prescribed at the upper tool according to the closing amplitude presented in Section 2.3. Gravity is applied for the entire model, implying a closing of the blank holder under its self-weight.
A forming result and resulting observations are displayed in Fig. 15  (b). First, gravity leads to a slag of the fabric, followed by initial slight compaction below the advancing blank holder. Hereafter, forming starts by initial contact of the upper tool with the slagging stack. This initial contact mainly occurs at the half dome centre and the edges which leads to a localized compaction in these areas. Further forming results in a thickness decrease at the flanges. Finally, the upper tool gets into largeareal contact with the stack, leading to final part thicknesses between t = 0.3 and 0.55 mm.
To verify the forming behaviour, a comparison with results of Boisse et al. [67] for a comparable material and double dome forming setup is displayed in Fig. 16. An initial fibre-orientation of ± 45 ∘ is used in both cases. A high qualitative agreement can be achieved in terms of shear angle distribution and outer shape, even though material locking angle and compaction forces are probably not identical. Still, if fibre reorientation is correctly taken into account, a good agreement can be achieved. In addition to the spatial comparison with Boisse et al. [67], Fig. 16 (b) provides a direct comparison between the optically evaluated and the predicted shear angle at two critical positions (P1,P2). The former are obtained from the above outlined experimental trials. A slight overestimation is reported for P1 and a good agreement is reached for P2. Thus, an overall suitable prediction of the in-plane shear angle is achieved.
Furthermore, an evaluation of thickness and corresponding fvc ϕ within every ply is enabled by the proposed approach. Fig. 17 (a) provides an insight into the predicted final thickness and fvc in the top and bottom plies. Starting with the thickness, some areas provide quite comparable though-thickness results. However, several critical areas such as the flanges (I), the half-dome centre or the corners (II) do not provide similar results throughout the stack thickness. Regarding the flange area (II), increased thickness compared to the surrounding areas is predicted for the top ply (0.55) but the opposite is true within the bottom ply (0.35). This contradictory relationship applies for several areas and also the fvc distribution on the right side of Fig. 17 (a). In addition to the thickness information, ϕ also takes the local in-plane shear deformation into account (cf. Fig. 15 (b)). Thus, the highest fvc (0.61) is predicted for that part of the flange where the highest shear angles arise (close to I). In general, high fvc are predicted for highly compacted areas. Unfortunately, elastic springback of the uncured stack in thickness direction prevents a direct validation of the predicted thickness and fvc values for the completely closed tool. The predicted location of the maximum fvc agrees with the experimental observations, but we assume a slight overestimation of the maximal predicted fvc at Zone I caused by tad overestimation of the maximal shear angle.
Based on these results, compaction is of high importance to accurately predict the local fvc, especially in areas with only minor shearing. In our case, pure compaction results in a maximal fvc changes of 0.25 to 0.5 (100%), whereas pure shearing from 0 to 50 • with constant thickness increases area weight and thus fvc only by 31% (cf. Section 2.2). This is supported by the fact, that for most of the part, overall spatial distribution of fvc is much more similar to the thickness distribution than to the shear angle distribution. Thus, a reasonable prediction of the local fvc should include both shearing and compaction as an interconnected source for a non-linear increase of the final local fvc.
Since compaction is now accurately considered in the model, this enables to predict reasonable reaction forces during forming. In Fig. 17   (b) reaction forces obtained by dry forming trials are compared with the numerically predicted ones. A very good agreement can be achieved for a substantial share of the tool stroke. Significant differences only arise for the first 10 mm of tool motion, because gravity induced slag of the stack prevents contact to the upper tool until a cavity of 50 mm is reached. Moreover, blank holder pins are modelled friction-less in the simulation. During experiments however, some slip-stick effects between blank holder pins and upper tool are observed and measured within the resulting reaction forces for the initial 15 mm of the tool stroke. After that, the free length of the pins is short enough to prevent such effects of sticking. Minor differences can further be observed for the beginning of the final compaction phase at remaining cavity height between 5 and 3 mm. Here, the simulation lags slightly behind. This can have several reasons. First, the scattering within the experiments is quite high in this phase. Moreover, bending forces play an important role for the overall reaction forces in this phase. Finally, some degree of contact penetrations cannot be prevented between the interfaces with leads to a slightly more compliant system for the sake of numerical stability. Nonetheless, reasonable reaction forces are obtained by the threedimensional forming simulation submodel. Moreover, the need for a suitable three-dimensional locking free formulation for forming simulations is demonstrated.

Combined forming and filling
In the following, the outlined simulation setup of the double dome geometry is enhanced with the proposed fluid submodel (cf. Section 3.3.1), following the introduced modelling strategy outlined in Section 3.1. Boundary conditions are extended with the initial fluid allocation within the stacks which is based on the experimental observation. For the 20mPas trials, a homogenous through-thickness pre-infiltration within the centre stack area is used, whereas only the two top plies are set as infiltrated during the 200mPas trials. The latter addresses the mostly superficial distribution prior to forming (cf. Fig. 6), which leads an increased initial pre-infiltrated area to achieve the required amount of fluid. A result, combining both forming and fluid-progression is shown in Fig. 18. Similar to the experiments, the fluid flow is mostly driven by gravity until final compaction and resulting in-plane squeeze flow takes place towards the end of the tool stroke.
A direct comparison between experimental results and numerical prediction is presented in Fig. 19. For the low viscosity (20 mPas, top row), a good prediction accuracy is achieved in terms of flow front progression and local resin amount. Leaked fluid prevents a proper evaluation of the bottom side at remaining cavities above 5 mm. Moreover, this leads to minor differences between experiment and simulation. Regarding the 200 mPas trials (bottom row), prediction accuracy is slightly lower, because the flow progression is both more superficial and strongly inhomogeneous in thickness direction. Sill, comparable tendencies of the fluid progression can still be achieved, especially with regard to thickness infiltration. But, the initial resin distribution leads to an overestimated flow progression on the top side, which is a direct consequence of the increased initial infiltration area which is needed to reach the required amount of fluid.
Some experimentally related sources of uncertainty are to be considered. Only four plies are used during the experiments, which leads to some scatter in terms of fluid progression because the transversal  permeability can differ in some areas due to gaps. Moreover, video analysis of the transparent tools is challenging by curvature-induced optical effects that are reduced where possible with calibrations. In addition, some scattering must be anticipated in the initial fluid distribution during the experiments.
The application on part level provides clear conclusions about the achievements (+) and current limitations (− ) of the presented formingdriven approach: + Combined modelling of forming and infiltration as captured by the presented approach is essential for an accurate WCM process model. + Initial results are promising and confirm the general suitability of the presented approach as well as the need for a three-dimensional forming model which captures for compaction. + The monolithic approach allows for an explicit and efficient solving of the model without iteratively coupled algorithms. + Gravity is considered during both forming and flow-progression, which is required when modelling WCM. − A suitable approach to consider the superficial flow (squeeze flow) is crucial to further exploit the capabilities of this approach, especially with regard to more complex geometries or unfavourable process parameter settings. − Heat conduction and curing needed to be implemented in forthcoming steps.
− The macroscopic, single-phase approach neglects air inside the mould as well as within the microstructure of the reinforcements, this currently prevents a prediction of air-traps and micro-saturation.

Conclusion
In the first part of this work, experimental trials are conduced, providing an insight into the physical mechanisms during WCM. Compaction and permeability behaviour of an undeformed and sheared woven fabric are determined and obtained results are found to be in good agreement with existing literature. Moving on, a novel transparent double dome tool is presented which facilitates in situ visualization of fluid progression during moulding. The results demonstrate the strong interaction between draping and fluid progression. Moreover, draping behaviour proved to be particularly important, because it provides the deformed porous network for the compression-induced infiltration.
In this light, and in contrast to existing form filling models, a forming-driven, three-dimensional WCM process simulation approach is presented containing two coupled forming and fluid-submodels. The forming submodel is able to account for non-linear transversal compaction and material damping without numerical locking. It is verified on part level via comparison with measured reaction forces during moulding and results from literature. Beyond that, the fluid model, based on Darcy's law, is able to account for both the interface flows and anisotropic permeability during large strains and for high element aspect ratios. To enable a more numerical efficient coupling C.T. Poppe et al. between the models, the fluid model is solved within the explicit solver in ABAQUS. Test cases are presented to demonstrate the validity of the fluid submodel. Finally, the combined forming-fluid process model is benchmarked with experimental results via in-mould fluid progression and fluid distribution during moulding. The strengths and drawbacks of the promising results are discussed.
Although demonstrated for woven fabrics, this approach can be adapted to other material classes such as uni-directional non-crimped fabrics (UD-NCF) or nonwovens as we will investigate in future work. Moreover, an explicit representation of the superficial fluid and implementation of curing are reasonable next steps. In addition, a strong fluid-structure-interaction with consideration of the fluid pressure within the draping model will be aspired, which should enable the prediction of more complex process effects such as flow-induced fibre displacement and viscous compaction forces.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.