A mixed-dimensional CutFEM methodology for the simulation of fibre-reinforced composites

We develop a novel unfitted finite element solver for composite materials with quasi-1D fibrous reinforcements. The method belongs to the class of mixed-dimensional non-conforming finite element solvers. The fibres are treated as 1D structural elements that may intersect the mesh of the embedding structure in an arbitrary manner. No meshing of the unidimensional elements is required. Instead, fibre solution fields are described using the trace of the background mesh. A regularised “cut” finite element formulation is carefully designed to ensure that analyses using such non-conforming finite element descriptions are stable. We also design a dedicated primal/dual operator splitting scheme to resolve the coupling between structure and fibrous reinforcements efficiently. The novel computational strategy is applied to the solution of stiff computational models whereby fibrous reinforcements may lose their bond to the embedding material above a certain level of stress. It is shown that the primal-dual 1D/3D CutFEM scheme is convergent and well-behaved in variety of scenarios involving such highly nonlinear structural computations.


Introduction
In this paper, we develop an efficient unfitted finite element solver for composite materials and structures featuring fibrous reinforcements. Such composites are commonly used in industry. In civil engineering, standard concrete structures are reinforced with steel lattices to circumvent the poor ability of concrete to keep its integrity under tensile loading. Short fibres may be also be added to concrete mixes to increase their strength and ductility. To characterise and certify engineering structures employing such reinforced composites, it is important to devise efficient and reliable finite element solvers to simulate the associated micro-mechanics.
Simulations involving "quasi-1D" embedded structural elements are difficult to carry in a standard finite element context. Indeed, representing such slender elements using finite elements that are dimensionality consistent with the embedding space leads to prohibitively large computational models. The 3D analysis of large-scale slender elements may at best be performed locally [23]. Away from regions of interest, asymptotic or numerical cross-sectional approximations of the mechanics of these structures need to be derived for numerical models to be tractable. In the context of quasi-1D embedded fibres, slender elements are typically represented as unidimensional structures with the properties of bars or beams. The mechanical coupling between these lower-dimensional elements and the surrounding domain needs to be handled with care when developing finite element solvers.
Following the taxonomy preposed in [23], three traditional finite element methodologies may be distinguished. In the diffuse approach (see e.g. [28]), stiffness is added to the element of the embedded domain that are intersected by reinforcements, leading to the desired local anisotropy. Such formulations are restricted to the analysis of arrays of slender elements that are geometrically regular and mechanically perfectly bonded to the surrounding material. In the discrete, conforming approach (see e.g. [28]), reinforcements are geometrically represented by edges of the embedding finite element mesh. The advantage is that slender elements may now possess their own kinematics, allowing to simulate sliding for example. However, elements of the embedding domain need to conform to the geometry of the fibrous reinforcement, leading to additional cost, and potentially to meshing difficulties. Finally, modern computational approach may use an embedded representation of fibres, whereby fibres are explicitly represented but may cross elements of the surrounding material in an arbitrary manner [5,17,23,27]. The present methodology belongs to the latter class of methods.
We focus in this paper on the simulation of 2D and 3D structures reinforced by 1D unidimensional structural elements. The reinforcements will be straight and we will assume that they have no bending energy (i.e. bars in tension and compression). In this specific but important engineering context, we propose to develop a novel primal/dual CutFEM approach to model the assembly of composite components. The fibres will possess their own degrees of freedom, as in [10,15]. This fundamental feature will be leveraged to represent sliding behaviours between the fibres and the surrounding material. Consistently with the CutFEM framework [6,18], no specific computational mesh will be created to represent the 1D reinforcements. Instead unidimensional mechanical fields will be described as the trace of intersected embedded elements. Such formulations have already been explored to solve Laplace-Beltrami equations over embedded manifolds [9,10], and to solve fluid flow problems in networks of cracks [10,14,16].
The novelty of the CutFEM approach proposed in the present paper is threefold. Firstly, we extend the "CutFEM on embedded manifolds" premise to the solution of advanced nonlinear solid mechanics problems involving unidimensional reinforcements. Secondly, we propose an efficient primal-dual numerical algorithm based on the LaTIn algorithm [3,21,22] to solve the 1D/3D coupled problem iteratively. This primal-dual formulation is a natural extension of our previous work on primal-dual CutFEM technologies for unilateral contact problems [12] (see [1,8,19,26] for closely related pieces of work from other research groups), and is relatively new in the context of CutFEM approaches, where consistent penalty formulations are usually chosen as a first step to the development of a coupling strategy [7,11,13]. Finally, CutFEM formulations need to be carefully regularised for system matrices to be well conditioned and, in the context of primal-dual algorithms, for inf-sup stability conditions to be satisfied. The regularisation technique proposed in this paper departs from earlier work [10] in that it is devised specifically for primaldual coupling conditions and that it relies on a consistent combination of ghost penalty regularisation [6] and added diffusion in the band of intersected embedding elements [10].
The proposed CutFEM methodology is to be applied to a family of nonlinear problems that are relevant to civil engineering design. We have two application targets in mind. Firstly, we wish to develop a solver for fracture in short-fibre reinforced concrete structures. Inclusions and inclusion/matrix nonlinear interface conditions will be handled efficiently using the primal-dual CutFEM algorithm presented in [11]. The proposed primal-dual algorithm will allow us to represent stiff slender reinforcements using a very similar algorithm, ensuring that the new developments proposed here are fully compatible with our existing computational library. Matrix failure will be represented by using the phase-field model and staggered solution algorithm presented in [25], without modification. Secondly, we wish to develop a methodology to simulate fibre pull-out. To this end, another novel element of the paper is the development of a new formulation for imperfectly bonded 1D fibres in 3D materials. More precisely, sliding will be allowed above a certain threshold, following a thermodynamically consistent approach, yielding a primal-dual stiff friction model (or a model of perfect plasticity with infinite stiffness for the 1D/3D interface) whose solution will be naturally handled by the LaTIn algorithm.
The article is organised as follows. In "Embedded fibre problem" section, we present the problem setting. Straight elastic unidimensional fibres are embedded in a background elastic (potentially damageable) material. The new imperfect bonding formulation is presented next. In "Stabilised enriched finite element formulation" section, we present the primal-dual cut finite element formulation that will be dedicated to the solution of composite materials with fibrous reinforcements. Stabilisation aspects are carefully detailed in this section. The LaTIn iterative scheme is detailed in "Operator-splitting LaTIn algo-

