1 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 \({\varvec{T}}\) and \({\varvec{S}}\) denote tensor-valued quantities of arbitrary order in the Cartesian basis system. The single tensor contraction is introduced as

$$\begin{aligned} {\varvec{T}}\cdot {\varvec{S}} =&\, [T_{ij\ldots kl}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l] \cdot [S_{mn\ldots op}\,{\varvec{e}}_m\otimes {\varvec{e}}_n\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p] \nonumber \\ =&\, T_{ij\ldots kl}\,S_{ln\ldots op}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_n \ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p \end{aligned}$$
(1)

following Einstein’s summation convention. Analogously, the double and triple contractions are defined as

$$\begin{aligned} {\varvec{T}}:{\varvec{S}} =&\, [T_{ij\ldots kl}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l] :[S_{mn\ldots op}\,{\varvec{e}}_m\otimes {\varvec{e}}_n\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p] \nonumber \\ =&\, T_{ij\ldots kl}\,S_{kl\ldots op}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p \end{aligned}$$
(2)

if \({\varvec{T}}\) and \({\varvec{S}}\) are of second or higher order and

$$\begin{aligned} {\varvec{T}}:\!\cdot \,{\varvec{S}} =&\, [T_{ij\ldots klm}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l\otimes {\varvec{e}}_m] :\!\cdot \,[S_{nop\ldots qr}\,{\varvec{e}}_n\otimes {\varvec{e}}_o\otimes {\varvec{e}}_p\ldots \otimes {\varvec{e}}_q\otimes {\varvec{e}}_r] \nonumber \\ =&\, T_{ij\ldots klm}\,S_{klm\ldots qr}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_q\otimes {\varvec{e}}_r \end{aligned}$$
(3)

if \({\varvec{T}}\) and \({\varvec{S}}\) are of at least third order.

For first-order tensors \({\varvec{a}}\) and \({\varvec{b}}\), the vector product is denoted under the use of the Levi-Civita tensor \(\varvec{\epsilon }\) so that

$$\begin{aligned} {\varvec{a}}\times {\varvec{b}} = a_ib_j\epsilon _{ijk}\,{\varvec{e}}_k. \end{aligned}$$
(4)

For a second-order tensor \({\varvec{T}}\), the product reads

$$\begin{aligned} {\varvec{T}}\times {\varvec{b}} = T_{ij}b_k\epsilon _{jkl}\,{\varvec{e}}_i\otimes {\varvec{e}}_l. \end{aligned}$$
(5)

The standard dyadic product is defined as

$$\begin{aligned} {\varvec{T}}\otimes {\varvec{S}} =&\, [T_{ij\ldots kl}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l] \otimes [S_{mn\ldots op}\,{\varvec{e}}_m\otimes {\varvec{e}}_n\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p] \nonumber \\ =&\, T_{ij\ldots kl}\,S_{mn\ldots op}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l \otimes {\varvec{e}}_m\otimes {\varvec{e}}_n\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p \end{aligned}$$
(6)

and two non-standard dyadic products are introduced in the form of

$$\begin{aligned} {\varvec{T}}\,{\overline{\otimes }}\,{\varvec{S}} =&\, [T_{ij\ldots kl}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l] \,{\overline{\otimes }}\,[S_{mn\ldots op}\,{\varvec{e}}_m\otimes {\varvec{e}}_n\ldots {\varvec{e}}_o\otimes {\varvec{e}}_p] \nonumber \\ =&\, T_{ij\ldots km}\,S_{ln\ldots op}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l \otimes {\varvec{e}}_m\otimes {\varvec{e}}_n\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p \end{aligned}$$
(7)

and

$$\begin{aligned} {\varvec{T}}\,{\underline{\otimes }}\,{\varvec{S}} =&\, [T_{ij\ldots kl}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l] \,{\underline{\otimes }}\,[S_{mn\ldots op}\,{\varvec{e}}_m\otimes {\varvec{e}}_n\ldots {\varvec{e}}_o\otimes {\varvec{e}}_p] \nonumber \\ =&\, T_{ij\ldots kn}\,S_{lm\ldots op}\,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l \otimes {\varvec{e}}_m\otimes {\varvec{e}}_n\ldots \otimes {\varvec{e}}_o\otimes {\varvec{e}}_p. \end{aligned}$$
(8)

The application of the gradient operator to a tensor \({\varvec{T}}\) of arbitrary order yields

$$\begin{aligned} \nabla _{\!\varvec{x}}{\varvec{T}} = \frac{\partial T_{ij\ldots kl}}{\partial x_m} \,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k\otimes {\varvec{e}}_l\otimes {\varvec{e}}_m \end{aligned}$$
(9)

and applying the divergence operator results in

$$\begin{aligned} \nabla _{\!\varvec{x}}\cdot {\varvec{T}} = \nabla _{\!\varvec{x}}{\varvec{T}} : {\varvec{I}} = \frac{\partial T_{ij\ldots kl}}{\partial x_l} \,{\varvec{e}}_i\otimes {\varvec{e}}_j\ldots \otimes {\varvec{e}}_k, \end{aligned}$$
(10)

where \({\varvec{I}}\) is the second-order identity tensor. The curl of a first-order tensor \({\varvec{a}}\) is denoted as

$$\begin{aligned} \nabla _{\!\varvec{x}}\times {\varvec{a}} = \epsilon _{ijk} \frac{\partial a_{k}}{\partial x_j} \,{\varvec{e}}_i. \end{aligned}$$
(11)

Different kinds of transpositions of tensors are defined depending on the order of the tensor. For a second-order tensor \({\varvec{T}}\), the transposition results in

$$\begin{aligned} {\varvec{T}}^\mathrm {t} = T_{ji} \,{\varvec{e}}_i\otimes {\varvec{e}}_j. \end{aligned}$$
(12)

For a fourth-order tensor \({\varvec{S}}\), two different types of transpositions are used as this work proceeds, namely

