A finite deformation isogeometric finite element approach to fibre-reinforced composites with fibre bending stiffness

It is a common technique in many fields of engineering to reinforce materials with certain types of fibres in order to enhance the mechanical properties of the overall material. Specific simulation methods help to predict the behaviour of these composites in advance. In this regard, a widely established approach is the incorporation of the fibre direction vector as an additional argument of the energy function in order to capture the specific material properties in the fibre direction. While this model represents the transverse isotropy of a material, it cannot capture effects that result from a bending of the fibres and does not include any length scale that might allow the simulation of size effects. In this contribution, an enhanced approach is considered which relies on the introduction of higher-gradient contributions of the deformation map in the stored energy density function and which eventually allows accounting for fibre bending stiffness in simulations. The respective gradient fields are approximated by NURBS basis functions within an isogeometric finite element framework by taking advantage of their characteristic continuity properties. The isogeometric finite element approach that is presented in this contribution for fibre-reinforced composites with fibre bending stiffness accounts for finite deformations. It is shown that the proposed method is in accordance with semi-analytical solutions for a representative boundary value problem. In an additional example it is observed that the initial fibre orientation and the particular bending stiffness of the fibres influence the deformation as well as the stress response of the material.


Introduction
In many industrial applications, fibre-reinforced composites are in high demand because of their advantageous properties compared to non-reinforced materials. Specifically speaking, a reinforcement by certain types of fibres can enhance the stiffness, durability and strength-to-weight ratio of materials to name only a few benefits, cf. [1,2]. Due to the high popularity of fibre-reinforced composites, efficient simulation methods The transposition of a second-order tensor T is denoted as Similarly, two types of transpositions are employed for a fourth-order tensor S, namely S t = S jikl e i ⊗ e j ⊗ e k ⊗ e l and S T = S kli j e i ⊗ e j ⊗ e k ⊗ e l .
For the differentiation of tensorial quantities, the gradient operator with respect to the (spatial) coordinates x, applied to a tensor T of arbitrary order, is introduced as On this basis, the divergence operator is obtained by the double contraction of the gradient with the second-order identity tensor I such that The curl of a first-order tensor a is defined as with denoting the third-order Levi-Civita tensor.

Fibre-reinforced composites possessing fibre bending stiffness
The enhanced approach for the modelling of fibre-reinforced composites, which takes not only the fibre direction but also a fibre bending stiffness into account by introducing an extended form of the stored energy function, has been presented in [9]. Along with new material parameters, additional kinematic quantities in terms of the fibre Fig. 1 Domain B in the reference configuration (B 0 ) and current configuration (B t ) assuming a couple-stress theory. Employed vectors and tensors: placement function ϕ, deformation gradient F, referential outward surface normal unit vector N, current outward surface normal unit vector n, referential fibre direction vector a 0 , current fibre direction vector a t , body force vector b, body couple vector c, vector of tractions t 1 acting on the surface ∂B t , vector of couples t 2 acting on the surface ∂B t , vector of displacements u prescribed on the surface ∂B t stretch gradient as well as the fibre curvature are considered, and a generalised continuum theory which includes higher-gradient contributions in the balance equations is assumed.

Finite strain kinematics
For a body B undergoing finite deformations, two configurations may be considered, namely a reference configuration at time t 0 as well as the current configuration which inherits the deformation of the body at time t. As shown in Fig. 1, the non-linear function ϕ : B 0 → B t describes the mapping of the position of a material point X of the reference configuration to its position x in the current configuration. The gradient of the mapping with respect to the reference position X is the deformation gradient F = ∇ X ϕ with J = det(F) > 0. The deformation gradient describes a linear relation between line elements dx in the current configuration and dX in the reference configuration, while its cofactor and determinant relate infinitesimal area as well as infinitesimal volume elements in a similar manner, i.e. dx = F · dX, n da = cof(F) · N dA, dv = J dV . (12) Vectors N and n denote the outward surface normal unit vectors in the reference and in the current configuration and the cofactor is defined as cof(F) = J F −t . For the specific modelling approach examined in this contribution and presented in [9], the second gradient of the placement function is considered further, so that is introduced. Within the scope of modelling fibre-reinforced composites, a unit vector field a 0 is commonly incorporated into the model as an additional directional quantity. It describes the fibre direction in the reference configuration for a material which is reinforced by a single family of fibres. Under the assumption that the fibres are embedded in the matrix material and convected with the deformation, the fibre direction vector in the deformed configuration is obtained as Since not only the direction but also the length of the fibres changes consequently to the deformation, the fibre stretch λ a occurs in (14).
Assuming that the embedded fibres are not perfectly flexible but instead exhibit a certain bending stiffness, the gradient of the fibre direction vector is incorporated into the enhanced material model. With respect to the reference configuration, this gradient takes the form cf. [9]. As discussed in detail in [20], the stress and energy contributions related to this gradient term do, in general, not vanish in the initial state. In order to obtain an initially stress-and energy-free condition, the fibres are assumed to be straight in the reference configuration, i.e. ∇ X a 0 = 0.