Embedded fibre problem
We consider the problem of an elastic body occupying domain ⊂ R d in dimension d ∈ {2, 3} (see Fig. 1). The material is reinforced by N f straight unidimensional fibres. The ith fibre occupies domain where P i 1 and P i 2 are two elements of R d . As a particular case, P i 1 and P i 2 are the extremities of fibre i if both points belong to domain . We define f = n f i=1 i f and we further assume that i f ∩ j f = ∅ if i = j (i.e. the fibres do not intersect). We define the director vector t i f of the ith fibre as t i f = . The norm used is the euclidean norm. Here, t f will denote the field defined over f whose value is t i f in i f . We also introduce basis n f (respectively (n 1 f , n 2 f )) for the line (respectively the plane) normal to t f , in the case d = 2 (respectively d = 3). The boundary of domain is denoted by ∂ and is additively split into a Dirichlet part ∂ d and a Neumann part ∂ n .
We denote by u m : → R d the displacement field of the matrix of the composite material. Here, u m will be searched in Sobolov space U m = H 1 ( ) of functions whose derivatives up to order one are square integrable. The displacement of the i th fibre is denoted by u i f : i f → R d , and we additionally introduce notation u f : f → R d to denote the field of f whose restriction to subdomain

Elastic composite with perfectly bonded embedded fibres
To ease the reader into the proposed numerical scheme, we first introduce a simple solid mechanics problem whereby the matrix and fibres are linear elastic, and where the bonding between these two phases is perfect (i.e. there is no displacement jump). We therefore seek a solution u := (u m , u f ) ∈ U m × U f minimising the energy functional where the potential energy functionals E m and E f are defined by and respectively. In the previous expressions, ∇ s . = 1 2 ∇ . + (∇ . ) T is the symmetric part of the displacement gradient. ∇ f is the projected gradient operator defined as Field f d : → R d is a volume source term, and t d : ∂ n → R d is a field of Neumann conditions. Finally, k f ∈ R is a fibre stiffness parameter, and C ∈ (R d ) 4 is the fourth-order Hooke tensor defined by its action, which in the case of isotropic materials is In the previous expression I d is the identity tensor, and λ and μ are the Lamé constants. Importantly, due to the action of the projector gradient operator, the fibres have no bending energy. The minimisation of J must be performed under constraints and The last equality enforces perfect bonding between matrix and elastic fibres. The problem of minimisation under constraint can be recast as the unconstrained extremisation of a Lagrangian. The corresponding variational principle is the following Euler-Lagrange system of coupled equations: with the additional condition that Dirichlet condition (7) must be satisfied. We have used notation Z = L 2 ( f ). Here, the Lagrange multiplier field λ ∈ Z can be interpreted as the opposite of a density of forces applied by the embedding material to the 1D fibres (choose δu m = 0 in (9) to see this).
For later use we define the jump operator by [u] := u m − u f , where we implicitly assume the existence of the cartesian product u between u m and u f .

