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

In the modelling of fibre-reinforced composites, it is well established to consider the fibre direction in the stored energy in order to account for the transverse isotropy of the overall material, induced by a single family of fibres. However, this approach does not include any length scale and therefore lacks in the prediction of size effects that may occur from the fibre diameter or spacing. By making use of a generalised continuum model including non-symmetric stresses and couple-stresses, the gradient of the fibre direction vector can be taken into account as an additional parameter of the stored energy density function. As a consequence, the enhanced model considers the bending stiffness of the fibres and includes information on the material length scale. Along with additional material parameters, increased continuity requirements on the basis functions follow in the finite element analysis. The isogeometric finite element method provides a framework which can fulfil these requirements of the corresponding weak formulation. In the present contribution, the method is applied to two representative numerical examples. At first, the bending deformation of a cantilever beam is studied in order to analyse the influence of the fibre properties. An increasingly stiff response is observed as the fibre bending stiffness increases and as the fibre orientation aligns with the beam’s axis. Secondly, a fibre-reinforced cylindrical tube under a pure azimuthal shear deformation is considered. The corresponding simulation results are compared against a semi-analytical solution. It is shown that the isogeometric analysis yields highly accurate results for the boundary value problem under consideration.


Introduction
In a wide range of applications, materials are reinforced with fibres in order to improve their mechanical behaviour, cf. [14,17,18,27]. Electrical and thermal properties of a material can also be manipulated by the reinforcement with certain types of fibres, see [1,34].
When fibre-reinforced composites are modelled numerically, the fibres are, in most cases, characterised by a direction vector only. This leads to a classic structural tensor approach under the assumption of perfectly flexible fibres, cf. [5,19,20,29,30,36]. Although this concept yields accurate results for many applications, it cannot account for size effects that follow from particular fibre properties such as their diameter or spacing. On these grounds, a new approach has been presented in [31]. It drops the assumption of perfectly flexible fibres and instead allows the incorporation of the fibre bending stiffness into the material model. Amongst others, this new material parameter is associated with the gradient of the fibre direction which introduces a length scale to the model. As a special case, only the directional derivative in the fibre direction is considered to account for fibre curvature effects.
The incorporation of the gradient of the fibre direction in the stored energy density function causes the requirement of a generalised continuum theory. Several approaches to generalised continua have been examined in the past. In [11,21] strain gradient elasticity has been studied. Micromorphic and micropolar concepts have been addressed in [10,13,15]. In the present contribution, the couple-stress theory is employed. The fundamental balance equations in this context have been derived in [22]. A partial differential equation of fourth order with, in general, non-symmetric stresses and couple-stresses follows accordingly.
For the solution of the partial differential equation within a finite element framework, enhanced continuity requirements regarding the underlying basis functions need to be considered, cf. [12,35]. A multi-field approach using Lagrangian basis functions has been elaborated in [2][3][4]. In [23,24], a formulation based on special types of elements, such as Hermite elements, is proposed to account for the higher continuity demand. In the present contribution, the isogeometric analysis, presented in [6,16], is employed. Therein, non-uniform rational B-Splines (NURBS) are used as basis functions. The fundamental theory of NURBS is provided in [26], and it has been shown that these functions can fulfil the particular continuity requirements. The isogeometric approach has been applied to gradient elasticity problems in [11,32], and fibre-reinforced solids have been modelled within the isogeometric concept in [28]. In the present contribution, the model of fibre-reinforced composites with fibre bending stiffness introduced in [31] will be analysed within an isogeometric finite element framework. A major advantage of this approach is the possibility to discretise and solve the fourth-order partial differential equation directly as opposed to the multi-field method where the problem is split into two partial differential equations of second order. This leads to less degrees of freedom in the isogeometric approach which lowers the computational effort.
We will start with the presentation of the kinematics as well as the balance equations for the generalised continuum approach in Sect. 2 and recapitulate the constitutive relations for fibre-reinforced composites including fibre bending stiffness. In Sect. 3, a weak form of the governing equation is developed and the discretisation is performed within a finite element framework. Section 4 provides details on the isogeometric analysis including higher-order derivatives of the NURBS basis functions and the corresponding coordinate mappings. The main focus lies on the fulfilment of the continuity requirements that follow from the fourth-order partial differential equation and the weak form associated with it. The presented finite element method is applied to representative numerical examples in Sect. 5. A beam subject to a bending deformation is studied, and a cylindrical tube is analysed subject to a pure azimuthal shear deformation. The results from the isogeometric analysis are evaluated with special regard to the influence of the fibre properties, and, particularly for the cylindrical tube, they are compared against a semi-analytical solution that is provided in [7].
Notation Let T and S denote tensor-valued quantities of arbitrary order in the Cartesian basis system. The single tensor contraction is introduced as if T and S are of at least third order. For first-order tensors a and b, the vector product is denoted under the use of the Levi-Civita tensor so that For a second-order tensor T , the product reads The standard dyadic product is defined as and two non-standard dyadic products are introduced in the form of and The application of the gradient operator to a tensor T of arbitrary order yields and applying the divergence operator results in where I is the second-order identity tensor. The curl of a first-order tensor a is denoted as Different kinds of transpositions of tensors are defined depending on the order of the tensor. For a second-order tensor T , the transposition results in For a fourth-order tensor S, two different types of transpositions are used as this work proceeds, namely Similarly, for a third-order tensor R, are employed.