Local balance equations
For the consideration of higher gradients of the placement function, an enhanced continuum model is required. In the present contribution, the couple-stress theory presented in [14] is applied. The spatial versions of the local balance equations will be summarised in the following for the quasi-static case and under the assumption of vanishing body forces and body couples, accordingly b = 0 and c = 0, cf. [9,20]. Considering a closed system, the local form of the balance of mass readṡ with spatial mass density ρ and with the material time derivative of • denoted as•.
The local form of the balance of linear momentum follows as with an, in general, non-symmetric Cauchy-type stress tensor σ . From the local form of the balance of angular momentum, its skew-symmetric part can be directly obtained as a function of the divergence of the couple-stress tensor m, to be specific By using the balance equations (17a) and (18a), the local form of the balance of energy follows as where U denotes the (mass specific) internal energy. For the isogeometric finite element analysis, the balance equations of linear and angular momentum are combined into one partial differential equation which is of fourth order since the couple-stress tensor, in general, may contain second-order derivatives of the placement function. Retaining the assumptions mentioned above, inserting (18b) into (17b) yields As discussed in detail in [14] and typical for the couple-stress theory, only the deviatoric part m = m − 1 3 tr(m) I of the couple-stress tensor contributes to the energy balance (19) because I : ∇ x ∇ x ×φ = 0 holds. The same applies to the partial differential equation (20), where the spherical part of the couple-stress tensor similarly vanishes as a result of the curl and divergence operation. Consequently, the spherical part of the couple-stress tensor does not have any impact on the balance equations relevant for the proposed simulation approach and accordingly remains undetermined as this work proceeds.

Constitutive model
The constitutive model for fibre-reinforced composites including fibre bending stiffness has been derived in [9]. The basis for this model is an extended list of invariants for the stored energy density function. Specifically speaking, the energy is considered as an isotropic function of invariants I i including three main arguments that have been introduced in Sect. 2.1, i.e.
This purely referential representation is based on the Cauchy-Green tensor C = F t · F, the referential fibre direction vector a 0 and the gradient Λ of the deformed fibre vector. In order to reduce the number of invariants, additional assumptions are made in [9] and adopted in this work. By employing only the directional projection κ 0 = Λ · a 0 , effects from fibre splay are neglected in addition to effects from fibre twist. Apart from that, the sense of the fibre orientation is not relevant from a physical point of view so that the fibre vector may appear in the stored energy only in even powers. From the general form of the stored energy density function in (21), the stress and couple-stress tensor are derived in [9] by considering the particular dependencies of the invariants. Accordingly, the symmetric part of the stress tensor follows as and the deviatoric part of the couple-stress tensor takes the form The specific form of the energy function used for the determination of the partial derivatives in (22) and (23) in this work will be presented in Sect. 4.

Isogeometric finite element approach
For the solution of the fourth-order partial differential equation in the simulation of fibre-reinforced composites with fibre bending stiffness, the isogeometric finite element method is employed. Since NURBS basis functions, used within the isogeometric analysis, provide C p−1 -continuity everywhere, except for the locations of repeated knots or control points, global continuity higher than C 0 can be realised. The NURBS basis functions are used in the discretised weak form of the balance equation which will be developed in the following section. For the application of the finite element method, a linearisation is performed afterwards so that the globally assembled system of equations can be solved by means of the Newton-Raphson method. Detailed information on the isogeometric analysis and on NURBS is provided in [27][28][29].