Embedded elastic fibres with perfectly plastic bonding
We now relax perfect bonding condition (8) and only impose it in directions that are normal to the fibre directors, i.e.
which may also be compactly written as . For later use, we split λ as follows: We assume that equilibrium Eq. (9) is satisfied. Due to the relaxation of the bonding conditions, λ f is now undetermined. In the following, we introduce a rate-independent plasticity-like evolution law for this quantity, following the classical framework of thermodynamically consistent materials.
The nonlinear time analysis will be conducted over pseudo-time interval T = [0 T ]. Field u m , u f and λ are now defined over product spaces U m × T , U f × T and Z × T , respectively. For the sake of simplicity, the initial conditions u m ( . , 0) and u f ( . , 0) will be set to zero.
At any time t ∈ T \{0}, the energy dissipation is defined by where the power of the external forces is (15) and the elastic energy of the composite is defined as Replacing the test functions by velocities in (9), the energy dissipation takes the following expression: where λ := t f · λ. We now postulate that solution (u m , u f , λ) satisfies the principle of maximum energy dissipation at all times, subject to the following inequality constraint: where Y is the unique positive scalar parameter of the friction model. The constrained maximisation principle will yield an expression for the evolution of λ , which is ensured by the convexity of f f . Introducing and extremising Lagrangian with respect to λ and Lagrange multiplier γ f , we obtain the usual Karush-Kuhn-Tucker conditions, which are composed of the three relations and of variational equality where f f ( . ) = sign( . ). The last equation yields the following expression for the Lagrange multiplier: We can finally state the system of equations that governs the evolution of the composite with imperfectly bonded fibres. The system is composed of equilibrium Eq. (9) perfect normal bonding Eq. (11) and the system of equations for the plastic flow of the fibre/matrix bond in the direction of the fibre director vectors The last inequality implies, in particular, that the dissipation is always positive, which means that the proposed model satisfies the second principle of thermodynamics. Whilst the two first equations allow us to calculate the amplitude of λ , the last one yields its sign. It is clear that λ acts in opposition to the motion of the matrix relative to the fibres.

Time discretisation
We use the implicit Euler scheme to discretise the time derivatives that appear in the plasticity model. T is split into n t time slabs of equal length T . The resulting discrete time grid is denoted by T t = {t 0 , t 1 , ... , t n t }. At time t n , n > 0, system (23) becomes: Notice that step length T does not appear in this set of rate-independent equations (i.e. pseudo-time dependency).