Fibre-reinforced composites with fibre bending stiffness
This section presents the modelling approach for fibre-reinforced composites. It is assumed that the fibres are convected with the matrix material and that they possess a bending stiffness. The governing equations have been derived in [31] for the general case of finite strains as well as for a linearised setting, which will be presented here.

Small strain kinematics
For the purpose of including the fibre bending stiffness in the model, a generalised continuum approach is considered. In the case of a small strain theory, it takes into account not only the first but also the second gradient of the unknown displacement field u. In particular, the analysis requires the small strain tensor as the symmetric part of the first gradient of the displacement field. For the linearised theory, all terms of order O(e n ) with e 1 and n ≥ 2 are neglected. Moreover, all derivatives of the displacement field are supposed to be of order O(e). Considering these properties, the vector is introduced. It incorporates the fibre direction which is described by the vector a with ||a|| = √ a · a = 1. The gradient of the vector β can be obtained as (17) under the assumption that the fibres are initially straight and homogeneously distributed. It becomes evident that this second-order tensor γ describes the projection of the second gradient of the displacement field onto the direction of the fibres and, as a consequence, includes the fibre curvature.

Balance equations
The consideration of second-order gradients of the displacement field in the modelling approach motivates the employment of a generalised continuum theory, such as the couple-stress theory. The corresponding balance equations have been derived in [22] and will be summarised in the following sections.

Balance of mass
For a body B ⊂ R 3 , as sketched in Fig. 1, the total mass is preserved assuming that B represents a closed system. Consequently, the balance of mass in integral form reads with the mass density denoted as ρ. The application of Reynolds transport theorem leads to the local form which can be denoted asρ with the velocity vector v =u and with• representing the time derivative of •.

Balance of linear momentum
The balance of linear momentum takes the body forces b that act within the body B, as well as the tractions t 1 on its surface ∂B into account. The sum of these force contributions equals the rate of change of linear momentum, such that the integral form of the balance equation follows as The boundary tractions can be rewritten in terms of Cauchy's theorem, accordingly where n represents the outward surface normal vector and σ is the stress tensor. Applying the divergence theorem to the integral form (20), the local form of the balance of linear momentum is obtained as

Balance of angular momentum
In analogy to the balance of linear momentum, the rate of change of angular momentum equals the sum of all couples that act on B. This can be expressed by the integral form of the balance of angular momentum, such that where r = x − x r is the distance between the spatial position x and a reference position x r . Vector c represents the body couples. Similar to the surface tractions t 1 , the couple-stress that act on the surface of the domain B are introduced as with the couple-stress tensor m. The application of the divergence theorem yields the local form of the balance equation, i.e.
A general outcome of the balance of angular momentum for generalised continua is that the skew-symmetric stresses are nonzero and can be derived using (25). It follows that the skew-symmetric part of the stress tensor can be expressed in terms of the divergence of the couple-stress tensor, accordingly Inserting this relation into the local form of the balance of linear momentum (22) leads to a partial differential equation of fourth order which can be denoted as This expression represents the governing equation for the couple-stress theory in the format that is also employed in [22].

Balance of energy
Finally, the balance of mechanical energy is considered. Contributions associated with the thermal energy are neglected. With the internal energy U , the corresponding balance equation reads The local form is obtained by using the balance equations of momentum in the form of (22) and (25), so that It becomes evident that the skew-symmetric part of the stress tensor does not occur in the balance equation.
Similarly, as further elaborated in [22], the spherical part of the couple-stress tensor does not contribute to the energy balance. As a consequence, this part remains indeterminate. As this work proceeds, the term couplestress tensor will refer to the deviatoric part of the tensor for the sake of simplicity, while the spherical part will be assumed to be zero.