$$\begin{aligned} {\varvec{S}}^\mathrm {t} = S_{jikl} \,{\varvec{e}}_i\otimes {\varvec{e}}_j\otimes {\varvec{e}}_k\otimes {\varvec{e}}_l\quad \hbox {and}\quad {\varvec{S}}^\mathrm {T} = S_{klij} \,{\varvec{e}}_i\otimes {\varvec{e}}_j\otimes {\varvec{e}}_k\otimes {\varvec{e}}_l. \end{aligned}$$
(13)

Similarly, for a third-order tensor \({\varvec{R}}\),

$$\begin{aligned} {\varvec{R}}^\mathrm {t} = R_{jik} \,{\varvec{e}}_i\otimes {\varvec{e}}_j\otimes {\varvec{e}}_k \quad \hbox {and}\quad {\varvec{R}}^\mathrm {T} = R_{kij} \,{\varvec{e}}_i\otimes {\varvec{e}}_j\otimes {\varvec{e}}_k \end{aligned}$$
(14)

are employed.

2 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.

2.1 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 \({\varvec{u}}\). In particular, the analysis requires the small strain tensor

$$\begin{aligned} \varvec{\varepsilon } = \frac{1}{2}\left[ {\varvec{h}} + {\varvec{h}}^\mathrm {t} \right] ~\text {with}~ {\varvec{h}} = \nabla _{\!\varvec{x}}{\varvec{u}} \end{aligned}$$
(15)

as the symmetric part of the first gradient of the displacement field. For the linearised theory, all terms of order \({\mathcal {O}}(e^n)\) with \(e\ll 1\) and \(n\ge 2\) are neglected. Moreover, all derivatives of the displacement field are supposed to be of order \({\mathcal {O}}(e)\). Considering these properties, the vector

$$\begin{aligned} \varvec{\beta } = {\varvec{a}} + \nabla _{\!\varvec{x}}{\varvec{u}}\cdot {\varvec{a}} \end{aligned}$$
(16)

is introduced. It incorporates the fibre direction which is described by the vector \({\varvec{a}}\) with \(\vert \vert {\varvec{a}}\vert \vert =\sqrt{{\varvec{a}}\cdot {\varvec{a}}}=1\). The gradient of the vector \(\varvec{\beta }\) can be obtained as

$$\begin{aligned} \varvec{\gamma } = \nabla _{\!\varvec{x}}\varvec{\beta } = \left[ \nabla _{\!\varvec{x}}\nabla _{\!\varvec{x}}{\varvec{u}}\right] \cdot {\varvec{a}} =: \varvec{\kappa } \cdot {\varvec{a}} \end{aligned}$$
(17)

under the assumption that the fibres are initially straight and homogeneously distributed. It becomes evident that this second-order tensor \(\varvec{\gamma }\) 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.

2.2 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.

2.2.1 Balance of mass

For a body \({\mathcal {B}} \subset {\mathbb {R}}^3\), as sketched in Fig. 1, the total mass is preserved assuming that \({\mathcal {B}}\) represents a closed system. Consequently, the balance of mass in integral form reads

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{{d}{{ t}}}} \int _{{\mathcal {B}}} \rho \,{\mathrm{d}}v = 0 \end{aligned}$$
(18)

with the mass density denoted as \(\rho \). The application of Reynolds transport theorem leads to the local form which can be denoted as

$$\begin{aligned} {\dot{\rho }}+\rho \,\nabla _{\!\varvec{x}}\cdot {\varvec{v}} = 0 \end{aligned}$$
(19)

with the velocity vector \({\varvec{v}} = \dot{{\varvec{u}}}\) and with \({\dot{\bullet }}\) representing the time derivative of \(\bullet \).

Fig. 1
figure 1

Domain of a body \({\mathcal {B}}\) in the couple-stress theory

Full size image

2.2.2 Balance of linear momentum

The balance of linear momentum takes the body forces \({\varvec{b}}\) that act within the body \({\mathcal {B}}\), as well as the tractions \({\varvec{t}}_1\) on its surface \(\partial {\mathcal {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

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}{{ t}}} \int _{{\mathcal {B}}} \rho \,{\varvec{v}} \,\mathrm{d}v = \int _{\partial {\mathcal {B}}} {\varvec{t}}_1 \,\mathrm{d}a + \int _{{\mathcal {B}}} \rho \,{\varvec{b}} \,\mathrm{{d}}v. \end{aligned}$$
(20)

The boundary tractions can be rewritten in terms of Cauchy’s theorem, accordingly

$$\begin{aligned} {\varvec{t}}_1 = \varvec{\sigma }^\mathrm {t} \cdot {\varvec{n}}, \end{aligned}$$
(21)

where \({\varvec{n}}\) represents the outward surface normal vector and \(\varvec{\sigma }\) 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

$$\begin{aligned} \rho \,\dot{{\varvec{v}}} = \nabla _{\!\varvec{x}}\cdot \varvec{\sigma }^\mathrm {t} + \rho \,{\varvec{b}}. \end{aligned}$$
(22)

2.2.3 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 \({\mathcal {B}}\). This can be expressed by the integral form of the balance of angular momentum, such that

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}{{ t}}}} \int _{{\mathcal {B}}} \rho \left[ {\varvec{r}} \times {\varvec{v}}\right] \mathrm{{d}}v&= \int _{\partial {\mathcal {B}}} \left[ {\varvec{r}} \times {\varvec{t}}_1 + {\varvec{t}}_2 \right] \mathrm{{d}}a + \int _{{\mathcal {B}}} \rho \left[ {\varvec{r}} \times {\varvec{b}} + {\varvec{c}} \right] \mathrm{{d}}v, \end{aligned}$$
(23)

where \({\varvec{r}}= {\varvec{x}}-{\varvec{x}}_{\mathrm{r}}\) is the distance between the spatial position \({\varvec{x}}\) and a reference position \({\varvec{x}}_{\mathrm{r}}\). Vector \({\varvec{c}}\) represents the body couples. Similar to the surface tractions \({\varvec{t}}_1\), the couple-stress that act on the surface of the domain \({\mathcal {B}}\) are introduced as