Discretised weak form
With regard to a finite element formulation of the governing equation (20), a weak form is derived. For this purpose, (17b) is multiplied by a test function η and integrated over the spatial domain B t , accordingly The application of the divergence theorem and of integration by parts yields Considering an additive decomposition of the stress tensor in the form of σ = σ sym + σ skw (26) and making use of the definition of the skew-symmetric stress part from (18b) with σ skw = − σ skw t , (25) can be rewritten as Therein, Cauchy's theorem is employed so that the traction vector t 1 = σ t · n is introduced. For the last term in (27), integration by parts as well as the divergence theorem are applied a second time so that is obtained. Vector t 2 = m t · n represents the couples acting on the surface such as the previously introduced vector t 1 accounts for the surface tractions. The terms which are related to the couple-stress tensor in (28) can be rewritten by using the definition of the curl operator, so that an alternative representation of the weak form reads This format allows an interpretation of the specific occurrences of the test function from the perspective of the principle of virtual power. The test function η itself may be regarded as a virtual velocity field which is related to the classic tensorial quantities in (29). In analogy, the term 1 2 ∇ x × η, which corresponds to the higher-gradient terms in (29), may be interpreted as a virtual spin vector, cf. [20].
From the weak form (28) that has been derived for the deformed configuration, the respective formulation in the reference configuration is obtained by making use of the relations specified in (12). From the third term on the right-hand side in (28), which includes the second-order gradient of the test function, two separate terms follow from the pull-back operation. The weak form in the reference configuration accordingly reads Following the isoparametric concept and employing a Bubnov-Galerkin interpolation scheme, domain B 0 , test function η and placement function ϕ are approximated by the same basis functions. Within the scope of this work, NURBS basis functions R are employed so that the discretised kinematic quantities are with the number of active basis functions on one element e denoted as n en . Inserting both relations into the referential weak form (30), a representation is obtained which includes derivatives of the basis functions up to second order. Accordingly, the internal force vector takes the discretised form where n el denotes the total number of elements. The external force vector respectively is Operator A establishes an assembly of the local force contributions to the global degrees of freedom.
From the integrability condition of the weak form, it follows that global C 0 -continuity is not sufficient for a proper approximation of the field variable and test function, since, e.g., products of second-order gradients occur in the second integral in (32). Within the examples presented in this work, an at least C 1 -continuous NURBS-based interpolation is employed throughout the whole domain.

Linearisation
The solution of the equation system (30) within the isogeometric finite element framework is obtained by employing a Newton-Raphson scheme. The residuum is determined by the internal and external force vectors discussed in detail in Sect. 3.1, i.e.
In the kth iteration step, the linearised form of the residuum reads Due to the particular dependencies of the residual function r h = r h (F h ( ϕ), Υ h ( ϕ)), two contributions are considered for the increment Δr h , to be specific with the global list of degrees of freedom ϕ, the tangent stiffness matrix K and with the partial derivatives in the discretised form. Under the assumption of dead loads, only the sensitivities of the internal force contributions need to be considered. The global system of equations is finally obtained as The general form of the tangent stiffness matrix K is derived in Appendix A.

Representative numerical examples
Two numerical examples are presented in this section to demonstrate the behaviour of a material reinforced by fibres possessing bending stiffness. The first example is used for a comparison of the numerical results obtained by the proposed isogeometric approach against the semi-analytical solution provided in [25]. The second model contains a more complex geometry and is more closely related to practical applications. For both examples, a plane strain condition is employed and the constitutive equations are based on the same energy function.