Cut finite element formulation Finite element spaces
We will solve the composite problem using an enriched finite element formulation. In particular, the embedded fibres will possess their own degrees of freedom defined via the finite element space of the embedding material restricted to the band of intersected elements. Let us introduce a triangulation T h of domain . Here, h is the diameter of the smallest sphere containing the elements of T h . We now define the finite element space U h ⊂ U m by Next, we define the set of all elements of T h that have a non-zero intersection with f (i.e. grey elements in Fig. 2) by The domain corresponding to this set is denoted by G h := K ∈G h K . Importantly, all fields defined over the embedded fibres will be represented by the restriction of finite Of key importance to CutFEM methods [7,12,18] is the introduction of regularisation terms to ensure that system matrices remain well-conditioned, notably by applying some form of ghost-penalty [6] stabilisation. To this end, we define the set of intersected element edges by Finite element formulation At time t n , n > 0, we look for a solution composed of fields Here, u m,h is the finite element displacement in the matrix phase. u f,h is a finite element field defined over the band of elements intersected by the fibers, and whose restriction to the 1D fiber domains will be our finite element approximation of the fiber displacements. λ h is the extended approximation of the Lagrange multiplier field, which will be sought for in the same space as u f,h . Fieldsw h andλ h are additional variables that need to be introduced to stabilise the CutFEM formulation.λ h is a fiber Lagrange multiplier field that belongs in L 2 ( f ). Finally,w h is a field that represents the jump of displacement between the fibers and the surrounding matrix material. This field will also be looked for in L 2 ( f ). Fieldsw h and λ h are not discretised a priori. Their numerical values will only be required at quadrature points defined over interface f (through the discrete evaluation of integrals). The 5 fields introduced here are fully defined by the set of equations provided below.
To avoid cluttered symbolic expressions, we will omit the dependency of solution fields to time. The only exception is for quantities at time t n−1 , which will be indicated via the superscript. Other notations, notably the use of superscript and jump symbol [ ] are meant to be consistent with the developments of the previous sections, and should be self-explanatory.
The enriched finite element formulation of the embedded fibre problem employs a twoscale formalism that requires splitting Z additively into the restriction of W h onto f and its supplementary space. Accordingly, the stabilised and extended L 2 -projection λ h ∈ W h of fieldλ h ∈ Z is defined by Here, term s ♥,h is a regularisation term that will be discussed in the next section.
This variational equality is complemented by the Signorini-type law for the. quantities locally in f , and by d − 1 scalar equality constraints Last, the following closure equation is introduced, following [12] The last equation regularises the. interface quantities by penalising the distance to their . h bulk counterpart. Parameter γ must be strictly positive. The bilinear and linear forms used in (30) are respectively defined by and Here, term s ,h is the second regularisation term of the formulation. It is discussed next.

Stabilisation
Without stabilisation (i.e. definition of s ,h , s ♥,h and mixed condition (33)), the 1D-3D CutFEM formulation suffers from several sources of numerical instability.
1. Elements may have a small intersection with a particular fibre, leading to a loss of control of finite element field u f,h (fibre displacement, which is extended to the band of intersected elements) and finite element field λ f,h . 2. The extended fibre kinematics (displacements and Lagrange multiplers) allows zero energy modes to be represented, leading to singular systems of equations. 3. Uncontrolled oscillatory modes may appear in the Lagrange multiplier fields due to the non-satisfaction of the LBB condition.
These problems will be circumvented via the introduction of the following regularisation term: where n F denotes the facet normal and the orthogonal gradient operator is defined by Ghost penalty parameter β is typically set to 0.1. The corresponding term extends the ellipticity provided by the fibres' elasticity operator to the entire band of intersected elements, thereby alleviating the ill-conditioning due to uncontrolled modes of u f,h · t f . Let us emphasise that only the tangential part of u f,h is controlled by the ghost-penalty stabilisation. This is consistent with the formulation of the deformation energy of the fibres, which develops in reaction to the projected displacement gradient. The normal part of the fibre displacement does not need to be regularised. It will simply follow the bulk displacement, which is trivially enforceable, the fibre having no bending energy to act against the corresponding constraint. This addresses point 1 above.
The second regularisation term eliminates zero energy modes from the system, which addresses point 2 above. It is a diffusion term which is integrated over the entire band of element, although a facet penalisation could have been used as well [10]. This term does not introduce any inconsistency (an extension of the exact solution in the normal direction may be constructed, which has vanishing orthogonal gradient). Typically, a small parameter value ζ = 10 −5 is sufficient to stabilise the formulation.
A similar term is proposed to stabilise the Lagrange multiplier field (see Points 1 and 3 above): where κ = is the harmonic mean of the fibre stiffness and the equivalent axial stiffness provided by the matrix. L is a non-dimensionalising parameter that represents the width of the matrix region that reacts to the deformations of a particular fibre. As shown in [12], the first term helps control the spurious oscillations that may appear in the Lagrange multiplier field. The second term eliminates zero energy modes. We showed in [12] that this type of stabilised primal/dual formulation could be interpreted as an alternative formulation for the stabilised augmented Lagrangian proposed in [8], in the spirit of the classical approach for unilateral contact problems first introduced in [2].