Material model
In the theory of fibre-reinforced composites with fibre bending stiffness, the stored energy takes into account not only the strain and the fibre direction vector, as in a classic structural tensor approach in transverse isotropy. It additionally considers the tensor γ , as defined in (17) for the linearised setting, that incorporates the fibre curvature. A representation of the stored energy density function in terms of n I invariants I i , which are dependent on {ε, a, γ }, is possible, i.e.
For this specific format of the stored energy density function, it follows from the derivations in [31] that the symmetric part of the stress tensor takes the form Together with the skew-symmetric counterpart, as specified in (26), the total stress tensor follows as The couple-stress tensor can be derived from the stored energy density function in a similar manner. It can be expressed as Remark 1 Although the vector β does not occur explicitly in the parameter list of the stored energy density function (31), it is included in the general expression (34) for the couple-stress tensor m. The related derivation is presented in [31, (9.9)], in analogy to the finite strain version, and a linear form is proposed [31, (9.12)]. The latter will be used as this work proceeds, cf. Sect. 3.2.

Finite element formulation
Within the finite element analysis, a discretised weak form of the fourth-order partial differential equation (27) is developed. In the following sections, the discretised kinematic quantities as well as the tangent operators are derived for the particular constitutive model and the resulting discretised weak form is presented.

Weak form
For the finite element formulation of the governing equation (27), a quasi-static case is assumed. Moreover, the body forces as well as the body couples and the couples acting on the boundary are neglected, i.e. b = c= t 2 = 0.
In order to obtain a weak form, the resulting equation is multiplied by a test function η and integrated over the domain B. After applying the divergence theorem twice, the weak form results in The notation ∇ 2 x • = ∇ x ∇ x • is employed for the second gradient.

Discretisation of the weak form
Within the concept of the finite element method, the weak form (35) is discretised by using basis functions R • . The discretisation of the test function and of the kinematic quantities for an element e yields wherein n en is the number of basis functions that have support on the particular element. Inserting the first relation into (35) and taking into account an assembly over the elements related to the domain B by the assembly operator A , the discretised weak form of the governing equation results in with the internal and external force vectors and Therein, the total number of elements is denoted as n el .
Based on the material model described in Sect. 2.3, the tangent operators for the finite element analysis can be derived. Recapitulating the assumption on higher-order terms made in Sect. 2.1, a linear relation between the tensors σ sym and ε as well as m and κ can be obtained. This results from the derivations in [31, (9.9), (9.12)] where, in particular, the higher-order terms that follow from the dyadic products in (34) are neglected. The tangent operators are accordingly With these formulations at hand, the symmetric part of the stress tensor as well as the couple-stress tensor can be expressed as which reflects the assumption of linear constitutive relations. Considering that, for the small strain setting, the internal force vector is a linear function of the displacement field, these relations can be used directly in the discretised weak form (39) together with the expressions (37) and (38), specifically It becomes evident that in the discretised equations, second-order derivatives of the basis functions are present. This leads to enhanced continuity requirements in the sense that global C 0 -continuity, as provided by classic Lagrangian basis functions, is not sufficient. These requirements are met by the isogeometric analysis. More details on this topic are provided in Sect. 4.

Isogeometric analysis
Due to the particular continuity requirements that follow from the incorporation of the fibre bending stiffness into the continuum model, the isogeometric method is employed for the finite element analysis. The fundamentals of this approach are described in detail in [6]. In [26], the topic of non-uniform rational B-Splines (NURBS), which serve as the basis functions in the isogeometric analysis, is addressed. In the following sections, the fundamentals of the isogeometric method are summarised with a special focus on higher-order derivatives of NURBS as well as their continuity properties.

Knot vectors and NURBS basis functions
B-Splines are defined on a parametric space which is spanned by knot vectors such as This format refers to a so-called open knot vector, since the first and last knots appear p + 1 times, where p represents the polynomial degree of the corresponding basis functions. In the two-dimensional case, one knot vector is required for each of the two parametric directions ξ and η. B-Splines are calculated from the knot vectors in a recursive manner depending on their polynomial degree. For p = 1, the expression for the one-dimensional B-Spline functions reads For all higher degrees, B-Splines are constructed recursively following Although B-Splines could be used as basis functions in a finite element analysis, the use of their rational form provides some advantages. In particular, circular structures, like in the example which will be addressed in Sect. 5.3, can be precisely modelled with NURBS as opposed to B-Splines. The difference between the two formulations is the consideration of weights w i, j , so that the two-dimensional NURBS basis functions are expressed as The numbers n and m denote the number of B-Spline basis functions in the respective parametric directions ξ and η with the polynomial degrees p and q.