Specific form of the energy function and stress tensors
In (21), the general form of the stored energy density function for the material model introduced in [9] has been provided. For the particular examples presented in this contribution, the energy takes an additive decomposition and, as in [20], specifically consists of three parts, namely The isotropic part resembles the behaviour of the matrix material without taking into account the reinforcement by fibres. It takes the form The Lamé-type constants are set to λ = 1.037 × 10 5 N mm −2 and μ = 4.4444 × 10 4 N mm −2 . The second part W λ a corresponds to the classic fibre stretch-related transversely isotropic contributions of the fibre-reinforced material.
It has already been studied extensively, e.g. in [7,30,31], and is thus not in the focus of the present contribution. As this work proceeds, the part W λ a is accordingly assumed to be constant. The last contribution incorporates the energy which is related to the higher-gradient terms that are specific for the material model proposed in [9]. From the particular form of the referential gradient of the fibre direction vector (15), it follows that invariant I 6 in (41) includes the gradient of the fibre stretch as well as the fibre curvature, see the discussion in [19]. In accordance with the latter, parameter c is associated with a bending stiffness of the fibres. However, the contributions of the fibre stretch gradient are not to be neglected when finite deformations are considered. The sense of the fibre orientation does not have any impact on the energy contribution in (41) since invariant I 6 is of even order in the fibre direction vector a 0 .
In the examples presented in this contribution, parameter c as well as the fibre direction will take different values in order to examine their influence on the simulation results.
Considering the stored energy density function specified in (39) and using the general derivations of the symmetric part of the stress tensor and the couple-stress tensor (22) and (23), their specific forms follow as and Recapitulating the structure of the stored energy density function (39), the first term in the above given form of the stress tensor represents the isotropic part, whereas the second term and the contributions of the couple-stress tensor correspond to the part that incorporates the fibre properties. For a vanishing fibre bending stiffness, i.e. c = 0, the model reduces to a classic neo-Hookean material.
In the small strain version of the model that has been discussed in detail in [9,22], the contributions from the fibre stretch gradient are neglected in the derivation of the stress and couple-stress tensor because they are of higher order. In contrast thereto, both contributions, namely the fibre curvature and the fibre stretch gradient, are included in the form of the stress tensors in (42) and (43).
In order to obtain the specific form of the tangent stiffness matrix for the considered energy function, the sensitivities of the stress and couple-stress tensor are derived and presented in Appendix A.

Cylindrical tube subject to a pure azimuthal shear deformation
The pure azimuthal shear deformation of a cylindrical tube with radially aligned fibres is considered in order to compare the numerical results with the semi-analytical solution provided in [25]. In [22], the same problem has been analysed by means of an isogeometric finite element method for a small strain setting. The model including finite deformations has been elaborated in [21] by employing a multi-field method.
For an accurate comparison of the results, dimensionless versions of the stress and couple-stress tensor are used, to be specific By using a cylindrical base system {e r , e ϕ , e z }, u ϕ represents the prescribed azimuthal displacement on the surface of a tube with inner radius r i and outer radius r o , see Fig. 2. The particular values for these quantities will be determined in Sect. 4.2.1.

Isogeometric model
The isogeometric model for the cylindrical tube with radii r i = 1 mm and r o = 2.5 mm is adopted from [22]. A polynomial degree of p = 4 is chosen. The employment of linear constraints for the location as well as displacement of repeated control points leads to global C 1 -continuity which is sufficient for the fulfilment of the integrability condition posed by (32). More details on the construction of the cylindrical tube model can be found in [22] and references cited therein.
In accordance with the calculations in [25], the tube is analysed under a pure azimuthal shear deformation which is prescribed on the outer radius by means of Dirichlet boundary conditions, whereas the inner radius is fixed. The azimuthal displacement is chosen in such a way that the outer surface undergoes a rotation of π/9 around the tube's middle axis. The outer radius of the tube is not changed by this deformation so that the applied boundary conditions yield the total volume of the tube to be conserved. Along with this observation, extensible fibres are assumed in order to obtain pure azimuthal shear in the finite deformation setting, cf. [26]. Surface tractions and surface couples are assumed to vanish at the outer cylinder radius, i.e. t 1 = 0 and t 2 = 0. The discretised model including the particular boundary conditions is shown in Fig. 2. The mesh contains n el = 832 elements and is based on n cp = 1260 control points.