Overview
We will develop a LaTIn solver that algorithmically decouples the matrix and fibre finite element problems. LaTIn solvers have proved to be particularly efficient when solving stiff Signorini-type problems without the need to further smoothen or regularise the discrete system of governing equations. The To close the system of equations, we require the iterate to satisfy the descent search direction given by The differential equations of the matrix and fibres are decoupled by this algorithmic construction. Positive scalar γ was first introduced in Eq. (33). It is such that γ =γ κL −1 , whereγ is a non-dimensional parameter.
Notice that because of closure Eq. (33) and extended projection property (29), the ascent search direction and Eq. (42) are also satisfied by the. k+ 1 2 quantities.

Linear stage: uncoupled systems of equations
The linear stage defined above consists in solving the following set of linear equations, which are independent for each phase of the composite: As detailed previously, the set of equations corresponding to the fibres are regularised by term s ,h .

"Inf-Sup-regularised" local stage: algorithmic solution
The regularised local stage defined previously can be solved by applying the following algorithmic procedure.
• Set In our example, we use the quadrature point corresponding to the integral of the product of two piecewise linear fields over f . This introduces a quadrature error, whose effect was studied and discussed in [12].

Numerical examples
All the numerical simulations presented in this section are performed using the LaTIn-CutFEM library [12] developed in finite element package FEniCS [4].