Derivatives of NURBS
For the analysis of fibre-reinforced composites with fibre bending stiffness, one main focus lies on the higherorder derivatives of the basis functions. In [26], the derivatives of NURBS with respect to the parametric coordinates ξ and η are derived as This general formulation yields the Kth (Lth) derivative with respect to ξ (η) and includes the derivatives and The derivatives of the one-dimensional B-Spline functions follow from with the coefficients specified from the knot vector. The derivatives of the B-Spline basis functions corresponding to the second parametric direction η are obtained analogously.

Coordinate mappings
In the discretised form of the governing equation (39), the integrals are solved numerically by means of Gaussian quadrature. Consequently, a mapping between the master element, on which the integration is performed, and the physical domain is necessary. In the isogeometric analysis, actually two different mappings are required for this purpose, cf. Fig. 2. In [6], both mappings are discussed in detail. At first, the master element is connected to the parametric space, on which the NURBS basis functions are defined. For an element e = [ξ i , ξ i+1 ] ⊗ [η j , η j+1 ], the corresponding Jacobian matrix is obtained by Secondly, a relation between the parametric and the physical space is established. To this end, a so-called control polygon, as sketched in Fig. 7, is considered. Matrix B contains the respective control points. As there are n cp = n·m control points, this matrix takes the dimension [n cp × n dm ] with n dm = 2 for the twodimensional case. Analogously, the NURBS basis functions will be stored in a vector R of dimension n cp . With these quantities at hand, the mapping between the vectors x = [x y] t and ξ = [ξ η] t can be written as The final form of the Jacobian determinant includes both contributions, accordingly With the relations (54) and (55), the derivatives of the NURBS basis functions, which have been determined in (49), can be mapped to the physical domain. Applying the chain rule, the first-and second-order derivatives with respect to the physical coordinates follow as and The expression for the third-order derivatives is obtained analogously and is denoted in "Appendix A". The third-order derivatives are not directly necessary for the analysis but only employed for the calculation of the skew-symmetric stresses for postprocessing purposes. In "Appendix A", all derivatives that are used in (57) and (58) are specified as well.

Continuity of NURBS basis functions
From the discretised equation (39), it has been shown that second-order derivatives of the basis functions are required within the finite element framework. Consequently, C 0 -continuity is not sufficient from a mathematical point of view. Using NURBS as basis functions ensures C ∞ -continuity within elements and across element boundaries. This motivates the use of the isogeometric finite element method for boundary value problems including higher gradient contributions. However, the continuity is decreased at the location of double knots as well as repeated control points. By additional procedures which will be presented in the following sections, C 1 -continuity is obtained at these locations. Global C 1 -continuity is sufficient for the approximation of the discretised form of the governing equation.

Mesh modification strategies
In isogeometric finite element approaches, the mesh and the corresponding basis functions are determined by the knot vectors. In particular, the number of nonzero knot spans specifies the number of elements in the mesh, whereas the multiplicity of reoccurring knots affects the continuity of the basis functions at these specific locations. Consequently, the knot vectors need to be adjusted in order to obtain the desired mesh characteristics. The corresponding procedures are derived in [26] and will be summarised in the following sections.

Knot insertion
Focussing on mesh refinement strategies at first, a similar procedure to those in classic finite element approaches can be found in the so-called knot insertion. Subsequently, new knots ξ ∈ [ξ k , ξ k+1 ) are inserted into the knot vector leading to an increasing number of elements and possibly to more accurate analysis results. When a knot vector is expanded due to knot insertion, the control polygon has to be adjusted. Let be the ith control point of the jth curve as defined in "Appendix B". For the calculation of the new control points and weights of this NURBS curve, the two-dimensional weighted control points are introduced. On the basis of this definition, the new set of control points and weights is calculated according to where For a two-dimensional NURBS surface, the algorithm is applied to each of the m NURBS curves. The same procedure is followed when inserting a knot into the knot vector for the second parametric direction η.