$$\begin{aligned} {\varvec{t}}_2 = {\varvec{m}}^\mathrm {t} \cdot {\varvec{n}} \end{aligned}$$
(24)

with the couple-stress tensor \({\varvec{m}}\). The application of the divergence theorem yields the local form of the balance equation, i.e.

$$\begin{aligned} {\varvec{0}} = \nabla _{\!\varvec{x}}\cdot {\varvec{m}}^\mathrm {t} + \rho \,{\varvec{c}} + \varvec{\epsilon }:\varvec{\sigma }. \end{aligned}$$
(25)

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

$$\begin{aligned} \varvec{\sigma }^{\mathrm{{skw}}}= \frac{1}{2} {\varvec{I}} \times \left[ \nabla _{\!\varvec{x}}\cdot {\varvec{m}}^\mathrm {t} + \rho \,{\varvec{c}} \right] . \end{aligned}$$
(26)

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

$$\begin{aligned} \rho \,\dot{{\varvec{v}}}&= \nabla _{\!\varvec{x}}\cdot \varvec{\sigma }^{\mathrm{{sym}}}+\rho \,{\varvec{b}} + \frac{1}{2} \nabla _{\!\varvec{x}}\times \left[ \nabla _{\!\varvec{x}}\cdot {\varvec{m}}^\mathrm {t} + \rho \,{\varvec{c}} \right] . \end{aligned}$$
(27)

This expression represents the governing equation for the couple-stress theory in the format that is also employed in [22].

2.2.4 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

$$\begin{aligned}&\frac{\mathrm{d}}{\mathrm{{d}{{ t}}}} \int _{{\mathcal {B}}} \rho \left[ \frac{1}{2} {\varvec{v}} \cdot {\varvec{v}} + U \right] \mathrm{{d}}v = \int _{\partial {\mathcal {B}}} \left[ {\varvec{t}}_1 \cdot {\varvec{v}} + \frac{1}{2} \, {\varvec{t}}_2 \cdot \nabla _{\!\varvec{x}}\times {\varvec{v}} \right] \mathrm{{d}}a + \int _{{\mathcal {B}}} \rho \left[ {\varvec{b}} \cdot {\varvec{v}} + \frac{1}{2}\,{\varvec{c}} \cdot \nabla _{\!\varvec{x}}\times {\varvec{v}}\right] \mathrm{{d}}v. \end{aligned}$$
(28)

The local form is obtained by using the balance equations of momentum in the form of (22) and (25), so that

$$\begin{aligned} \rho \, {\dot{U}} = \varvec{\sigma }^{\mathrm{{sym}}}: \nabla _{\!\varvec{x}}{\varvec{v}} + \frac{1}{2}{\varvec{m}}^\mathrm {t} : \nabla _{\!\varvec{x}}\nabla _{\!\varvec{x}}\times {\varvec{v}}. \end{aligned}$$
(29)

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 couple-stress tensor will refer to the deviatoric part

$$\begin{aligned} \overline{{\varvec{m}}} = {\varvec{m}}-\frac{1}{3}{\text {tr}}({\varvec{m}})\,{\varvec{I}} \end{aligned}$$
(30)

of the tensor for the sake of simplicity, while the spherical part will be assumed to be zero.

2.3 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 \(\varvec{\gamma }\), as defined in (17) for the linearised setting, that incorporates the fibre curvature. A representation of the stored energy density function in terms of \(\mathrm {n_I}\) invariants \(I_i\), which are dependent on \(\{\varvec{\varepsilon },{\varvec{a}},\varvec{\gamma }\}\), is possible, i.e.

$$\begin{aligned} W\left( I_i\left( \varvec{\varepsilon },{\varvec{a}},\varvec{\gamma }\right) \right) . \end{aligned}$$
(31)

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

$$\begin{aligned} \varvec{\sigma }^{\mathrm{{sym}}}= \sum _{n=1}^{\mathrm{{n_I}}} \frac{\partial W}{\partial I_{n}} \frac{\partial I_n}{\partial \varvec{\varepsilon }}. \end{aligned}$$
(32)

Together with the skew-symmetric counterpart, as specified in (26), the total stress tensor follows as

$$\begin{aligned} \varvec{\sigma } = \varvec{\sigma }^{\mathrm{{sym}}}+ \varvec{\sigma }^{\mathrm{{skw}}}. \end{aligned}$$
(33)

The couple-stress tensor can be derived from the stored energy density function in a similar manner. It can be expressed as

$$\begin{aligned} \overline{{\varvec{m}}}^\mathrm {t} = \frac{2}{3} \, \varvec{\epsilon } :\! \sum _{n=1}^{\mathrm{{n_I}}} \frac{\partial W}{\partial I_{n}} \! \left[ \varvec{\beta } \otimes \frac{\partial I_n}{\partial \varvec{\gamma }} + \!\left[ \frac{\partial I_n}{\partial \varvec{\gamma }}\right] ^{\!\mathrm {t}} \! \otimes \varvec{\beta } \right] \!. \end{aligned}$$
(34)

Remark 1