2D structure with a hole and two embedded 1D reinforcements
The first numerical example that we wish to discuss is represented in Fig. 3. The structure is two-dimensional, and undergoes plane strain deformations. The computational domain is the unit square = [0 1] × [0 1]. A vertical load is applied in the form of a Dirichlet field u d at the top of the structure, i.e. ∂ d = {x ∈ | x · e 2 = 1}. The Young's modulus E in the circle of diameter r = 0.15 with centre (0.5 0.5) is set to 1.0. Outside the circle, the Young's modulus is set to 0.1. Hence, the circular inclusion is made of a comparatively hard material. Poisson's ratio for both materials are set to 0.3. Two straight reinforcements are introduced, as represented in Fig. 3. The reinforcements are stiff compared to the matrix, with k f = 1.0. The friction coefficient Y is relatively low at Y = 0.2, which will result in early yielding.
The prescribed displacement field is uniform over ∂ d , and of amplitude A(t) as reported in Fig. 5 (top subfigure). The structure is loaded in traction, then unloaded and forced into a compressive state, before being reloaded in traction. Figure 5 also shows the 26 time steps that will be used by the implicit time integrator. There is no source term, i.e. f d = 0 and no imposed tractions, i.e. t d = 0.
A coarse computational mesh is represented in Fig. 3. The set of intersected cells G h over which the finite element fibre displacements and Lagrange multipliers are defined are also represented. The fibres do not conform with the regular background mesh. The circle does not conform to mesh edges either. In order to handle the mismatch in material properties, kinematic enrichment is performed using the CutFEM technology, including ghost penalty stabilisation and the primal/dual coupling introduced in [12].
We now examine Fig. 4. The figure presents the field of Lagrange multipliers after the third time step (framed in the top subfigure of Fig. 5). The colour scale indicates the magnitude of the Lagrange multipliers, which are piecewise linear in the band of intersected elements. We emphasise the fact the all our fibre fields are traces of d-dimensional standard finite element fields defined over these bands of intersected elements. Fibres do not possess their own meshes. We see that the multipliers are aligned with the director vec-   tors. This is normal, but not enforced directly in our approach. Fibres do not have bending energy. Therefore, they can only generate forces in the direction of their respective director vector, thereby resisting to deformations of the matrix phase in these directions only. In our experience, erroneous stabilisation of the fibre fields generate tangential forces that pollute the overall finite element results. In Fig. 5, we report the evolution of these forces with time. Towards the end of the increasing phase of the loading amplitude A(t), these forces saturate. Indeed, they are constrained not to exceed friction coefficient Y . Instead, plastic slip develops to accommodate the boundary conditions. At the end of the compression phase, the forces are also fully saturated, but take the opposite signs, as they now react to compressive loading. The last picture is interesting. It is representative of transition phases that appear when the direction of the load changes. In that case, forces at the extremities of the fibres immediately change directions, while due to adhesion, the forces in the middle of the fibres carry residual forces in the opposite direction, creating a mixed traction/compression stress state in the matrix as well. Notice that the degrees of freedom corresponding to the Lagrange multipliers are well-behaved throughout this analysis.
A sequence of computational meshes is represented in Fig. 6. The solution is for the 20th of the computational times shown in Fig. 5. The structure is globally in tension, but residual compressive stresses exist in the matrix due to the past global compression states. Looking at the finest of the finite element solutions, we see that stress concentrations develop at the extremities of the fibers. Away from these extremities, the stiff fibres restrict the value of the tensile stress that develops in the embedding material. It is clear that the stress field shown in Fig. 6 converges with mesh refinement.
In Fig. 7, the four meshes are warped using the displacement of the matrix phase, and both fibres are warped with their own displacement fields. The solution inside the circle is not represented. The absence of noticeable difference between the 4 displacement fields indicates that the method converges with mesh refinement. The symmetries observed in the displacement of the fibres are also qualitative indicators that the proposed strategy and Fig. 8 Global reaction force of the structure as a function of the prescribed extension. The hysteretic behaviour is due to the friction law. Convergence of this response curve can be seen when the mesh is progressively refined. The meshes used to produce these four curves are represented in Fig. 7 its numerical implementation are correct. Notice that meshes now appear to conform to the geometry of the circular inclusion. These are not computational meshes but auxiliary meshes obtained by subdividing elements that are cut by the inner circle, which is done for integration and visualisation purposes only. The finite element meshes used to describe the displacement field in the matrix are the regular grids displayed in Fig. 6.
In Fig. 8, we present quantitative mesh convergence results for the example described above. We report the evolution of the reaction force at the Dirichlet boundary as a function of the loading amplitude. In the previous expression, n = e 2 is the outer normal to the Dirichlet domain boundary, and e 2 is the second canonical vector of R 2 . Time response R h (t) is computed four times using the sequence of meshes represented in Fig. 7. Clearly, the nonlinear behaviour of the interface between fibres and matrix creates an hysteresis that is akin to that of plastic materials undergoing cycling deformations. The area bounded by the upper and lower part of the cycle is the energy that has been dissipated during the cycle by friction. The response curve converges with mesh refinement. They are also free to rotate. The embedded interface between the four smaller parallelepipeds and the large beam is unilateral frictionless contact. When A(t) is negative, the top two smaller parallelepipeds loose contact with the beam, and the structure is at rest

