Momentum conserving dynamic variational approach for the modeling of ﬁber-bending stiﬀness in ﬁber-reinforced composites

: Rotor-dynamical systems made of 3D-ﬁber-reinforced composites which are subjected to dynamical loads exhibit an increased ﬁber bending stiﬀness in numerical simulations. We propose a numerical modeling approach of ﬁber-reinforced composites that treats this behavior accurately. Our model uses a multi-ﬁeld mixed ﬁnite element formulation based on a dynamic variational approach, as demonstrated in [4], to perform long-term dynamic simulations that yield numerical solutions with increased accuracy in eﬃcient CPU-time. We extend a Cauchy continuum with higher-order gradients of the deformation mapping as an independent ﬁeld in the functional formulation, as suggested in [2], to model the bending stiﬀness of ﬁbers accurately. This extended continuum also takes into account the higher-order energy contributions including the ﬁber curvature along with popular proven approaches that avoid the numerical locking eﬀect of the ﬁbers eﬃciently. We apply the proposed approach on a cantilever beam with a hyperelastic, transversely isotropic, polyconvex material behavior in a transient dynamic analysis. The beam is subjected to a bending load with a strong dependence of the overall stiﬀness on the ﬁber orientation. The spatial and temporal convergence as well as the conservation properties are analyzed. It is observed that the model needs an improved numerical treatment to conserve total momenta as well as total energy.


INTRODUCTION
The finite element method for dynamical problems has received much attention over the last two decades, and approaches to solve them are still computationally demanding and timeconsuming. The extensive development of new materials like fiber-reinforced composites constantly creates a need for more generalized algorithms for numerical simulations. The approach to avoid locking behavior in finite elements has significantly improved the accuracy and efficiency of almost any modern finite element code. Nevertheless, the usage of the same is not widely understood in the dynamic regime. Moreover, any modification of the standard continuum to better the accuracy of the numerical solution has to satisfy corresponding physical balance laws. Recent developments in the energy-momentum scheme provide a better opportunity to address these problems in a dynamic regime. The motivation of this work is to enhance the application of the lightweight design of rotating systems made of fiber-reinforced composites by taking advantage of the outcomes of the research mentioned above. In [1], the authors investigated the deformation patterns on mesoscopic level induced by the fibers in fiber-reinforced composites. Their results from the three-point bending test point out that these deformations eventually influence the bending stiffness of the composite material on the macroscopic level. Unfortunately, a standard Cauchy continuum is not well suited to capture Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain Figure 1: Transversely isotropic continuum with fibers oriented in direction of a 0 and element-wise independent deformation gradient (cp. [3]).
these effects in fiber-reinforced models in numerical terms. Various approaches that have been proposed to solve this issue are limited to static problems [2,8]. These drew our attention to capture the fiber bending stiffness in dynamical problems, which would help us reduce the unaccounted out-of-plane bending rigidity of an arbitrary geometry.
As a first step, to capture the fiber bending stiffness, we begin with assuming a constitutive model, where the strain energy function takes not only the strain and fiber direction vector into account, but also the information of fiber-curvature. A transversely isotropic continuum B 0 is considered with fibers at each point of the continuum oriented in direction of vector a 0 in material configuration. In contrast to [2], we introduced a deformation gradient F as an element-wise independent field in our Hu-Washizu based internal energy functional in [7]. Similarly, Γ is introduced as an independent mixed field for ∇ X [ F ] to capture the fiber curvature effects. In this work, we propose an additive split of strain energy function in terms of C and Λ as (see Figure 1), Ψ total (I i ( C, Λ)) = Ψ iso (I 1 ( C), I 2 ( C), I 3 ( C)) + Ψ aniso (I 4 , I 5 ) + Ψ hg (I 6 ( Λ)), which is in line with the variation of [2], where I 6 := k 0 · k 0 , k 0 := ( Λ · a 0 ) and Λ is an independent mixed field for Λ := F t · G which is the pure referential representation of G. G is defined as the referential gradient of the spatial fiber direction vector a t =λ Fā t andλ F is the fiber stretch. Thus, I 1 , I 2 , I 3 , and I 4 , I 5 are the usual isotropic and anisotropic principal invariants based on the right Cauchy green tensor C, which is a mixed field variable for C = [ F ] t · F . For simplicity, anisotropic part of the strain energy due to I 4 , I 5 is assumed to be constant and its effects are neglected within the framework of this work. In this way, we frame our new extended continuum.
Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain

PRINCIPLE OF VIRTUAL POWER
In the second step, we formulate the power functional for the extended continuum. The Equation (3) shows the internal power functional with new independent mixed field variables F and C, energetically conjugated with independent first Piola-Kirchhoff stress tensor P and second Piola-Kirchhoff stress tensor S, respectively. Similarly, Γ and Λ are energetically conjugated with independent B and Λ, respectively.
Here we represent triple contraction of tensors by 3 . The mass-specific body load B and a traction loadT on the Neumann boundary ∂ T B 0 are considered as external forces. Further, algorithmic stress tensorsS andĀ are introduced in the external power functional to derive energy-momentum time integration. More details on this topic can be found in [5].φ denotes the prescribed boundary displacement with respect to the reaction force R as its associated Lagrange multiplier in the Dirichlet boundary ∂ ϕ B 0 . These yield to the following external power functional, The algorithmic stress tensors are defined as, And finally, the kinetic power functional with mass density ρ 0 , velocity v and linear momentum p is defined by, Book of Extended Abstracts of the 6 th ECCOMAS Young Investigators Conference 7 th -9 th July 2021, Valencia, Spain

WEAK FORMULATION
In the next step, to derive a weak formulation for extended continuum, we apply virtual power principle to the total energy balance of the system leading to the following equation, The variation of all power functionals is performed with respect to their dependencies to derive weak forms from the corresponding virtual powers. As in [4], the symbol δ * is used in the sense of variation performed with respect to both temporally continuous time rate fields and temporally discontinuous Lagrange multiplier fields. The resulting integrals of the weak forms of the extended Cauchy-Boltzmann continuum with fiber curvature is expressed in their continuous form in this paper for simplicity. The weak mechanical momentum equation is obtained as, To solve equation (9), the first Piola-Kirchhoff stress is required and determined from its weak form, Similarly, to solve the above equation, independent strains tensors are obtained from their corresponding weak strain equations: and stress tensors from their corresponding weak stress equations: Weak curvature-strain equations are expressed as, Linear momentum needed in the weak mechanical momentum equation (9) can be obtained by dissolving weak velocity equation into weak momentum equation, We discretize the weak forms spatially and temporally on the elemental level by Gaussian quadrature using Lagrangian ansatz funtions. The time rate variable fields (•) e,n i are approximated on the n-th time step by k + 1-th order Lagrange polynomials corresponding to the normalized time α ∈ [0, 1] on each time step [t n , t n+1 ] by and the stress fields as well as Lagrange multiplier fields are (•) e,n i are approximated on the n-th time step by k-th order Lagrange polynomials by Similarly e-th finite element are approximated in space using standard local shape functions N A (ξ), A = 1, · · · , n node defined on the reference domain. The resulting tangent matrix is condensated to pure displacement form by staggering the solution of globally discontinuous mixed fields on the elemental level. We implement this in our In-house finite element code 'fEMcon' and the resulting linear systems of equations are solved using PARDISO solver [6].

BALANCE LAWS
With the introduction of new independent field variables, the extended standard Cauchy continuum has to fulfill physical balance laws. To conserve total momentum and energy at every discrete time step entails doing a special numerical treatment.
Following the steps in [5], suitable test functions δ * φ , δ * ˙ F and δ * P are employed in (9) for an arbitrary axial vector c = constant and eliminating the first Piola-Kirchhoff tensor yields a time-integrator that eventually conserves total angular momentum, Similarly, employing different choice of suitable test functions for δ * φ , δ * ˙ F and δ * P and eliminating the first Piola-Kirchhoff tensor, we derive the total energy conserving time-integrator as below,