Although the vector \(\varvec{\beta }\) 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 \(\overline{{\varvec{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.

3 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.

3.1 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. \({\varvec{b}}\!=\!{\varvec{c}}\!=\!{\varvec{t}}_2\!=\!{\varvec{0}}\). In order to obtain a weak form, the resulting equation is multiplied by a test function \(\varvec{\eta }\) and integrated over the domain \({\mathcal {B}}\). After applying the divergence theorem twice, the weak form results in

$$\begin{aligned} \int _{\partial {\mathcal {B}}} \varvec{\eta } \cdot {\varvec{t}}_1 ~\mathrm{{d}}a - \int _{{\mathcal {B}}} \nabla _{\!\varvec{x}}\varvec{\eta } : \varvec{\sigma }^{\mathrm{{sym}}}~\mathrm{{d}}v - \int _{{\mathcal {B}}} \left[ \left[ -\frac{1}{2}\varvec{\epsilon }\right] : \nabla _{\!\varvec{x}}^2\varvec{\eta }\right] : {\varvec{m}}^\mathrm {t} ~\mathrm{{d}}v = 0. \end{aligned}$$
(35)

The notation \({\nabla _{\!\varvec{x}}^2\,\bullet =\nabla _{\!\varvec{x}}\nabla _{\!\varvec{x}}\,\bullet }\) is employed for the second gradient.

3.2 Discretisation of the weak form

Within the concept of the finite element method, the weak form (35) is discretised by using basis functions \(R^{\bullet }\). The discretisation of the test function and of the kinematic quantities for an element e yields

$$\begin{aligned} \varvec{\eta }^e&= \sum _{A=1}^{\mathrm{{n_{en}}}} \varvec{\eta }^{eA} R^A, \end{aligned}$$
(36)
$$\begin{aligned} {\varvec{h}}^e&= \sum _{A=1}^{\mathrm{{n_{en}}}} {\varvec{u}}^{eA} \otimes \nabla _{\!\varvec{x}}R^A, \end{aligned}$$
(37)
$$\begin{aligned} \varvec{\kappa }^e&= \sum _{A=1}^{\mathrm{{n_{en}}}} {\varvec{u}}^{eA} \otimes \nabla _{\!\varvec{x}}^2R^A, \end{aligned}$$
(38)

wherein \(\mathrm {n}_{\mathrm{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 \({\mathcal {B}}\) by the assembly operator , the discretised weak form of the governing equation results in

(39)

with the internal and external force vectors

(40)

and

(41)

Therein, the total number of elements is denoted as \(\mathrm {n}_{\mathrm{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 \(\varvec{\sigma }^{\mathrm{{sym}}}\) and \(\varvec{\varvec{\varepsilon }}\) as well as \(\overline{{\varvec{m}}}\) and \(\varvec{\kappa }\) 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

(42)

With these formulations at hand, the symmetric part of the stress tensor as well as the couple-stress tensor can be expressed as

(43)

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

(44)

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.

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.

4.1 Knot vectors and NURBS basis functions

B-Splines are defined on a parametric space which is spanned by knot vectors such as

$$\begin{aligned} \Xi = \{\xi _1=\cdots =\xi _{p+1},\xi _{p+2},\ldots , \xi _{\text {n}},\xi _{\text {n}+1}=\cdots =\xi _{\text {n}+p+1} \}. \end{aligned}$$
(45)

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 \(\xi \) and \(\eta \). 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

$$\begin{aligned} N_{i}^0(\xi ) = {\left\{ \begin{array}{ll} 1 &{} \text {for} ~ \xi _i \le \xi < \xi _{i+1},\\ 0 &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$
(46)

For all higher degrees, B-Splines are constructed recursively following

$$\begin{aligned} N_{i}^p(\xi )&= \frac{\xi -\xi _i}{\xi _{i+p}-\xi _i} N_{i}^{p-1}(\xi ) + \frac{\xi _{i+p+1}-\xi }{\xi _{i+p+1}-\xi _{i+1}} N_{i+1}^{p-1}(\xi ). \end{aligned}$$
(47)

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

$$\begin{aligned} R_{i,j}^{p,q}(\xi ,\eta )&= \frac{N_{i}^{p}(\xi ) M_{j}^{q}(\eta ) w_{i,j}}{\sum _{{\hat{i}}=1}^{\mathrm{n}} \sum _{{\hat{j}}=1}^{\mathrm{m}} {N_{{\hat{i}}}^{p}(\xi ) M_{{\hat{j}}}^{q}(\eta ) w_{{\hat{i}},{\hat{j}}}}} =: \frac{{\widetilde{R}}_{i,j}^{p,q}(\xi ,\eta )}{W(\xi ,\eta )}. \end{aligned}$$
(48)

The numbers \(\mathrm{n}\) and \(\mathrm{m}\) denote the number of B-Spline basis functions in the respective parametric directions \(\xi \) and \(\eta \) with the polynomial degrees p and q.

4.2 Derivatives of NURBS

For the analysis of fibre-reinforced composites with fibre bending stiffness, one main focus lies on the higher-order derivatives of the basis functions. In [26], the derivatives of NURBS with respect to the parametric coordinates \(\xi \) and \(\eta \) are derived as

$$\begin{aligned} R_{i,j}^{(\mathrm {K},\mathrm {L})}&= \frac{1}{W}\bigg [ {\widetilde{R}}_{i,j}^{(\mathrm {K},\mathrm {L})}- \sum _{k=1}^{\mathrm{{K}}} \begin{pmatrix} \mathrm {K} \\ k \end{pmatrix} W^{(k,0)}R_{i,j}^{(\mathrm {K}-{ k},\mathrm {L})} - \sum _{l=1}^{\mathrm{{L}}} \begin{pmatrix} \mathrm {L} \\ l \end{pmatrix} W^{(0,l)}R_{i,j}^{(\mathrm {K},\mathrm {L}-{ l})} \nonumber \\&\quad - \sum _{k=1}^{\mathrm{{K}}} \begin{pmatrix} \mathrm {K} \\ k \end{pmatrix} \sum _{l=1}^{\mathrm{{L}}} \begin{pmatrix} \mathrm {L} \\ l \end{pmatrix} W^{(k,l)}R_{i,j}^{(\mathrm {K}-{ k},\mathrm {L}-{ l})} \bigg ]. \end{aligned}$$
(49)

This general formulation yields the \(\mathrm {K}{{\mathrm{th}}}\) (\(\mathrm{L}{{\mathrm{th}}}\)) derivative with respect to \(\xi \) (\(\eta \)) and includes the derivatives

$$\begin{aligned} {\widetilde{R}}_{i,j}^{(\mathrm{K},\mathrm{L})} = w_{i,j}~\frac{\mathrm{d}^K}{\mathrm{d}\xi ^K}N_i^p(\xi )~\frac{\mathrm{d}^L}{\mathrm{d}\eta ^L}M_j^q(\eta ) \end{aligned}$$
(50)

and

$$\begin{aligned} W^{(k,l)} = \sum _{i=1}^{\mathrm{n}} \sum _{j=1}^{\mathrm{m}} w_{i,j}~\frac{\mathrm{d}^{{ k}}}{\mathrm{d}\xi ^{{ k}}} N_i^p(\xi )~\frac{\mathrm{d}^{{ l}}}{\mathrm{d}\eta ^{{ l}}}M_j^q(\eta ). \end{aligned}$$
(51)

The derivatives of the one-dimensional B-Spline functions follow from

$$\begin{aligned} \frac{\mathrm{d}^{{ k}}}{\mathrm{d}\xi ^{{ k}}} N_{i}^p(\xi ) = \frac{p!}{\left[ p-k\right] !} \sum _{j=0}^{k} a_{k,j} N_{i+j}^{p-k}(\xi ) \end{aligned}$$
(52)

with the coefficients

$$\begin{aligned} a_{0,0}&= 1, \nonumber \\ a_{k,0}&= \frac{a_{k-1,0}}{\xi _{i+p-k+1}-\xi _i}, \nonumber \\ a_{k,j}&= \frac{a_{k-1,j}-a_{k-1,j-1}}{\xi _{i+p+j-k+1}-\xi _{i+j}} ~\text {for}~ j=1,\dots ,k-1, \nonumber \\ a_{k,k}&= \frac{-a_{k-1,k-1}}{\xi _{i+p+1}-\xi _{i+k}} \end{aligned}$$
(53)

specified from the knot vector. The derivatives of the B-Spline basis functions corresponding to the second parametric direction \(\eta \) are obtained analogously.

4.3 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.

Fig. 2
figure 2

Coordinate mapping

Full size image

At first, the master element is connected to the parametric space, on which the NURBS basis functions are defined. For an element \(\Omega _e = [\xi _{i},\xi _{i+1}]\otimes [\eta _{j},\eta _{j+1}]\), the corresponding Jacobian matrix is obtained by

$$\begin{aligned} {\varvec{J}}_1 = \begin{bmatrix} \frac{\partial \xi }{\partial {\widetilde{\xi }}} &{} 0 \\ 0 &{} \frac{\partial \eta }{\partial {\widetilde{\eta }}} \end{bmatrix}\! =\! \begin{bmatrix} \frac{1}{2} \left[ \xi _{i+1}-\xi _{i}\right] &{} 0 \\ 0 &{} \frac{1}{2} \left[ \eta _{j+1}-\eta _{j}\right] \end{bmatrix} \end{aligned}$$
(54)

in the two-dimensional case.

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 \({\varvec{B}}\) contains the respective control points. As there are \(\mathrm{{n}_{\mathrm{cp}} = n\!\cdot \!m}\) control points, this matrix takes the dimension \([\mathrm {n}_{\mathrm{cp}} \times n_{\mathrm{dm}}]\) with \(\mathrm n_{dm} = 2\) for the two-dimensional case. Analogously, the NURBS basis functions will be stored in a vector \({\varvec{R}}\) of dimension \(\mathrm {n}_{\mathrm{cp}}\). With these quantities at hand, the mapping between the vectors \({\varvec{x}}= \left[ x \quad y \right] ^{\mathrm{{t}}}\) and \(\varvec{\xi } = \left[ \xi \quad \eta \right] ^{\mathrm{{t}}}\) can be written as

$$\begin{aligned} {\varvec{J}}_2 = \frac{\partial {\varvec{x}}}{\partial \varvec{\xi }} = {\varvec{B}}^\mathrm {t} \cdot \frac{\partial {\varvec{R}}}{\partial \varvec{\xi }}. \end{aligned}$$
(55)

The final form of the Jacobian determinant includes both contributions, accordingly

$$\begin{aligned} \text {det}({\varvec{J}}) = \text {det}({\varvec{J}}_1)\, \text {det}({\varvec{J}}_2). \end{aligned}$$
(56)

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

$$\begin{aligned} \frac{\partial {\varvec{R}}}{\partial {\varvec{x}}}&= \frac{\partial {\varvec{R}}}{\partial \varvec{\xi }} \cdot \frac{\partial \varvec{\xi }}{\partial {\varvec{x}}} \end{aligned}$$
(57)

and

$$\begin{aligned} \frac{\partial ^2 {\varvec{R}}}{\partial {\varvec{x}}\partial {\varvec{x}}}&= \frac{\partial {\varvec{R}}}{\partial \varvec{\xi }} \cdot \frac{\partial ^2 \varvec{\xi }}{\partial {\varvec{x}}\partial {\varvec{x}}} + \left[ \left[ \frac{\partial \varvec{\xi }}{\partial {\varvec{x}}}\right] ^{\mathrm {t}} \cdot \left[ \frac{\partial ^2 {\varvec{R}}}{\partial \varvec{\xi }\partial \varvec{\xi }} \right] ^\mathrm {t} \cdot \frac{\partial \varvec{\xi }}{\partial {\varvec{x}}} \right] ^\mathrm {t}\!. \end{aligned}$$
(58)

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.

4.4 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^\infty \)-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.

4.5 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.

4.5.1 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 \({\overline{\xi }} \in [ \xi _k,\xi _{k+1} )\) are inserted into the knot vector \(\Xi \) 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

$$\begin{aligned}{}[{\varvec{B}}_i]_j = [x_i,y_i] \end{aligned}$$
(59)

be the \(i{\mathrm{th}}\) control point of the \(j{\mathrm{th}}\) 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

$$\begin{aligned} {\varvec{B}}_i^{{\varvec{w}}} = [w_ix_i,w_iy_i,w_i]^\mathrm {t}, ~ i=1,\dots ,\mathrm {n} \end{aligned}$$
(60)

are introduced. On the basis of this definition, the new set of control points and weights is calculated according to

$$\begin{aligned} \overline{{\varvec{B}}}_i^{{\varvec{w}}} =\alpha _{i}{\varvec{B}}_i^{{\varvec{w}}} + \left[ 1-\alpha _{i}\right] {\varvec{B}}_{i-1}^{{\varvec{w}}}, \end{aligned}$$
(61)

where

$$\begin{aligned} \alpha _i= {\left\{ \begin{array}{ll} 1 &{} \text {for}~ 1 \le i \le k-p,\\ \frac{{\overline{\xi }}-\xi _i}{\xi _{i+p}-\xi _i} &{} \text {for}~ k-p+1 \le i \le k,\\ 0 &{} \text {for}~ k+1 \le i \le \text {n}+p+2.\\ \end{array}\right. } \end{aligned}$$
(62)

For a two-dimensional NURBS surface, the algorithm is applied to each of the \(\mathrm {m}\) NURBS curves. The same procedure is followed when inserting a knot into the knot vector for the second parametric direction \(\eta \).

4.5.2 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^{\infty }\)-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 \(\xi = \xi _r\ne \xi _{r+1}\) of initial multiplicity s is removed t times, the new control points in the \(t{\mathrm{th}}\) removal step can be found as

$$\begin{aligned}&\overline{{\varvec{B}}}_{i,t}^{{\varvec{w}}} =\frac{1}{\alpha _{i}} \left[ {\varvec{B}}_{i,t-1}^{{\varvec{w}}} -(1-\alpha _{i}){\varvec{B}}_{i-1,t}^{{\varvec{w}}} \right] \quad \hbox {for}\quad r-p-t+1\le i \le \tfrac{1}{2}\,[2r-p-s-t]\nonumber \\&\hbox {and}\quad \overline{{\varvec{B}}}_{j,t}^{{\varvec{w}}} = \frac{1}{1-\alpha _{j}} \left[ {\varvec{B}}_{j,t-1}^{{\varvec{w}}} -\alpha _{j}{\varvec{B}}_{j+1,t}^{{\varvec{w}}} \right] \quad \text {for}\quad \tfrac{1}{2}\,[2r-p-s+t+1]\le j \le r-s+t-1 \end{aligned}$$
(63)

with the constant coefficients

$$\begin{aligned} \alpha _i = \frac{{\overline{\xi }}-\xi _i}{\xi _{i+p+t}-\xi _i}, \quad \alpha _j = \frac{{\overline{\xi }}-\xi _{j-t+1}}{\xi _{j+p+1}-\xi _{j-t+1}}. \end{aligned}$$
(64)

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].

4.6 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.

$$\begin{aligned}{}[{\varvec{B}}_1^{{\varvec{w}}}]_j = [{\varvec{B}}_{\mathrm{n}}^{{\varvec{w}}}]_j ~\forall \, j. \end{aligned}$$
(65)

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 \({\varvec{B}}_i^{{\varvec{w}}}\) of each curve specified in (60).

In order to achieve \(C^1\)-continuity, a second condition must be fulfilled, namely

$$\begin{aligned} {\varvec{B}}_2^{{\varvec{w}}} - {\varvec{B}}_1^{{\varvec{w}}} = k \left[ {\varvec{B}}_{\mathrm{n}}^{{\varvec{w}}} - {\varvec{B}}_{\mathrm{n-1}}^{{\varvec{w}}} \right] . \end{aligned}$$
(66)

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 \({\varvec{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 \({\varvec{B}}_1\) can be expressed in terms of the neighbouring points according to

$$\begin{aligned} {\varvec{B}}_1 = \frac{1}{2} \left[ {\varvec{B}}_2 + {\varvec{B}}_{\mathrm{n-1}} \right] \end{aligned}$$
(67)

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.

4.7 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 \({\varvec{x}}\), the relation can be denoted as

$$\begin{aligned} {\varvec{x}}(\xi ,\eta ) = {\varvec{R}}(\xi ,\eta ) \cdot {\varvec{B}}. \end{aligned}$$
(68)

Considering the \(\mathrm {n}\) control points \(\widehat{{\varvec{B}}}\) on the outer boundary of the domain only and choosing a number of \(\mathrm {n}_{\mathrm{sp}}\) interpolating points \(\widehat{{\varvec{x}}}\) on the outer mesh surface, the expression is expanded to

$$\begin{aligned} \widehat{{\varvec{x}}} = \widehat{{\varvec{R}}} \cdot \widehat{{\varvec{B}}} \end{aligned}$$
(69)

with the matrices

$$\begin{aligned}&\widehat{{\varvec{x}}}= \begin{bmatrix} x_1 &{} y_1 \\ x_2 &{} y_2\\ \vdots &{} \vdots \\ x_{\mathrm{{n}}_{\mathrm{sp}}} &{} y_{\mathrm{{n}}_{\mathrm{sp}}} \end{bmatrix},\quad \widehat{{\varvec{R}}} = \begin{bmatrix} R_1^1 &{} \dots &{} R_1^{\mathrm{{n}}} \\ \vdots &{} \ddots &{} \vdots \\ R^1_{\mathrm{{n}}_{\mathrm{sp}}} &{} \dots &{} R^{\mathrm{{n}}}_{\mathrm{{n}}_{\mathrm{sp}}} \end{bmatrix} \quad \hbox {and}\quad \widehat{{\varvec{B}}} :=({\varvec{B}})_{\mathrm{{m}}} = \begin{bmatrix} B_x^1 &{} B_y^1 \\ B_x^2 &{} B_y^2 \\ \vdots &{} \vdots \\ B_x^{\mathrm{{n}}} &{} B_y^{\mathrm{{n}}} \end{bmatrix}. \end{aligned}$$
(70)

Since the isogeometric analysis comprises the isoparametric concept, the same structure holds for the displacements \({\varvec{u}}\) that are associated with either the control points or the mesh points, i.e.

$$\begin{aligned} \widehat{{\varvec{u}}}_\text {sp} = \widehat{{\varvec{R}}} \cdot \widehat{{\varvec{u}}}_\text {cp}. \end{aligned}$$
(71)

For the direct solution of the equation system in (71), matrix \(\widehat{{\varvec{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 \(\widehat{{\varvec{u}}}_{\mathrm{cp}}\) are then imposed to the respective control points as inhomogeneous Dirichlet boundary conditions and will lead to the desired displacements \(\widehat{{\varvec{u}}}_\text {sp}\) of the actual boundary points of the body.

5 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.

5.1 Specific form of the energy function

As proposed in [3], an additive decomposition of the stored energy density function according to

$$\begin{aligned} W = W^{\mathrm{iso}} + W^{\lambda } + W^{\varvec{\kappa }} \end{aligned}$$
(72)

is considered. For the isotropic part, resembling the matrix material, a quadratic form in the small strain tensor is assumed, specifically speaking

$$\begin{aligned} W^{\mathrm{iso}} = \frac{1}{2} ~\lambda ~I_1^2 + \mu ~I_2 \end{aligned}$$
(73)

with the invariants

$$\begin{aligned} I_1 = \text {tr}(\varvec{\varepsilon }),\quad I_2 = \text {tr}(\varvec{\varepsilon }^2). \end{aligned}$$
(74)

The second part \(W^{\lambda }\) 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

$$\begin{aligned} W^{\varvec{\kappa }} = c\, I_6 \end{aligned}$$
(75)

with the fibre bending stiffness parameter c, and with the invariant

$$\begin{aligned} I_6 = \left[ \varvec{\gamma } \cdot {\varvec{a}}\right] \cdot \left[ \varvec{\gamma } \cdot {\varvec{a}}\right] \end{aligned}$$
(76)

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

(77)

and

(78)

The resulting expressions are

$$\begin{aligned} \varvec{\sigma }^{\mathrm{{sym}}}= \lambda \,\text {tr}(\varvec{\varepsilon })\,{\varvec{I}} + 2\,\mu \,\varvec{\varepsilon }\end{aligned}$$
(79)

and

(80)

In the present examples, the Lamé constants are specified to \(\lambda = 1.037\times 10^5\, \hbox {N}\,\hbox {mm}^{-2}\) and \(\mu = 4.444\times 10^4\,\hbox {N}\,\hbox {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.

$$\begin{aligned} d_f = \frac{8}{3}c\,. \end{aligned}$$
(81)

In further accordance with [7], the non-dimensional quantity

$$\begin{aligned} \lambda ^* = \frac{d_f}{2\,\mu \,R_i\, [R_a-R_i]} \end{aligned}$$
(82)

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 \(\lambda ^*\) on the stress and couple-stress components is examined for the numerical examples.

5.2 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.

5.2.1 Geometric description

A rectangular beam with a height of \(10\,\hbox {mm}\) and a length of \(100\,\hbox {mm}\) is considered, so that the ratio of its dimensions is 1/10. For the isogeometric model of the beam, the polynomial degree of the underlying NURBS basis functions is chosen to be \(p=q=4\). The corresponding knot vectors \(\Xi = \{\xi _1,\ldots ,\xi _{\text {n}+p+1} \}\) and \({\mathcal {H}} = \{\eta _1,\ldots ,\eta _{\text {m}+q+1} \}\) are the same for both parametric directions and given by

$$\begin{aligned} \Xi = {\mathcal {H}} = \{0,0,0,0,0,1,1,1,1,1\}. \end{aligned}$$
(83)

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.

Fig. 3
figure 3

Control polygon with \(p=4\) and \(\mathrm{{n_{cp}}=25}\)

Full size image
Fig. 4
figure 4

Boundary conditions and refined mesh with \(\mathrm{{n_{cp}}=1456}\) and \(\mathrm{{n_{el}}=1000}\)

Full size image

5.2.2 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 \(\Vert \overline{{\varvec{t}}}\Vert =30\,\hbox {N}\,\hbox {mm}^{-2}\) which is uniformly distributed over the right surface and accounted for in the analysis by means of Neumann boundary conditions.

5.2.3 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 \(\alpha = 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_y^P\) 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\in \left[ 0,10^7\,\hbox {N}\right] \), 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.

Fig. 5
figure 5

Tip displacement in dependence on the fibre bending stiffness parameter c

Full size image

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\,\hbox {N}\). In particular, the fibre angle \(\alpha =\arccos ({\varvec{a}}\cdot {\varvec{e}}_x)\), i.e. the angle between the fibre vector \({\varvec{a}}\) and the horizontal direction, takes values between \(\alpha =0\) and \(\alpha =\pi /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 \(\alpha =\pi /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.

Fig. 6
figure 6

Tip displacement in dependence on the fibre angle \(\alpha \)

Full size image

5.3 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”.

5.3.1 Geometric description

In contrast to the representation with non-rational B-Splines, a circle can be exactly described by using NURBS. The main idea is to select a set of control points and to calculate the corresponding weights by ensuring that the requirement

$$\begin{aligned} x^2+y^2=R^2 \end{aligned}$$
(84)

is fulfilled for every point (x, y) of the circle with radius R. A full circle can be constructed from a quarter or half circle representation by connecting either four quarter or two half circle curves by means of multiply occurrring knots, see [25].

For a polynomial degree of \(p=4\), the knot vector for a complete circle follows as

$$\begin{aligned} \Xi = \{0,0,0,0,0,\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{2},\frac{1}{2}, \frac{1}{2},\frac{1}{2}, \frac{3}{4},\frac{3}{4},\frac{3}{4},\frac{3}{4},1,1,1,1,1\}. \end{aligned}$$
(85)

For the two-dimensional isogeometric analysis, the knot vector for the second parametric direction is added accordingly

$$\begin{aligned} {\mathcal {H}} = \{0,0,0,0,0,1,1,1,1,1\}, \end{aligned}$$
(86)

so that the control polygon in Fig. 7 is obtained. The corresponding control points and weights are listed in “Appendix B”.

Fig. 7
figure 7

Control polygon with \(p=4\) and \(\mathrm {n}_{\mathrm{cp}} = 85\)

Full size image

The physical mesh results from the application of the NURBS basis functions to the control points and is shown in Fig. 8. In this case, the strategies for mesh refinement and increase in continuity have been applied as described in detail in Sect. 4.5. In order to obtain accurate results near the inner and outer radii of the cylinder, a small mesh size has been chosen for these areas.

Fig. 8
figure 8

Boundary conditions and refined mesh with \(\mathrm {n}_{\mathrm{cp}} = 8004\) and \(\mathrm{{n_{el}}=6912}\)

Full size image

5.3.2 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

$$\begin{aligned} u_{\varphi }(r=R_i)&= 0 \end{aligned}$$
(87)
$$\begin{aligned} u_{\varphi }(r=R_a)&= {\overline{u}}_\varphi \end{aligned}$$
(88)
$$\begin{aligned} m_{rz}(r=R_i)&= 0 \end{aligned}$$
(89)
$$\begin{aligned} m_{rz}(r=R_a)&= 0 \end{aligned}$$
(90)

using cylindrical coordinates. The latter two conditions are in accordance with the assumption of a vanishing surface couple-stress from Sect. 3.1.

5.3.3 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\,\hbox {mm}\) and inner radius \(R_i=1.0\,\hbox {mm}\) is discretised by \(\mathrm {n_{el}}=6912\) finite elements. The prescribed azimuthal displacement is \({\overline{u}}_\varphi =0.025\,\hbox {mm}\). The simulation is performed for different values of the non-dimensional parameter \(\lambda ^*\) introduced in (82), namely \({\lambda ^*\in \{0.0,0.005,0.03,0.1,\pi \}}\). 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

$$\begin{aligned} \zeta = \frac{R_a}{R_i}, \end{aligned}$$
(91)

the dimensionless stress contributions follow as

$$\begin{aligned} \varvec{\sigma }^* = \frac{R_i}{{\overline{u}}_\varphi \, \mu } \frac{\zeta ^2-1}{\zeta }\, \varvec{\sigma } \quad \hbox {and}\quad {\varvec{m}}^* = \frac{1}{{\overline{u}}_\varphi \, \mu } \frac{\zeta ^2-1}{\zeta } \, {\varvec{m}}. \end{aligned}$$
(92)

From the dimensionless expressions for the radius and azimuthal displacement

$$\begin{aligned} r^* = \frac{r}{R_i}, \quad u_\varphi ^* = \frac{u_\varphi }{{\overline{u}}_\varphi } \end{aligned}$$
(93)

the respective form of the fibre slope \({[u_\varphi ]'}\) can be derived, so that

$$\begin{aligned}{}[u^*_\varphi ]' :=\frac{\partial u_\varphi ^*}{\partial r^*} = \frac{R_i}{{\overline{u}}_\varphi } \frac{\partial u_\varphi }{\partial r}. \end{aligned}$$
(94)

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 \(\lambda ^*\), which is associated with the fibre bending stiffness, leads to less bending of the fibres. In fact, for large values such as \(\lambda ^*=\pi \), the fibres almost remain straight.

In Table 1, the fibre slope at the inner boundary, i.e. \(r=R_i=1.0\,\hbox {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 \(\lambda ^*\).

Fig. 9
figure 9

Deformation of a radial fibre

Full size image
Table 1 Fibre slope at the inner boundary \(r=R_i\)
Full size table
Fig. 10
figure 10

Distribution of the dimensionless stress component \([\sigma _{r \varphi }^{\mathrm{sym}}]^*\) and couple-stress component \([m_{rz}]^*\)

Full size image
Fig. 11
figure 11

Distribution of the dimensionless stress component \({\overline{t}}_{(r\theta )}\), corresponding to \([\sigma _{r \varphi }^{\mathrm{sym}}]^*\), and couple-stress component \({\overline{m}}_{rz}\), corresponding to \([m_{rz}]^*\). First published in Journal of Mechanics of Materials and Structures in Vol. 6 (2011), No. 1-4, published by Mathematical Sciences Publishers [7]

Full size image

Stress and couple-stress components   For the two-dimensional boundary value problem of pure azimuthal shear, the relevant stress and couple-stress components are \(\sigma _{r\varphi }\), \(\sigma _{\varphi r}\) and \(m_{rz}\). In addition to the components of the total stress tensor, the values for its symmetric part, namely \(\sigma _{r \varphi }^{\mathrm{{sym}}}=\sigma _{\varphi r}^{\mathrm{{sym}}}\), 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_{rz}\) is given in Fig. 12. For \(\lambda ^*=0\), the couple-stress remains zero, because the higher-order terms are not activated. For all other values of \(\lambda ^*\), the couple-stress 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_{rz}\) 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 \(\sigma _{r \varphi }\) and \(\sigma _{\varphi 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 \(\sigma _{r \varphi }\) is quite similar for all fibre bending stiffness values. For \(\sigma _{\varphi r}\) on the contrary, a monotonous decrease is obtained for perfectly flexible fibres and very small values of \(\lambda ^*\), 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.

Fig. 12
figure 12

Distribution of the dimensionless couple-stress component \([m_{rz}]^*\)

Full size image
Fig. 13
figure 13

Distribution of the dimensionless stress components \([\sigma _{r \varphi }]^*\) and \([\sigma _{ \varphi r}]^*\)

Full size image
Fig. 14
figure 14

Distribution of the dimensionless stress components \({\overline{t}}_{r\theta }\) and \({\overline{t}}_{\theta r}\), corresponding to \([\sigma _{r \varphi }]^*\) and \([\sigma _{ \varphi r}]^*\). First published in Journal of Mechanics of Materials and Structures in Vol. 6 (2011), No. 1-4, published by Mathematical Sciences Publishers [7]

Full size image

6 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.