Knot removal
The issue of a decreased basis function continuity at the locations of repeated knots is addressed by the procedure of knot removal. The general property of NURBS states that they show C ∞ -continuity everywhere, including inter-element locations. However, an exception is observed at points where a knot occurs multiple times in the knot vector. At these points, the continuity is decreased to C p−k with k being the multiplicity of the knot. The aim of knot removal in this contribution is the increase in continuity at repeated knot locations in a way that at least C 1 -continuity is ensured. In this way, the continuity requirements for the analysis of fibre-reinforced composites with fibre bending stiffness, discussed in Sect. 4.4, are fulfilled.
For the process of removing knots from the knot vector of a NURBS curve, the procedure of knot insertion from the previous section is reversed. If a knot ξ = ξ r = ξ r +1 of initial multiplicity s is removed t times, the new control points in the tth removal step can be found as with the constant coefficients Due to the local support of NURBS, the remaining control points are not affected by the removal process and do, as a consequence, not change.
Remark 2 It is, in general, not known a priori if a knot is removable from the knot vector without changing the geometry. In every step of the knot removal process, this has to be examined individually. Detailed information on this method is provided in [26,Sect. 5.4].

Linear constraints
In order to obtain a C 1 -connection at repeated control points, the approach described in [11] is employed. Control points are especially repeated when a closed geometry such as a circle is modelled using NURBS. In this case, the first and the last control points of each curve are equal, i.e.
This yields a C 0 -continuous connection between the starting and end point of the NURBS curve. Equally as in the knot refinement and removal procedures, the condition generally applies to the weighted control points B w i of each curve specified in (60). In order to achieve C 1 -continuity, a second condition must be fulfilled, namely Therein, the constant factor k is determined by the corresponding knot vector. For the example presented in Sect. 5.3, the knot vector as well as the control points follow a symmetric pattern so that k = 1 is found. This is similar to the derivation in [6] for the special case of Bézier curves. Moreover, the weights that correspond to the control points employed in (66) are all equal to 1 so that, instead of the weighted control points, the condition can be applied directly to the respective points B i of the control polygon.
Since C 1 -continuity implies that the connection is also C 0 -continuous, (65) can be inserted into (66). Consequently, the control point B 1 can be expressed in terms of the neighbouring points according to in order to ensure C 1 -continuity at the location of this repeated control point. Within the isogeometric finite element framework, the relation (67) is implemented as a linear constraint following the procedure described in [33]. In particular, the constraints are applied to the displacement values at the respective control points and are also accounted for in their initial positions.

Inhomogeneous Dirichlet boundary conditions
Unlike Lagrange polynomials, NURBS basis functions are not necessarily interpolatory at the entire boundary of the domain. Accordingly, a specific technique is required to accurately impose inhomogeneous boundary conditions. In [9], two possibilities to prescribe inhomogeneous Dirichlet boundary conditions have been introduced. The procedures are based on the relation between the control points and respective points of the physical mesh. For one mesh point x, the relation can be denoted as Considering the n control points B on the outer boundary of the domain only and choosing a number of n sp interpolating points x on the outer mesh surface, the expression is expanded to with the matrices Since the isogeometric analysis comprises the isoparametric concept, the same structure holds for the displacements u that are associated with either the control points or the mesh points, i.e.
For the direct solution of the equation system in (71), matrix R must be quadratic and invertible. To this end, the number of interpolating points on the mesh must be chosen to be equal to the number of control points, and an adequate distribution of the mesh points has to be ensured, cf. [9]. The calculated displacement values u cp are then imposed to the respective control points as inhomogeneous Dirichlet boundary conditions and will lead to the desired displacements u sp of the actual boundary points of the body.

Representative numerical examples
In the following section, two numerical examples are presented. In both cases, the model for fibre-reinforced composites that has been elaborated in this contribution is employed. At first, the bending deformation of a fibre-reinforced beam is analysed in order to examine the influence of the fibre bending stiffness and of the fibre orientation on the mechanical behaviour of a representative body under load. In a second step, focus is on the azimuthal shear deformation of a cylindrical tube for which a semi-analytical solution can be derived. A comparison between the latter and the numerical results obtained from the isogeometric finite element method is made in order to validate the proposed isogeometric finite element formulation.