Self-healing concrete block with temperature-activated elastic fibres
Description of the computational model One of the motivating elements for the present work is the simulation of advanced concrete structures. In recent years, considerable work has been carried out in the development of self-healing cementitious materials. A particular application, currently investigated at Cardiff University, is the use of shrinkable fibres in the development of a crack closure system (see Fig. 9). The fibres are made of polyethylene terephthalate (PET), a shape memory polymer (SMP) that shrinks when heated at temperatures above 60 Celsius [20]. Under restrained conditions and upon thermal activation these SMP fibres can develop significant shrinkage stresses [20,24]. Preliminary results from an experimental study showed that SMP fibres with end-anchorages and embedded in a cementitious matrix were effective in closing macro-cracks [24]. Nevertheless, to achieve an optimum mix design for these advanced cement based composites, further work focused on the optimisation of key parameters (e.g. fibre geometry, fibre content) is required. To this end, virtual experiments carried out with numerical models that are able to capture the underlying mechanisms would be very advantageous. We now propose to show how experimental devices such as the one represent in Fig. 9 may be simulated using the proposed embedded FE technology, which we hope will eventually enable experimentalists to explore material design spaces virtually. The computational model is represented in Fig. 10. It is composed of one main square parallelepiped (the beam), and four additional smaller square parallelepiped designed to replicate the loading conditions of typical testing devices used in civil engineering. We apply zero Dirichlet conditions along a line of each of the two smaller parallelepipeds located underneath the structure. This allows these two support to rotate freely (the mesh conforms with those lines). The two smaller parallelepipeds at the top of the structure are loaded in a similar fashion, but with a downward prescribed motion of the action lines, with an amplitude that is proportional to A(t), i.e. we reuse the loading function introduced in the previous example to simplify the presentation. The interface between the large beam and the four smaller parallelepipeds is represented by a unilateral frictionless contact model. We finally introduce two fibres as represented in Fig. 10, the interface between the fibres and the matrix being modelled by the friction-type model derived previously.
The object composed of the 5 parallelepipeds is meshed conformingly to its boundary. The fibres are not meshed by the meshing software. Instead, they are introduced in a non-conforming manner using the technique proposed in this paper. In a similar fashion, we handle the four unilateral contact conditions described above nonconformingly using the primal/dual LaTIn-CutFEM algorithm [12], i.e. the four contact interfaces do not conform with the facets of the mesh and are represented by the zero isoline of a finite element level-set function.
Numerical results The displacement of the system at the end of the first loading phase is represented in Fig. 10. The set of elements that are intersected by the fibres is also displayed. We see that the fibres and the band of elements they are embedded in follow the global deformation of the structure, as expected.
The system presented here will not behave like the 2D system presented in the previous section. When A(t) becomes negative at time t = 1.5, the two parallelepipeds onto which the prescribed motion is applied move upward without adhering to the remainder of the structure (traction forces cannot be transmitted by the unilateral contact interfaces). The beam will stay at rest until t = 2.7, after which it will undergo bending deformations again as induced by the downward motion of the two parallelepipeds. This resting time will allow us to study residual stresses, i.e stresses that have developed due to inelastic phenomena (i.e. friction), and do not go back to zero once the structure is fully unloaded.
These residual stresses are shown in Fig. 11. We see that when the structure is undergoing progressively increasing bending conditions, the bottom fibre is acting against the extension of the beam, while the top fibre acts against compression. During rest (subplot (b) in Fig. 11), due to friction, these stress states do not subside completely, and a small amount of permanent bending deformations remains without external action.
Effect of the numerical stabilisers In Fig. 12, we show the effect of the stabilisers introduced to make sure that the solution of finite element problems with nonconforming 1D embedded elastic elements is well-posed. The top picture is the representation of an uncontrolled deformation mode that appeared when setting ζ to a value that is very close to zero. The Fig. 11 Lagrange multipliers of the fibre friction model: a development of fibre forces in reaction to bending deformations, limited in amplitude by plastic slip and b residual stresses at rest, once the two top paralellepipeds have lost contact with the beam struture. Stress σ 11 is the longitudinal stress measured along the horizontal axis of the figure zero energy modes explode in amplitude. The restriction of the solution to the fibre is still correct, but carrying such large values of the degrees of freedom through long time analyses is of course not advisable. Setting ζ to a small but non-vanishing value, as suggested previously, ensures that such uncontrolled sets of degrees of freedom are eliminated from the system. This results in fully controlled deformations of bands of intersected elements, as shown in Fig. 10 (last two subplots).
The same remarks can be made regarding Lagrange multiplier fields. When setting ζ λ to a vanishingly small value, uncontrolled modes develop, as seen in Fig. 12a, b. A small regularisation parameter ensures that the problem is well-posed, which is exemplified by the beautiful purely axial distribution of the extended Lagrange multiplier field shown Fig. 12c, d.