Numerical results
In the present example, the fibre bending stiffness is represented by the non-dimensional parameter cf. [22,25]. Within the simulations performed by means of the proposed isogeometric approach, different values have been employed for this parameter, namely λ * ∈ {0.0, 0.005, 0.03, 0.1}. In all cases, the fibres are initially aligned in radial direction. In Fig. 3, the deformation pattern of the tube is shown. In the case where λ * = 0, no fibres are present in the material. As the energy contributions belonging to the fibres become active by employing non-zero values for the parameter, a resistance against bending, which increases in accordance with higher values of λ * , can be observed.
In Fig. 4a, the results for the dimensionless stress contributions [σ r ϕ ] * and [σ ϕr ] * are presented as a function of the cylinder radius. Fig. 4b shows the respective results for the symmetric stress part [σ sym r ϕ ] * as well as for the only nonzero couple-stress contribution [m r z ] * . The presented solutions for the finite strain setting are not only in accordance with the results from [21] produced by a multi-field method by using the same constitutive model, but also with those obtained by the linearised formulation documented in [22]. For a prescribed rotation up to π/9, as employed in this contribution, both formulations show quantitatively coinciding results. Also the semi-analytical solution, which has been obtained in [25] from the application of a power series method, is in quantitative agreement with the results as shown in [22]. Whereas these observations hold for the prescribed rotation of π/9 of the outer tube radius, a deviation from the small strain results in [22,25] is observable when significantly larger deformations are considered. In Fig. 5a-d, simulation results corresponding to rotations of π/9, π/6 and 2π/9 around the tube's middle axis are presented for λ * = 0.03. With an increasingly large deformation, the difference between the linearised and the non-linear modelling approach becomes more pronounced in all stress and couple-stress contributions considered.

Tensile test of a notched plate
A more complex example is represented by the fibre-reinforced notched plate shown in Fig. 6. By employing an offset between the notches, a bending deformation mode is activated as a tensile test is performed. The influence of the fibre bending stiffness parameter c as well as of the initial fibre direction field will be elaborated in this example.

Isogeometric model
In Fig. 6, the mesh for the notched plate is shown and the dimensions of the model are presented. The width of the inner part of the plate is w = 12 mm and the total length is l = 120 mm. Between the centre points of the circular notches with radii r = 3 8 w = 4.5 mm, an offset of d = 3 2 w = 18 mm is employed. The discretisation of the notched plate is obtained by n el = 990 elements. For the approximation of the placement field and the test function, NURBS basis functions with a polynomial degree of p = 4 are used in accordance with the previous example. The control polygon thus consists of n cp = 1410 control points. Within the whole domain, a global continuity of C p−1 , accordingly C 3 , is ensured due to the characteristic properties of NURBS, cf. [27][28][29].
In analogy to a uni-axial tensile test, the plate is clamped at both ends, meaning that displacements in e 2 -direction of the left and right boundary nodes are prevented. In e 1 -direction, a uniform displacement is prescribed on both sides. The total elongation of the plate is set to 1/5 of its length so that u = 12 mm. Similar to the previous example, it is assumed that no surface tractions or surface couples are present, i.e. t 1 = 0 and t 2 = 0.