Specific form of the energy function
As proposed in [3], an additive decomposition of the stored energy density function according to is considered. For the isotropic part, resembling the matrix material, a quadratic form in the small strain tensor is assumed, specifically speaking with the invariants The second part W λ describes the transversely isotropic behaviour of the material reinforced by a single family of fibres. However, this part is assumed to be constant because the focus lies on the examination of the fibre curvature and fibre bending stiffness. These properties are considered in the last part of the above given stored energy density function. Specifically, this part is determined as with the fibre bending stiffness parameter c, and with the invariant that includes the fibre curvature as the projection of the second gradient of the fibre vector onto the fibre direction.
From this particular form of the stored energy density function, the symmetric part of the stress tensor as well as the couple-stress tensor can be obtained by using the linear constitutive relations (43) together with the specific tangents and The resulting expressions are and In the present examples, the Lamé constants are specified to λ = 1.037 × 10 5 N mm −2 and μ = 4.444 × 10 4 N mm −2 .
In comparison with the model proposed in [7] with reference to the particular example of a cylindrical tube with radially aligned fibres, the constant factor in (80) correlates with the therein defined material parameter d f , i.e.
In further accordance with [7], the non-dimensional quantity is defined. It becomes evident that a length scale is introduced since the inner and outer radii R i and R a of the cylinder occur in the definition of this material parameter which is connected with the curvature of the fibres. In Sects. 5.2.3 and 5.3.3, the influence of the respective fibre bending stiffness parameter c or λ * on the stress and couple-stress components is examined for the numerical examples.

Bending of a fibre-reinforced beam
A two-dimensional beam, as shown in Fig. 4, is bended under the action of a force at its free end. The beam is reinforced by a single family of fibres under the assumption that the fibres are embedded in the matrix material and that they resist bending. Overall, a linearised setting and a plane strain state are assumed. The purpose of this example is the investigation of the fibre properties. Specifically speaking, the influence of the fibre bending stiffness and of the orientation of the fibre direction vector is examined. Since only two distinct knot values are considered, the coarsest mesh for the beam simply consists of one finite element. The corresponding control points and weights are provided in "Appendix B", and Fig. 3 represents the respective control polygon. For the simulation, a refined mesh with additional control points is used. It is shown in Fig. 4. The refinement is performed by means of knot insertion as described in detail in Sect. 4.5.1.

Boundary conditions
As shown in Fig. 4, the beam under consideration is fixed on its left end. To be more precise, the displacement in both coordinate directions as well as the couple-stress on this surface are assumed to be zero. On the opposite side, the beam is loaded by a traction t = 30 N mm −2 which is uniformly distributed over the right surface and accounted for in the analysis by means of Neumann boundary conditions.

Numerical results
The aim of applying the isogeometric analysis to the fibre-reinforced beam including fibre bending stiffness is to elaborate the influence of the particular properties of the fibres. In a first step, the fibres are assumed to be aligned with the beam's axis, i.e. the fibre angle is α = 0. For this setting, the fibre bending stiffness is varied. In a second step, a constant value of the fibre bending stiffness parameter is chosen, while different values for the fibre angle are examined.
Influence of the fibre bending stiffness The fibre bending stiffness is represented by the parameter c in the part of the energy function (75) that is related with the curvature of the fibres. For the limiting case c = 0, the underlying model coincides with the case of perfectly flexible fibres. In Fig. 5, the vertical displacement u P y of the tip of the beam, represented by the point P in Fig. 4, is depicted with a logarithmic scale for the parameter c. It is observed that the displacement takes its maximum value for a vanishing fibre bending stiffness. By a variation of c in the range of c ∈ 0, 10 7 N , an S-shaped curve is obtained. With an increasing bending stiffness of the fibres, the bending of the beam decreases monotonously. Moreover, it is observed that the slope of the curve decreases for high values of c and a that second limiting value is approached. Influence of the fibre orientation Complementing the previous example, different fibre orientations are analysed in the following, while the fibre bending stiffness parameter is c = 10 5 N. In particular, the fibre angle α = arccos(a · e x ), i.e. the angle between the fibre vector a and the horizontal direction, takes values between α = 0 and α = π/2. A monotonous increase in the vertical tip displacement can be observed when the fibre angle is increased, cf. Fig. 6. Similar to the results with varying fibre bending stiffness, an S-shaped curve is obtained. The minimum resistance against bending is found for α = π/2 so that the fibres lie vertically in the beam, i.e. in the direction of the applied force. On the other hand, when the fibres are aligned with the beam's axis, the impact of the fibre bending stiffness is most noticeable. In particular, a reduction of the vertical displacement at point P of about 18 % as compared to a beam that is reinforced with fibres in vertical direction is observed.