NUMERICAL EXAMPLE
In order to understand the anisotropic behavior exhibited by the fibers, we apply our proposed approach on a simple cantilever beam of length 15cm, width 2cm, and height 1cm. The numerical model is assumed to be reinforced with a single family of extensible fibers submerged in the matrix material and exhibiting resistance to bending. 20-noded tri-quadratic serendipity elements are used to discretize the beam into 24 finite elements. A Gaussian quadrature scheme with 27 quadrature points is employed to evaluate the integrals numerically. The left end is chosen as the Dirichlet boundary, such that the displacement of nodes at this boundary are fixed in all three directions e 1 , e 2 , e 3 . As a Neumann boundary condition, on the opposite  free end of the beam, a deformation-dependent transient pressure loadp = 200f is prescribed, which always creates traction in the direction parallel to this surface. Standard Neo-Hookean type material ansatz is chosen for the isotropic part of the strain energy function and cI 6 for the higher-order gradient energy part as in [2], Ψ total (I 1 , I 3 , I 6 ; λ, µ, c) = λ For this setup, the simulation is performed for following test cases: 1. the fiber bending stiffness material parameter c is varied assuming the fibers are oriented with beam's axis, i.e. a 0 = e 1 2. the fiber orientation a 0 is varied for a constant value of the fiber bending stiffness material parameter c The material and simulation parameters chosen as per the table above are in SI units.

Influence of the fiber bending stiffness
For c = 0, the numerical model exhibits the behavior of a non-reinforced beam. As expected, our results show that increasing the stiffness parameter value stiffens the overall response of the composite material. With the increasing values of c, the onset of the higher-gradient part of the energy function is more pronounced. However, fiber stiffness has no significant influence beyond a certain range of c for a chosen load. Figure 4 presents the trend of displacement of a point P at the top edge of the Neumann boundary in e 2 direction.

Influence of the fiber orientation
For c = 10 6 , the orientations of the fiber reinforcements are varied to understand the behavior of the beam. The fiber angle takes values between α = 0 and α = π/2. It is observed in Figure  5 that the cantilever is less stiffer for the fiber orientation a 0 = [1 1 1] t , and as expected, the composite is out of its longitudinal plane. What is surprising is the fact that the degree of stiffening achieved for α = π/6 and α = π/4 is more than for the orientation α = 0, which is counter-intuitive. However, from Figure 5 we can understand that the trend of the plot is independent of inertia. Despite this, we can still state that our proposed time integrator conserves total momenta and total energy.

Different type of energy function
Complementing the above fact on total energy conservation in previous test cases, we also studied the effect of the time integrator with a non-linear energy function based on the quadratic of curvature measure cI 2 6 . It is evident from Figure 6 that for the chosen linear and non-linear type of anisotropic energy functions with respect to higher-order gradient of the deformation gradient, the total energy is conserved.

CONCLUSIONS
To sum up our work, we demonstrated the influence of fiber curvature on the bending stiffness of the cantilever beam as a numerical example. We introduced an independent field variable for spatial fiber direction vector in the continuum using Hu-Washizu's principle to capture the curvature effect. We have succeeded in combining an energy-momentum scheme with the principle of virtual power for the proposed mixed element formulation to preserve the time evolution of energy functions. In this way, the spurious errors arising from fibers are significantly reduced in numerical simulations. In addition to that, our energy-momentum scheme guarantees to obtain the desired accuracy with larger time steps and therefore reduced total CPU time.
The presented contribution has highlighted the importance of the curvature measure through the invariant I 6 in bending-dominated problems in dynamic scenarios. The maximum bending stiffness has been obtained with increased fiber bending stiffness parameter. For the simulated coarse mesh, it is observed that the bending stiffness is maximum when the fibers tend to align with the beam's axis except for smaller angles. To further our research, we intend to extend our formulations to thermomechanical problems.