Numerical results
Within the isogeometric finite element analysis of the notched plate, two different values of the fibre bending stiffness parameter are employed and two initial fibre orientations are considered. Specifically speaking, a reinforcement with fibres that are aligned in e 1 -direction as well as diagonally aligned fibres are taken into account for fibre bending stiffness parameters c ∈ {2 × 10 4 N, 3 × 10 4 N}. In the case of diagonally aligned fibres, the fibre direction vector exhibits an angle of α = π/4 to the e 1 -direction.
In Fig. 7, the deformed configuration of the notched plate is shown after the analysis under the above-mentioned boundary conditions without the presence of fibres. Due to the particular boundary conditions discussed in Sect. 4.3.1, the plate is stretched in e 1 -direction and the notches change their shape significantly. The middle part of the plate undergoes a bending deformation in consequence of the offset between the two notches. Figure 8 provides a more detailed view on the deformed mesh and compares the deformation pattern for c = 0 and c = 3 × 10 4 N for fibres aligned with the e 1 -direction, i.e. α = 0. The largest influence of the fibres on the overall deformation is obtained in the region of the notches. Due to the bending deformation, the fibre curvature part in the energy contribution (41) is activated. For the different fibre orientations, Fig. 9 presents the deformed mesh in more detail. It is shown that for fibres which are aligned with the e 1 -direction in the initial configuration, the deformation of the vertical element edges is rather similar to the case of a non-reinforced material. However, if a fibre angle of α = π/4 is employed in the initial state, a different deformation pattern is obtained in the region of the circular notches. In particular, the vertical element edges are bent into one preferred direction near the lower boundary, cf. Fig. 9b. 2 .5 π/9 π/6 2 π/9 radius r    example the skew-symmetric stresses take values of the same order as their symmetric counterparts so that both contributions have a significant impact on the total stress values. This relation, however, depends on the prescribed fibre bending stiffness parameter. For rather small values, the symmetric part of the stress tensor becomes more prominent in comparison to the skew-symmetric part, whereas for an increasing fibre bending stiffness, the skewsymmetric stresses exceed the symmetric contributions. For a fibre bending stiffness parameter of c = 3 × 10 4 N, as employed in the results in Fig. 14, the skew-symmetric part of the stress contribution σ 12 is dominant. Figure 15 presents the integrated reaction force in e 1 -direction over the prescribed displacement value for the different cases of initial fibre alignment and fibre bending stiffness. It can be observed that the curvature of the graphs, which show an overall non-linear behaviour, is slightly different depending on the initial fibre direction and fibre bending stiffness. Especially in the last simulation steps an alignment in e 1 -direction yields the highest reaction force. For higher values of the fibre bending stiffness parameter, the material response becomes stiffer and increasingly distinguishable from the case of a non-reinforced material.
Remark 1 In Sect. 4.1, the energy contribution belonging to a transversely isotropic material behaviour has been assumed to be constant in order to examine the influence of the higher-gradient contributions, respectively of the fibre bending stiffness, exclusively. If the influence of this energy part had additionally been taken into account instead, a more significant difference in the deformations and in the reaction force would have been expected when comparing the non-reinforced with the fibre-reinforced material.

Concluding remarks
In extension to the isogeometric finite element approach to fibre-reinforced composites possessing fibre bending stiffness developed in [22] subject to a small strain setting, an isogeometric approach to the general finite strain version of the problem has been presented in this contribution.
On the basis of the framework introduced in [9], a referential form of the fibre direction gradient has been included in the stored energy density function through an additional invariant. In particular, this gradient includes the fibre curvature as well as the gradient of the fibre stretch. For the incorporation of these higher-gradient contributions into the continuum model, the couple-stress theory has been employed and a partial differential equation of fourth order has been obtained.
Due to the continuity properties of NURBS basis functions, an isogeometric approach has been used to solve a discretised version of this equation within a finite element framework. To this end, a global continuity of at least C 1 has been realised for both examples presented in this contribution.
The analysis of a cylindrical tube under pure azimuthal shear with radially aligned fibres has shown that the results obtained by the presented method are in accordance with the finite element results from [20,22] as well as with the semi-analytical solution from [25] up to a certain amount of prescribed shear. For significantly large deformations, a deviation from the results obtained by the linearised formulation can, however, be observed. By the analysis of a plate with circular offset notches, a more detailed view on the impact of the fibre bending stiffness and initial fibre direction has been achieved. Especially in the region of the inhomogeneous deformation, the higher-gradient contributions to the stored energy density function significantly influence the deformation pattern as well as the stress and couple-stress distribution in the material. Depending on the value of the fibre bending stiffness parameter, the skew-symmetric and symmetric stress contributions can take similar values. Thus, both parts contribute significantly to the total stresses.
Overall, this contribution established a basis for the modelling of complex boundary value problems including fibre-reinforced composites with fibre bending stiffness in a general finite deformation setting. In future works, this basic approach shall be used for advanced constitutive modelling such as for materials with initially curved fibres and extended by employing additional deformation modes like fibre twist. Also complex three-dimensional structures shall be analysed with the proposed method.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. + J −1 I ⊗ I : [F · Λ · a 0 ⊗ a 0 ] · I ⊗ I : The sensitivities of the kinematic quantities with respect to the deformation gradient as well as to the second gradient of the placement function are ∂Λ ∂ F = I ⊗ a 0 · I ⊗ I : Υ : I ⊗ I + I ⊗ [F · ∇ X a 0 ] :