Fibre-reinforced cylindrical tube under pure azimuthal shear
As a more advanced example, a cylindrical tube subject to a pure azimuthal shear deformation is analysed under the assumptions of plane strain and a linearised setting. In [7], the boundary value problem has been presented in detail and a semi-analytical as well as a numerical solution has been provided for different properties of the fibre-reinforced material. We will focus here on composites in which the fibres are convected with the matrix material and aligned in the radial direction. For this case, it has been found in [7] that azimuthal shear strain and radial stretching become uncoupled. This observation holds regardless of the assumptions made on the extensibility of the fibres and the compressibility of the matrix. In this contribution, the problem of pure azimuthal shear on a cylindrical tube is addressed in terms of the isogeometric analysis and by applying the material model specified in Sect. 5.1. In the following, the geometry model for the cylinder is described and the boundary conditions for the finite element analysis are outlined. The numerical results are compared to the semi-analytical solution provided in [7] in order to evaluate the accuracy of the isogeometric finite element approach. Additionally, a convergence study has been performed and is provided in "Appendix C".

Boundary conditions
For the solution of the fourth-order partial differential equation (27), two different sets of boundary conditions are presented in [7]. They are associated with a pure azimuthal shear deformation of the cylindrical tube. In both cases, the inner radius of the cylinder is kept fixed, while the points on the outer radius are subjected to a circumferential displacement resulting in a shear deformation. Besides these Dirichlet conditions, the couple-stress is forced to vanish on the outer boundary. On the inner boundary, either the same assumption is employed or the fibre slope is set to zero. For the results that will be presented in the following section, the first option is considered, so that the complete set of boundary conditions can be summarised as using cylindrical coordinates. The latter two conditions are in accordance with the assumption of a vanishing surface couple-stress from Sect. 3.1.

Numerical results
The isogeometric analysis of the fibre-reinforced cylindrical tube under pure azimuthal shear has been performed on the geometry shown in Fig. 8. The tube with outer radius R a = 2.5 mm and inner radius R i = 1.0 mm is discretised by n el = 6912 finite elements. The prescribed azimuthal displacement is u ϕ = 0.025 mm. The simulation is performed for different values of the non-dimensional parameter λ * introduced in (82), namely λ * ∈ {0.0, 0.005, 0.03, 0.1, π}. A semi-analytical solution to the boundary value problem has been obtained by means of a power series method in [7].
Dimensionless quantities In order to compare the results of the isogeometric analysis with the semi-analytical results published in [7], dimensionless quantities are employed. The conversion factors that will be used in the following have also been specified in [4]. Using the ratio the dimensionless stress contributions follow as  Fibre deformation and slope By applying an azimuthal displacement on the outer radius of the cylinder while keeping the inner boundary fixed, a deformation of the fibres, which are convected with the matrix material, is observable. Figure 9 shows the deformation pattern of a radial fibre that was initially aligned with the horizontal direction. It becomes evident that an increase in the parameter λ * , which is associated with the fibre bending stiffness, leads to less bending of the fibres. In fact, for large values such as λ * = π, the fibres almost remain straight.
In Table 1, the fibre slope at the inner boundary, i.e. r = R i = 1.0 mm, is examined. The values corresponding to the semi-analytical solution are taken from [7] for comparison. It is found that the deviation is smaller than 0.4% for all considered values of λ * .
Stress and couple-stress components For the two-dimensional boundary value problem of pure azimuthal shear, the relevant stress and couple-stress components are σ r ϕ , σ ϕr and m r z . In addition to the components of the total stress tensor, the values for its symmetric part, namely σ sym r ϕ = σ sym ϕr , are analysed in this section and shown in Fig. 10 together with the couple-stress. Comparing the results of the isogeometric analysis against the semi-analytical values, taken from [7] and presented in Fig. 11, no significant difference can be observed. Accordingly, a high accuracy of the isogeometric approach can be deduced for the particular boundary value problem studied.
A more detailed view on the couple-stress component m r z is given in Fig. 12. For λ * = 0, the couple-stress remains zero, because the higher-order terms are not activated. For all other values of λ * , the couple-stress  [7] increases with increasing values of the non-dimensional material parameter. Due to the particular boundary conditions discussed in Sect. 5.3.2, the couple-stress m r z vanishes at the boundaries of the tube.
The total stress tensor includes not only the symmetric but also a skew-symmetric part. Figure 13 shows how the resulting stress components σ r ϕ and σ ϕr are distributed over the radius of the cylinder. In Fig. 14, the results from [7] are employed. When regarding the results from the numerical as well as the semi-analytical calculation, a similar profile along the radius can be observed. Both figures show that the stress component σ r ϕ is quite similar for all fibre bending stiffness values. For σ ϕr on the contrary, a monotonous decrease is obtained for perfectly flexible fibres and very small values of λ * , whereas for higher values, a parabolic shape is observable near the inner boundary. In further comparison of Figs. 13 and 14, slightly different values occur especially at the inner radius of the cylinder. Since the deviation is smaller than 5%, this might be caused by numerical errors. It is further noted in [8] that one term is missing in the expression of the skew-symmetric stresses in [7]. However, the results presented therein are considered to be very accurate.