Fracture of quasi-brittle composites with randomly distributed short fibres
To further illustrate the capabilities of the proposed approach, we extend it to the simulation of damage in cement-based composite materials reinforced with randomly distributed Fig. 12 Effect of stabilisation terms on extended fibre fields: a uncontrolled extended displacement mode, b uncontrolled extended Lagrange multiplier modes, c regularised extended Lagrange multiplier and d restriction of the regularised extended Lagrange multiplier field to the 1D fibre domain short fibres. The overall mechanical behaviour as well as the crack propagation in these composites are largely gouverned by the crack-bridging action of fibres.
The problem setting is graphically described in Fig. 13. The interface between fibres and matrix is now rigid (i.e. infinite friction coefficient Y ). Uniform linear Dirichlet boundary conditions in the e 2 direction are now applied on the entire boundary ∂ . We continue using function A(t) to described the overall evolution of the load with time. The perfect kinematic coupling between hard elastic inclusions and matrix is handled via the pri- Further stability and convergence tests would be required to study the effect of the damage field on the overall stability of the unfitted finite element formulation. This is left to future investigations.
We note that a particular drawback of the proposed unfitted finite element approach is that fictitious domains corresponding to disconnected physical elements (i.e. two separate inclusions, two different fibres, a fibre and an inclusion), should not overlap for the numerical results to be sensible. This sets a lower bound on the element size, which may preclude their use in preliminary engineering studies where very coarse meshes are used to inexpensively study the overall behaviour of a complex system.

Conclusion
We have presented a stabilised cut finite element solver to simulate the micro-mechanics of fibre-reinforced composites. The fibre elements have their own degrees of freedom, allowing complex interface conditions and efficient iterative solvers to be devised. No specific mesh is required for the slender structural elements. The unidimensional mechanical fields are represented as 1D traces of background 2D or 3D finite element fields.
Starting from this elegant numerical strategy, we developed a primal-dual finite element formulation to simulate friction-restricted sliding between the fibres and the embedding material. We also designed an efficient LaTIn algorithm to solve the resulting stiff nonlinear system of equations.
Just like any CutFEM solver, ours needs to be stabilised to fully eliminate numerical instabilities resulting from arbitrary cutting and embedding of finite element meshes. To this end, we have proposed a natural combination of ghost-penalty regularisation to stabilise the projected gradient, and background mesh penalisation to eliminate uncontrolled gradients in the normal direction. This approach has also been used to stabilise the dual fields. We have shown that the method is stable, convergent, and versatile. It handles high levels of nonlinearities very well, which was shown in examples of featuring 1D/3D imperfect bonding and nonlinear damage mechanics.
This work is an important stepping stone towards the development of flexible virtual testing platforms, in particular for cement-based composite materials. Avenues for further research and developments are manifold. For instance, we would like to investigate how the proposed method can handle structural elements with bending energy, where locking is very likely to cause numerical issue. Similarly, the present investigations are limited to straight fibres, and an obvious extension would be to allow for curved fibres to be represented. Finally, we would like to be able to simulate pull-out, i.e. the reduction of added stiffness as fibres are progressively loosing contact areas with embedding materials when cracks progressively open. This requires to move away from the assumption of small perturbation, and to continuously update the geometrical configuration of the embedded elements in the evolving embedding structure.

A Unregularised LaTIn local stage
We can now check whether the previously made hypothesis is consistent with equation and proceed to the smoothing step of the local stage. Otherwise, our hypothesis is not true: there is sliding, which means that necessarily, We now make the hypothesis that the sign of the relative velocity between fibre and matrix is positive. In this case, we have that