Concluding remarks
The present contribution focusses on the simulation of fibre-reinforced composites where the fibres exhibit a certain bending stiffness. The material model has been developed in accordance with the derivations presented in [31]. For a small strain setting, the stress and couple-stress tensor have been specified and the balance equations for a couple-stress theory have been derived. An isogeometric finite element framework has been developed for solving the resulting partial differential equation of fourth order because the particular continuity requirements can be fulfilled by the use of NURBS as basis functions.
With the isogeometric framework, two numerical examples have been examined. In the first example, a fibre-reinforced beam subject to a bending deformation has been simulated in order to study the influence of the fibre properties. It has been confirmed that an increasing fibre bending stiffness parameter leads to a stiffer response. In addition, the influence of the fibre orientation has been studied in detail. The minimum bending has been obtained for the fibres being aligned with the beam's axis.
In the second example, a cylindrical tube under pure azimuthal shear has been analysed. Due to the assumption of radially aligned fibres, the equations for the displacement in radial and azimuthal direction become uncoupled. Under appropriate boundary conditions, the model has been simulated by the isogeometric finite element approach. In order to evaluate the respective simulation results, they have been compared against the semi-analytical solution provided in [7]. It has been found that the isogeometric analysis yields highly accurate results in the modelling of fibre-reinforced composites including fibre bending stiffness. In particular for the symmetric stress components as well as for the couple-stress components, no significant difference to the semi-analytical results is observable when plotted over the radius of the cylindrical tube.
Moreover, it has been shown that the continuity properties of the basis functions have a significant impact on the simulation. By the knot removal procedure as well as the imposition of linear constraints for repeated control points, the continuity requirements for the analysis of materials including higher gradient contributions can be fulfilled globally within the isogeometric analysis.
In future research, the developed isogeometric finite element framework including higher gradients of the displacement field shall be extended to finite deformations by employing the respective constitutive relations specified in [31]. In addition, more complex boundary value problems, particularly in the three-dimensional space, shall also be studied.
Funding Open Access funding provided by Projekt DEAL.
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/.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.

Appendix A: Derivatives
In addition to the expressions derived in Sect. 4.3, the remaining derivatives up to third order are presented in the following. The third derivative of one NURBS basis function R with respect to the physical coordinates reads Therein, the relations between x and ξ are required. The mapping ∂ξ /∂ x can be obtained by the inversion of the Jacobian matrix J 2 specified in (55). The remaining relations are and as well as and The derivatives and are obtained in analogy to (55).

Appendix B: Control points and weights
As stated in Sect. 4.3, the control points used to construct a control polygon are stored in the list The geometry model for the beam contains n = m = 5 control points in each direction so that the total number of control points is n cp = n·m = 25. The respective coordinates are and all weights are equal to 1.
For the cylindrical tube, the total number of control points is n cp = n·m = 85 since the control polygon consists of m = 5 curves containing n = 17 control points each. The specific control points [B i ] j in the jth curve are given by and are the same for every curve.

Appendix C: Convergence studies
The results for the cylindrical tube, presented in Sect. 5.3.3, have been obtained with a discretisation based on n el = 6912 elements. The corresponding physical mesh is shown in Fig. 8. In order to investigate the convergence behaviour of the isogeometric finite element method with regard to the particular boundary value problem under consideration, two additional discretisations have been taken into account. In Fig. 15, the mesh consists of n el = 1856 elements with an element size that is approximately twice as large in both parametric directions compared to the previous mesh. In Fig. 16, an even coarser mesh with n el = 272 elements is considered.  Fig. 21 for all three meshes in one plot. The dimensionless material parameter which is related to the length scale is chosen to be λ * = 0.1 in this case. By comparing the different graphs, it can be observed that the solution converges for a sufficiently fine discretisation. From Fig. 21, it is concluded that the graphs for the two finer meshes almost coincide, while the results for the mesh with n el = 272 elements differ slightly in the boundary regions of the circular domain. For this discretisation, a non-smooth behaviour can be noted especially at the inner radius, which vanishes as the mesh size decreases. In the inner region of the cylindrical tube, the analysis yields an accurate solution already for a coarse mesh.