1 Introduction

Fiber optical sensors, such as fiber Bragg gratings embedded into a surrounding matrix material, are often used to measure the strain at the sensor position inside the fiber. Such measurements can be used, for instance, to infer the magnitude of residual stresses in the matrix material. This can be achieved through an appropriate inverse problem, based on a forward deformation simulation. However, the numerical simulation of the stresses and strains inside components with an embedded measurement fiber under mechanical loadings is challenging due to the difference in typical length scales between the fiber diameter and part geometry. A potential remedy is to simulate the matrix material in the absence of the fiber, and to incorporate the effect of the latter only a posteriori.

The first attempt to analytically model the stress transfer from a uniaxially loaded matrix material to an embedded fiber was made by [5], leading to the emergence of the research field referred to as shear-lag theory nowadays (see for instance [17, 19, 20]), with various applications to fiber optical sensors described, e. g., in [13, 14, 34]. Here, the uniaxial stress in fiber direction is related to the shear stress at the fiber matrix interface, which is recognized as the dominating mechanism of stress transfer from the matrix to the fiber material.

In order to deduce the full strain state inside the material surrounding the fiber instead of only two stress components, one can alternatively employ the strain transfer principle (STP) described in [10, 12]. The STP postulates a linear relationship between the strain tensor inside the sensor and the strain tensor of the far field of the surrounding matrix material (i. e., as though there was no fiber present). This linear relationship can be expressed analytically in the form of the strain transfer tensor and it is valid for orthotropic matrix as well as orthotropic fiber materials; see [10]. An extension of this model to coated fibers and temperature differences between the matrix and the fiber is also available in [32]. The predictions of the analytical STP have been confirmed by various experimental works, for instance [11, 16, 33].

The STP yields particularly good results when the material properties of the matrix and the fiber are similar, or when the fiber material is softer than the matrix. In these cases the fiber does not restrain the deformations of the matrix material in fiber direction under a certain load. The strain in fiber direction in the matrix material is transferred verbatim to the strain inside the fiber in fiber direction and vice versa. In case the matrix material is softer than the fiber, however, the fiber itself may restrain deformations of the entire part/matrix material, and the strain in fiber direction no longer carries over verbatim from the matrix material. The magnitude of this effect also depends on the fiber diameter, part dimensions as well as load conditions under consideration.

To overcome this drawback one has to take into consideration the entire geometry of the part, the embedded fiber as well as the load conditions. However, since the fiber diameter is usually small compared to the part dimensions, a fully resolved finite element (FE) model is often impractical. Instead, we propose an extension of the STP. Our method combines the practical benefits of the STP with the improved accuracy of a fully resolved finite element model. In particular, we may also simulate the deformation of the matrix material in the absence of the fiber.

We refer to our proposed approach as the extended STP. We apply the classical STP to deduce all strain components except the strain in fiber direction from a FE model of the matrix material without the fiber. By contrast, the strain in fiber direction is derived from the solution of a modified FE model. The latter is obtained by superimposing the elastic properties of the bulk matrix material with a one-dimensional model of the fiber. This does not require the fiber to be resolved in the computational mesh.

Models for the simulation of embedded fibers were previously developed in the context of reinforcement of concrete structures (see for instance [6, 7, 24, 25]). They often require an alignment of the reinforcement with the mesh. More recent approaches like [4, 9, 31] use a coupling between the background mesh and the one-dimensional fiber. By using interpolation of quantities between the 1d and 3d mesh, no alignment of the fiber to the background mesh is required. In particular [31] requires the solution of a difficult saddle point problem. The saddle point problem is iteratively solved using a penalty method in order to obtain an approximate solution.

The drawback for using non-aligned meshes (so called “smeared elements”) is that their coupling approach is only valid down to a background mesh element size of about the cross-section diameter of the fiber. Smeared elements may also lead to stability issues. In particular, mortar methods as employed by [31], which are also standard for treatment of conformally meshed interfaces, can result in over-constrained and locking behavior, which causes a loss of energy convergence and unstable interfacial fluxes as stated by [26].

Since the focus of this paper is primarily to improve the original strain transfer principle by including the strain in fiber direction, we use a very simple method for obtaining this strain, which does not introduce additional degrees of freedom and requires only the solution of a linear system of equations. However, in contrast to [31], our method requires that edges of the mesh are aligned with the path of the fiber, which results in a more difficult meshing procedure. Our method for computing the strain in fiber direction can be replaced any time, i. e., by one of the methods mentioned above, given their suitability for thin fibers.

The paper is structured as follows. Section 2 states the linear elasticity problem for a component manufactured from matrix material with an embedded fiber. Section 3 introduces the strain transfer principle and recalls the results from the existing analytical theory. Section 4 presents our extension to this theory. In Sect. 5 we derive the variational form of our elasticity model, and Sect. 6 presents detailed numerical results of our extended STP in comparison with the original STP and with fully resolved finite element computations.

Nomenclature We denote by \(A : B\) the double contraction of a rank-4 tensor A with a matrix B, i. e., \((A : B)_{ij} = \sum _k \sum _\ell A_{ijk\ell } B_{k\ell }\). We also use \(A : B\) for the double contraction of two matrices, i. e., \(A : B = \sum _i \sum _j A_{ij} B_{ij}\). Moreover, \(\mathbb {I}\) denotes the identity on rank-4 tensors, i. e., \(\mathbb {I} : A = A\) holds for any matrix A of appropriate dimensions. \(A \otimes B\) denotes the outer product between matrices A and B, i. e., \((A \otimes B)_{ijk\ell } = A_{ij} B_{k\ell }\). Finally \(a \cdot b\) denotes the usual dot product between vectors a and b, i. e., \(a \cdot b = \sum _i a_i b_i\).

2 Linear elasticity with embedded fiber

Let \(\varOmega \subset \mathbb {R}^d\) denote the domain (\(d=3\)) occupied by the part under consideration. Furthermore, let \(\mathbb {C}:\varOmega \rightarrow L({\mathrm{Sym}}(d))\) denote the stiffness tensor field on \(\varOmega \), i. e., for every material point \(x \in \varOmega \), \(\mathbb {C}(x)\) is a linear mapping between symmetric \(d \times d\) strain matrices and symmetric \(d \times d\) stress matrices. The function \(f :\varOmega \rightarrow \mathbb {R}^d\) denotes a force density field, e. g., due to gravity. Let \(U_{\mathrm{ad}}\) denote the set of all admissible displacements of the part which satisfy given boundary conditions such that (2.3) has a unique solution.

The elastic deformation energy of a displacement field \(u \in U_{\mathrm{ad}}\) on the domain without embedded fiber is given by

$$\begin{aligned} E_{\mathrm{matrix}}(u)= \frac{1}{2} \int _{\varOmega } \varepsilon (u) : \mathbb {C} : \varepsilon (u) \,\mathrm{d} x- \int _{\varOmega } f \cdot u \,\mathrm{d} x, \end{aligned}$$
(2.1)

where \(\varepsilon (u)\) denotes the symmetrized displacement gradient

$$\begin{aligned} \varepsilon (u) = \frac{1}{2} \left( \nabla u + (\nabla u)^{\mathrm{T}} \right) . \end{aligned}$$
(2.2)

We define the equilibrium solution, which solves the following minimization problem, as

$$\begin{aligned} u_{\mathrm{nf}}= \min _{u \in U_{\mathrm{ad}}} E_{\mathrm{matrix}}(u), \end{aligned}$$
(2.3)

where the subscript denotes the absence of the fiber. We refer the reader to [3] for an account of the mathematical theory.

Let \(\gamma :[0, L] \rightarrow \varOmega \) denote the arc-length parameterization of a curve which models the center line of the fiber. Consequently, L denotes the total length of the fiber. Let us assume that each point of the fiber \(\gamma (t)\) is tied to the corresponding point in the domain \(\varOmega \), i. e., we do not consider slip between fiber and matrix material. Given a deformation field \(u \in U_{\mathrm{ad}}\) on the domain, the energy of a one-dimensional fiber generally consists of three parts,

$$\begin{aligned} E_{\mathrm{fiber}}(u) = E_{\mathrm{stretch}}(u) + E_{\mathrm{bend}}(u) + E_{\mathrm{twist}}(u), \end{aligned}$$
(2.4)

modelling the stretching energy \(E_{\mathrm{stretch}}(u)\), the bending energy \(E_{\mathrm{bend}}(u)\) and the twisting energy \(E_{\mathrm{twist}}(u)\), respectively. Following [29, 30], the stretching energy is given by

$$\begin{aligned} E_{\mathrm{stretch}}(u) =\frac{1}{2} \int _{\gamma } E A \varepsilon _\gamma (u)^2 \,\mathrm{d} s, \end{aligned}$$
(2.5)

where E denotes the effective elastic modulus of the fiber material defined below, A is the cross-sectional area and \(\varepsilon _\gamma (u)\) is the strain in fiber direction. Similarly, the bending energy is given by

$$\begin{aligned} E_{\mathrm{bend}}(u) = rac{1}{2} \int _{\gamma } E I \kappa _\gamma (u)^2 \,\mathrm{d} s, \end{aligned}$$
(2.6)

where I denotes the area moment of inertia and \(\kappa _\gamma \) denotes the curvature of \(\gamma \) (i. e., the derivative of the bending angle). For the twisting energy we have

$$\begin{aligned} E_{\mathrm{twist}}(u) = \frac{1}{2} \int _{\gamma } G I_T \delta _\gamma (u)^2 \,\mathrm{d} s, \end{aligned}$$
(2.7)

where G denotes the shear modulus, \(I_T\) is the torsion constant for the section and \( \delta _\gamma (u)\) denotes the derivative of the torsion angle of \(\gamma \). We consider only homogeneous fibers for which E, A, I, G and \(I_T\) are constant. However, the curvature \(\kappa _\gamma \) may vary along the fiber.

Usually, the radius of curvature of the fiber, \(1/\kappa _\gamma (u)\), is much larger than the radius R of the fiber itself. Similarly, the length over which the fiber twists by a full turn \(2\pi /\delta _\gamma (u)\) is usually much larger than R. Furthermore, we have \(A \propto R^2\), \(I \propto R^4\) and \(I_T \propto R^4\). This permits us to conclude \(I \kappa _\gamma (u)^2 \ll R^2\) and \(I_T \delta _\gamma (u)^2 \ll R^2\). As a consequence, the bending and twisting energy terms can be neglected compared to the stretching term for most applications involving only small deformations of the domain \(\varOmega \).

Therefore, we neglect the bending and twisting energies and consider

$$\begin{aligned} E_{\mathrm{fiber}}(u) = E_{\mathrm{stretch}}(u) = \frac{1}{2} \int _{\gamma } E A \varepsilon _\gamma (u)^2 \,\mathrm{d} s. \end{aligned}$$
(2.8)

The equilibrium solution, which solves the following minimization problem is defined as

$$\begin{aligned} u_{\mathrm{f}}= \min _{u \in U_{\mathrm{ad}}} E_{\mathrm{matrix}}(u) + E_{\mathrm{fiber}}(u), \end{aligned}$$
(2.9)

where the subscript denotes the presence of the fiber in the deformation energy.

We return to the definition of the effective stretching Young’s modulus of the fiber material. In order to compensate for the existing matrix material in the volume occupied by the fiber, E is defined as

$$\begin{aligned} E = \max \{0, E_{\mathrm{f}}- E_{\mathrm{m}}\}, \end{aligned}$$
(2.10)

where \(E_{\mathrm{f}}\) is the stretching Young’s modulus of the fiber and \(E_{\mathrm{m}}\) is the stretching Young’s modulus of the matrix in the local fiber direction \(v = \dot{\gamma }(s)\). \(E_{\mathrm{m}}\) is calculated from the matrix compliance as

$$\begin{aligned} E_{\mathrm{m}}= (\varPi _{\gamma } : \mathbb {C}^{-1} : \varPi _{\gamma })^{-1}, \end{aligned}$$
(2.11)

where \(\varPi _{\gamma }= v v^{\mathrm{T}}\). Equation (2.11) recovers the classical Young’s modulus for an isotropic matrix material [2, eqns. (27) and (30)]. For the case \(E_{\mathrm{f}}\le E_{\mathrm{m}}\) which indicates soft fiber material, this results in \(u_{\mathrm{f}}= u_{\mathrm{nf}}\). Otherwise, \(u_{\mathrm{f}}\) and \(u_{\mathrm{nf}}\) are different.

The fiber cross-section is assumed to be of circular shape, such that

$$\begin{aligned} A = \pi R^2 \end{aligned}$$

holds, and the strain in fiber direction \(\varepsilon _\gamma (u)\) is given by

$$\begin{aligned} \varepsilon _\gamma (u) = \varPi _{\gamma } : \varepsilon (u). \end{aligned}$$

These relations allow us to evaluate and minimize the total deformation energy in (2.9). Notice that this does not require to resolve the fiber in the computational mesh.

3 Strain transfer principle

The STP states that there exists a linear relationship between the strain at the center line of the fiber and the strain inside the matrix material. In other words, there exists \(T\in L({\mathrm{Sym}}(d))\) such that

$$\begin{aligned} \varepsilon (u_{\mathrm{f}}) \approx T : \varepsilon (u_{\mathrm{nf}}) \quad \text {along } \gamma . \end{aligned}$$
(3.1)

\(T\) depends only on the geometry of the fiber described by its radius R and the material parameters of the matrix and fiber and possibly the fiber orientation if any of the materials are anisotropic. Equation (3.1) is exact only if \(u_{\mathrm{nf}}\) and \(u_{\mathrm{f}}\) are identical outside of the fiber. In general, this is approximately fulfilled if the fiber has only little influence on the overall deformation \(u_{\mathrm{f}}\). This is for instance the case if \(E_{\mathrm{f}}\ll E_{\mathrm{m}}\), i. e., the fiber is a relatively soft inclusion in the matrix material.

Analytical representation There exists an analytical representation of the strain transfer principle by [10], which considers a homogeneous, uncoated orthotropic fiber with elliptic cross-section embedded into a homogeneous orthotropic fiber reinforced composite with coinciding fiber directions. Under the assumptions of perfect bonding between fiber and matrix material, small deformations, and a uniform stress distribution in the fiber, the displacement and stress continuity conditions at the fiber matrix interface are evaluated using expressions from [12] to describe the linearly elastic effect of the elliptic inclusion in the matrix material. Residual strains in the fiber are neglected and the fiber is assumed to have infinite length. This leads to a 1:1 relation between \((\varepsilon (u_{\mathrm{f}}))_1\) and \((\varepsilon (u_{\mathrm{nf}}))_1\), for \(\varepsilon (u_{\mathrm{f}})\) and \(\varepsilon (u_{\mathrm{nf}})\) in Voigt notation, i. e.,

$$\begin{aligned} \varepsilon (u) = \begin{pmatrix} \varepsilon (u)_{11}&\varepsilon (u)_{22}&\varepsilon (u)_{33}&2 \varepsilon (u)_{23}&2 \varepsilon (u)_{13}&2 \varepsilon (u)_{12} \end{pmatrix}^{\mathrm{T}} . \end{aligned}$$

We assume here that the first axis of the coordinate system is aligned with the fiber direction. We additionally assume constant temperatures, i. e., \(\varDelta T=0\). The detailed derivation in [10] then results in a sparse strain transfer matrix T with the non-zero entries expressed by the relations

$$\begin{aligned}&(\varepsilon (u_{\mathrm{f}}))_1 = (\varepsilon (u_{\mathrm{nf}}))_1, \end{aligned}$$
(3.2a)
$$\begin{aligned}&\begin{pmatrix} (\varepsilon (u_{\mathrm{f}}))_2 \\ (\varepsilon (u_{\mathrm{f}}))_3 \\ (\varepsilon (u_{\mathrm{f}}))_4 \\ \varTheta (u_{\mathrm{f}}) \end{pmatrix}= \big ({O-UL^{-1}N}^{-1}\big )\nonumber \\&\quad \left[ U L^{-1} W \begin{pmatrix} \big ({C_{21}^{\mathrm{f}}- C_{21}}\big ) (\varepsilon (u_{\mathrm{nf}}))_1\\ \big ({C_{31}^{\mathrm{f}}- C_{31}}\big ) (\varepsilon (u_{\mathrm{nf}}))_1\\ 0 \end{pmatrix}\right. \nonumber \\&\qquad \left. +\big ({K-U L^{-1} W Q}\big ) \begin{pmatrix} (\varepsilon (u_{\mathrm{nf}}))_2\\ (\varepsilon (u_{\mathrm{nf}}))_3\\ (\varepsilon (u_{\mathrm{nf}}))_4 \end{pmatrix}\right] , \end{aligned}$$
(3.2b)
$$\begin{aligned}&(\varepsilon (u_{\mathrm{f}}))_5 = \frac{b+\sqrt{\frac{C_{55}}{C_{66}}} a}{b+\frac{C^{\mathrm{f}}_{55}}{\sqrt{C_{66}C_{55}}}a} (\varepsilon (u_{\mathrm{nf}}))_5, \end{aligned}$$
(3.2c)
$$\begin{aligned}&(\varepsilon (u_{\mathrm{f}}))_6 = \frac{a+\sqrt{\frac{C_{66}}{C_{55}}} b}{a+\frac{C^{\mathrm{f}}_{66}}{\sqrt{C_{66}C_{55}}}b} (\varepsilon (u_{\mathrm{nf}}))_6 . \end{aligned}$$
(3.2d)

Here \(\varTheta (u_{\mathrm{f}})\) is the angular displacement of the sensor, which we ignore here. Moreover, we denote the entries of the stiffness tensor \(\mathbb {C}\) in Voigt notation by \(C_{ij}\), \(1 \le i,j \le 6\) and, similarly, the entries of the compliance tensor \(\mathbb {S}=\mathbb {C}^{-1}\) in Voigt notation by \(S_{ij}\), \(1 \le i,j \le 6\). The superscript \(\cdot \,^{\mathrm{f}}\) denotes material parameters of the fiber and no superscript denotes matrix material parameters. Furthermore, \(a = b = R\) denotes the lengths of the semi-axes of the fiber’s cross-section. The matrices in (3.2b) are given by

$$\begin{aligned} K= & {} \begin{pmatrix} a &{} 0 &{} 0\\ 0 &{} 0 &{} \frac{b}{2}\\ 0 &{} 0 &{} \frac{a}{2}\\ 0 &{} b &{} 0 \end{pmatrix} , \quad H = \begin{pmatrix} 0\\ b\\ -a\\ 0 \end{pmatrix} , \quad W = \begin{pmatrix} 1 &{} 0 &{} 0\\ 0 &{} 1 &{} 0\\ 0 &{} 0 &{} 1\\ 0 &{} 0 &{} 1 \end{pmatrix}, \\ Q= & {} \begin{pmatrix} C_{22} &{} C_{23} &{} 0\\ C_{32} &{} C_{33} &{} 0\\ 0 &{} 0 &{} C_{44} \end{pmatrix} , \quad U = \begin{pmatrix} \delta _1 &{} -\delta _2 &{} \delta _1 &{} \delta _2 \\ \delta _2 &{} \delta _1 &{} -\delta _2 &{} \delta _1 \\ \delta _3 &{} -\delta _4 &{} -\delta _3 &{} -\delta _4 \\ \delta _4 &{} \delta _3 &{} \delta _4 &{} -\delta _3 \end{pmatrix}, \end{aligned}$$

with

$$\begin{aligned} \begin{aligned} \delta _1&= 2 \beta _{23} + 2 \beta _{22} (\mu _R^2 - \mu _I^2) , \quad \delta _2 = 4 \beta _{22} \mu _R \mu _I , \quad \\ \delta _3&= 2 \mu _R \left( \beta _{23} + \frac{\beta _{33}}{\mu _R^2+\mu _I^2} \right) ,\\ \quad \delta _4&= 2 \mu _I \left( \beta _{23} - \frac{\beta _{33}}{\mu _R^2+\mu _I^2} \right) \end{aligned} \end{aligned}$$

and

$$\begin{aligned}&L = \begin{pmatrix} \frac{2 \mu _I}{b} &{} \frac{2 \mu _R}{b} &{} \frac{2 \mu _I}{b} &{} -\frac{2 \mu _R}{b}\\ \frac{2}{a} &{} 0 &{} \frac{2}{a} &{} 0\\ 0 &{} -\frac{2}{b} &{} 0 &{} -\frac{2}{b}\\ -\frac{2 \mu _R}{a} &{} \frac{2 \mu _I}{a} &{} \frac{2 \mu _R}{a} &{} \frac{2 \mu _I}{a} \end{pmatrix}, \quad O=\begin{pmatrix} K&H \end{pmatrix},\\&N=\begin{pmatrix}WQ_{\mathrm{f}}&\varvec{0} \end{pmatrix}, \end{aligned}$$

where \({\varvec{0}}\) denotes the zero vector \({\varvec{0}}\in \mathbb {R}^{4 \times 1}\). The matrix \(Q_{\mathrm{f}}\) is similar to Q but with fiber material parameters \(C^{\mathrm{f}}_{ij}\) instead of \(C_{ij}\). Finally, by [12], the parameters \(\beta _{ij}, \mu _R\) and \(\mu _I\) are related to the entries of \(\mathbb {S}\) via

$$\begin{aligned} \begin{aligned} \beta _{ij}&= S_{ij} - \frac{S_{i1}S_{j1}}{S_{11}} , \quad \mu _R = \sqrt{\sqrt{\frac{\beta _{33}}{4 \beta _{22}}} - \frac{2 \beta _{23} + \beta _{44}}{4 \beta _{22}}} , \quad \\ \mu _I&= \sqrt{\sqrt{\frac{\beta _{33}}{4 \beta _{22}}} + \frac{2 \beta _{23} + \beta _{44}}{4 \beta _{22}}} . \end{aligned} \end{aligned}$$

4 Extension of the strain transfer principle

We primarily consider the case \(E_{\mathrm{f}}> E_{\mathrm{m}}\), which is relevant for fiber Bragg grating applications. A typical case is that of a glass fiber embedded in either an isotropic or fiber reinforced plastic. For the sake of simplicity, we ignore here the soft protective coating of the fiber. We propose the following extension of the STP (in tensor notation)

$$\begin{aligned} \varepsilon (u_{\mathrm{f}}) \approx (\mathbb {I}- \varPi _{\gamma }\otimes \varPi _{\gamma }) : \left( T: \varepsilon (u_{\mathrm{nf}}) \right) + \varPi _{\gamma }\varepsilon _\gamma (u_{\mathrm{f}}) \quad \text {on } \gamma .\nonumber \\ \end{aligned}$$
(4.1)

Equation (4.1) lets us recover the strain tensor \(\varepsilon (u_{\mathrm{f}})\) using the strain field \(\varepsilon (u_{\mathrm{nf}})\) computed from the solution \(u_{\mathrm{nf}}\) of (2.3) (without a fiber), and using the strain in fiber direction \(\varepsilon _\gamma (u_{\mathrm{f}})\) computed from the solution \(u_{\mathrm{f}}\) of (2.9) (taking the stiffening due to the fiber into account). Note that since we consider a one-dimensional fiber, the full \(\varepsilon (u_{\mathrm{f}})\) tensor is not directly computable from \(u_{\mathrm{f}}\) in a meaningful way, but only the component \(\varepsilon _\gamma (u_{\mathrm{f}})\) in fiber direction is available. We remark that (4.1) is equivalent to the original STP, but with a corrected strain in fiber direction.

5 Variational formulation and numerical discretization

Problem (2.9) will be solved using the method of finite elements. The displacement \(u \in U_{\mathrm{ad}}\) minimizing the energy in (2.9) is characterized by the variational formulation

$$\begin{aligned}&\int _{\varOmega } \varepsilon (u) : \mathbb {C} : \varepsilon (v) \,\mathrm{d} x+ \int _{\gamma } E A \, \varepsilon _\gamma (u) \, \varepsilon _\gamma (v) \,\mathrm{d} s\nonumber \\&\quad = \int _{\varOmega } f \cdot v \,\mathrm{d} x\quad \text {for all } v \in V, \end{aligned}$$
(5.1)

where \(U_{\mathrm{ad}}= \{u \in H^1(\varOmega )|u = u_0 \text { on } \varGamma _D\}\) and \(\varGamma _D\) is the boundary of \(\varOmega \) with imposed Dirichlet boundary conditions \(u = u_0\) for the displacement. The corresponding test space V is given by \(V = \{u \in H^1(\varOmega )|u = 0 \text { on } \varGamma _D\}\). For the finite element discetization, the domain \(\varOmega \subset \mathbb {R}^3\) is approximated by a tetrahedral mesh consisting of a set of tetrahedrons \(\mathcal {T}(\varOmega )\) such that the fiber \(\gamma \) is approximated by a set of edges \(\mathcal {E}(\gamma )\) of the mesh. The set of admissible displacements \(U_{\mathrm{ad}}\) is approximated by functions \(U_{{\mathrm{ad}},h}\) which are piecewise linear on each tetrahedron, globally continuous and satisfy the Dirichlet boundary conditions. Similarly, the set of test functions is approximated by functions \(V_h\) which are piecewise linear on each tetrahedron, globally continuous and are zero on the boundary \(\varGamma _D\). Due to the linearity, \(\varepsilon (u)\) and \(\varepsilon (v)\) are constant on each tetrahedron. In this discrete setting, the variational form (5.1) becomes

$$\begin{aligned}&\sum _{t \in \mathcal {T}(\varOmega )} \int _{t} \varepsilon (u) : \mathbb {C} : \varepsilon (v) \,\mathrm{d} x+ \sum _{e \in \mathcal {E}(\gamma )} \int _{e} E A \varepsilon _\gamma (u) \varepsilon _\gamma (v) \,\mathrm{d} s\nonumber \\&\quad = \sum _{t \in \mathcal {T}(\varOmega )} \int _{t} f \cdot v \,\mathrm{d} x\qquad \text {for all } v \in V_h . \end{aligned}$$
(5.2)

For an edge \(e \in \mathcal {E}(\gamma )\) and its incident vertices \(p_1\) and \(p_2\), the strain in fiber direction can be calculated as

$$\begin{aligned} \varepsilon _\gamma (u) = \frac{(u(p_2) - u(p_1)) \cdot (p_2 - p_1)}{|p_2 - p_1|^2}, \end{aligned}$$
(5.3)

where \(u(p_1)\) and \(u(p_2)\) denote the nodal displacements in \(p_1\) and \(p_2\).

6 Numerical demonstration

We test the proposed extension of the STP on an example of intermediate complexity. The geometry is a \({100}\,\mathrm{mm} \times {100}\,\mathrm{mm} \times {10}\,\mathrm{mm}\) plate with a \({32}\,\mathrm{mm}\) diameter bore located at \((x,y) = ({60}\,\mathrm{mm}, {40}\,\mathrm{mm})\), as shown in Fig. 1a. As matrix material we use a fiber reinforced plastic, for which we compute the effective stiffness tensor by homogenization, described below in Sect. 6.2. The sensor fiber is made of glass. All FE computations were performed using Dolfin/FEniCS 2019.1; see [1, 15]. The strain transfer described in (3.2) is also computed numerically. Further details are given in the following sections.

Fig. 1
figure 1

Plate with hole (a), path of the sensor fiber (b) and generated mesh (c), (d)

Full size image

6.1 Domain and mesh generation

The geometry shown in Fig. 1a was created using FreeCAD.Footnote 1 Notice that this geometry also contains a fiber of diameter \({4}\mathrm{mm}\) and length \(s= {170}\,{\mathrm{mm}}\) along the path shown in Fig. 1b, which begins at \((x,y,z) = ({-10}\mathrm{mm}, {20}\,{\mathrm{mm}}, {5}\,{\mathrm{mm}})\) and ends at \((x,y,z) = ({80}\,{\mathrm{mm}}, {110}\,{\mathrm{mm}}, {5}\,{\mathrm{mm}})\). The sole purpose of resolving the relatively thick fiber in the mesh is to create a reference finite element solution to compare to the results obtained by the STP. The fiber cross-section itself was split into 4 quadrants such that after meshing there will be an edge following the center line of the fiber. For meshing, the geometry was exported from FreeCAD as a STEP file. The STEP file was then loaded into Gmsh using the OpenCASCADE plugin. The boundaries of the different regions were associated using the “Coherence” function in Gmsh and different subdomains, boundaries and the fiber center line were labeled. The characteristic length used for meshing was calculated from curvature, such that the mesh is more refined near the fiber. The resulting mesh at the boundary and inside of the domain can be seen in Fig. 1c and d, respectively. The mesh contains 161, 476 nodes and 963, 272 tetrahedral elements. The Gmsh mesh including subdomains, boundaries and paths was then loaded and converted to the XDMF format using meshio.Footnote 2

As was mentioned above, the mesh with the three-dimensional fiber resolved is used for the purpose of computing reference solutions. However, the same mesh was also used when computing the strains for the embedded one-dimensional fiber using our extended STP and also for the solution without an embedded fiber. In these cases, the material inside of the meshed fiber was set equal to the matrix material. We followed this procedure to avoid the influence of different meshes on the solutions. In practice, there would be no need to refine the mesh close to the fiber, nor to resolve the fiber in the mesh.

6.2 Homogenization of matrix material

For the demonstration, we used glass fiber reinforced polypropylene as matrix material. For glass we use a Young’s modulus of \(E_g = {73}\,{\mathrm{GPa}}\) and a Poisson ratio of \(\nu _g = 0.18\). For polypropylene we use \(E_{pp} = {1.665}\,{\mathrm{GPa}}\) and \(\nu _{pp} = 0.36\). The fibers are assumed to be parallel and of infinite length in z-direction with a diameter which results in a fiber volume fraction of \(50\%\). The effective stiffness matrix was computed using fibergen from [21], which employs a Lippmann-Schwinger approach (see [18]) with a conjugate gradient solver on a staggered grid described in [8, 28]. Using a laminate mixing rule at the interfaces, see [27], only a resolution of \(55 \times 32 \times 1\) voxels is required for the representative volume element to achieve a sufficient accuracy. The method also allows the computation of effective material properties for other fiber distributions, e. g., for injection molded components. For our FE calculations, we assume the reinforcement fibers to be oriented in the x-direction of the plate. In this instance, one has to swap the x- with the z-axis of the homogenized matrix material stiffness returned by fibergen, which in Voigt notation reads

$$\begin{aligned} \mathbb {C}= \begin{pmatrix} 38.61&{}\quad 2.43 &{}\quad 2.43 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 2.43 &{}\quad 6.34 &{}\quad 3.03 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 2.43 &{}\quad 3.03 &{}\quad 6.34 &{}\quad 0 &{} \quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 1.65 &{} \quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 1.75 &{} 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 1.75 \end{pmatrix} {\mathrm{GPa}}\nonumber \\ \end{aligned}$$
(6.1)

and enters the computation of the strain transfer matrix in the following section.

Fig. 2
figure 2

Diagonal strain components of the solutions along the path of the sensor fiber. The plots on the left show the results for a matrix fiber volume fraction of \(0\%\). The plots on the right are for a matrix fiber volume fraction of \(50\%\)

Full size image

6.3 Computation of strain transfer matrices

For the embedded fiber, we assume the same properties of glass as above, i. e., \(E_{\mathrm{f}}= {73}\,{{\mathrm{GPa}}}\) and a Poisson ratio of \(\nu _{\mathrm{f}}= 0.18\). The strain transfer matrix between strain tensors in Voigt notation representing the respective strain transfer tensor T in (3.1) can then be computed analytically as demonstrated in Sect. 3. However this approach is limited to the case where the fiber axis coincides with one of the axes of orthotropy of the surrounding material. Alternatively, we can evaluate it numerically, which does not impose this limitation. Therefore we employed the numerical approach using fibergen [21] in a similar fashion as for the homogenization. As representative volume element for the latter, we chose a \(1 \times 1\) box with a disc of diameter 0.05 placed at the center representing the fiber, while the remaining domain represents the matrix material. Due to use of periodic boundary conditions for the displacements, the diameter of the disc has to be sufficiently small and the resolution sufficiently large (in our case \(512 \times 512\) voxels). For the identification of the strain transfer matrix, six linearly independent load cases with prescribed strain \(E^{(i)}\) are required. The prescribed strain represents the far field or matrix strain at the fiber position. The computed strain field \(\varepsilon ^{(i)}\) for prescribed strain \(E^{(i)}\) is evaluated at the center of the domain to obtain the strain inside the fiber \(\varepsilon ^{(i)}_{\mathrm{f}}= \varepsilon ^{(i)}(0.5, 0.5)\). The transfer matrix T is then given by the relation

$$\begin{aligned} \begin{pmatrix} \varepsilon ^{(1)}_{\mathrm{f}}| \dots | \varepsilon ^{(6)}_{\mathrm{f}}\end{pmatrix} = T \begin{pmatrix} \varepsilon ^{(1)}_{\mathrm{m}}| \dots | \varepsilon ^{(6)}_{\mathrm{m}}\end{pmatrix}, \end{aligned}$$
(6.2)

where the strains in each column are given in Mandel notation (the notation is only important to interpret the numerical values below). In our instance, the fiber is oriented in z-direction, i. e., parallel to the reinforcement fibers of the matrix material. In order to compute the strain transfer matrix for instances where the sensor fiber has an angle \(\alpha \) to the reinforcement fibers we keep the sensor fiber oriented in z-direction but rotate the matrix material around the x-axis. The rotated \(\mathbb {C}^R\) (in full tensor notation) is then given by

$$\begin{aligned} \mathbb {C}^R_{ijk\ell } = R_{im} R_{jn} R_{kp} R_{\ell q} \mathbb {C}_{mnpq}, \end{aligned}$$

where the sum is carried out over all free indices (using Einstein summation) and \(R=R_x\) represents the rotation matrix for rotation by the angle \(\alpha \) around the x-axis

$$\begin{aligned} R_x = \begin{pmatrix} 1 &{} 0 &{} 0 \\ 0 &{} \cos \alpha &{} -\sin \alpha \\ 0 &{} \sin \alpha &{} \cos \alpha \end{pmatrix}. \end{aligned}$$

Similarly as for \(\mathbb {C}\) of our plate, one has to swap the x- with the z-axis of T to obtain the strain transfer matrix in the correct coordinate system for our example. Furthermore, if the fiber orientation (in the \(x-y\)-plane) has an angle \(\beta \ne 0\) to the x-axis one has to rotate T around the z-axis by angle \(\beta \), i. e.,

$$\begin{aligned} T^R_{ijk\ell } = R_{im} R_{jn} R_{kp} R_{\ell q} T_{mnpq}, \end{aligned}$$

where again the sum is carried out over all free indices (Einstein summation) and \(R=R_z\) represents the rotation matrix for a rotation by the angle \(\beta \) around the z-axis

$$\begin{aligned} R_z = \begin{pmatrix} \cos \beta &{} -\sin \beta &{} 0 \\ \sin \beta &{} \cos \beta &{} 0 \\ 0 &{} 0 &{} 1 \end{pmatrix}. \end{aligned}$$

Note that T also has to be blown up to a full 4-tensor and then converted back to a matrix in Mandel notation.

In the first horizontal section of the path of the sensor fiber in our example (see Fig. 1b) we have \(\beta = 0\) and the computed (and properly rotated) strain transfer matrix as defined in (6.2) is given by

$$\begin{aligned} T = \begin{pmatrix} 1.00 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ -\,0.15 &{}\quad 0.10 &{}\quad 0.02 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ -\,0.15 &{}\quad 0.02 &{}\quad 0.10 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.08 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.11 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.11\end{pmatrix}. \end{aligned}$$

In the middle of the arc section at an angle of \(\beta = \frac{\pi }{4}\), the strain transfer matrix is given by

$$\begin{aligned} T = \begin{pmatrix} 0.67 &{}\quad 0.16 &{}\quad 0.01 &{}\quad 0 &{}\quad 0 &{}\quad 0.25 \\ 0.03 &{}\quad 0.30 &{}\quad 0.01 &{}\quad 0 &{}\quad 0 &{}\quad 0.27 \\ -\,0.11 &{}\quad -\,0.06 &{}\quad 0.11 &{}\quad 0 &{}\quad 0 &{}\quad -\,0.12 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.10 &{}\quad 0.03 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.03 &{}\quad 0.13 &{}\quad 0 \\ 0.21 &{}\quad 0.38 &{}\quad -\,0.01 &{}\quad 0 &{}\quad 0 &{}\quad 0.63\end{pmatrix} \end{aligned}$$

and finally in the vertical section (\(\beta = \frac{\pi }{2}\)) we have

$$\begin{aligned} T = \begin{pmatrix} 0.55 &{}\quad -\,0.14 &{}\quad 0.02 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 1.00 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ -\,0.08 &{}\quad -\,0.14 &{}\quad 0.11 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.10 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.11 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0.11\end{pmatrix}. \end{aligned}$$

Note that in all sections the strain in fiber direction is always transferred verbatim from the matrix material, as recognized, e. g., from the unit diagonal entry in the first and the last fiber sections.

Fig. 3
figure 3

Mixed strain components of the solutions along the path of the sensor fiber. The plots on the left show the results for a matrix fiber volume fraction of \(0\%\). The plots on the right are for a matrix fiber volume fraction of \(50\%\)

Full size image

6.4 Solution of linear elasticity problems

The discretization and solution of (5.1), its counterpart coming from (2.1), and the reference solution are performed using Dolfin/FEniCS  2019.1; see [1, 15]. Meshes, subdomains and the fiber path are loaded from the XDMF files generated as described in Sect. 6.1. The plate is clamped on the lower boundary (\(y={0}\,{\mathrm{mm}}\), see Fig. 1a) and a fixed displacement of \((1,0,0)^{\mathrm{T}}\) mm is enforced on the upper (\(y={100}\,{\mathrm{mm}}\)) boundary. The volume force f is set to zero. For the solution of the arising linear systems we use the conjugate gradient method together with an AMG preconditioner from the Dolfin PETSc backend. Three solutions are obtained: the displacement field of the reference solution \(u_{\mathrm{r}}\) (with three-dimensionally resolved fiber); the solution \(u_{\mathrm{nf}}\) with the fiber neglected by setting the material inside the fiber subdomain to the matrix material; and the solution \(u_{\mathrm{f}}\) from the superimposed one-dimensional fiber model (5.1), where the material inside the fiber subdomain is also set to the matrix material but the fiber’s stiffness enters through the stretching energy term.

6.5 Evaluation

In Figs. 2 and 3 we plot all components of the corresponding strain tensors along the path of the fiber. The strain for the reference solution \(u_{\mathrm{r}}\) is denoted by \(\varepsilon ^{\mathrm{r}}\), the strain for \(u_{\mathrm{nf}}\) is given by \(\varepsilon ^{\mathrm{nf}}\) and the strain \(\varepsilon ^{\mathrm{f}}\) for \(u_{\mathrm{f}}\) is obtained from our extended STP given in (4.1). For comparison, the strain obtained by the original STP (3.1) is denoted by \(T\varepsilon ^{\mathrm{nf}}\). In addition to a matrix fiber volume fraction of \(50\%\) (right plots of Figs. 2 and 3), we also performed the same simulations with a fiber volume fraction of \(0\%\) (left plots of Figs. 2 and 3) representing a very soft and isotropic matrix material.

Ideally, the STP solutions \(T\varepsilon ^{\mathrm{nf}}\) (red triangles) and \(\varepsilon ^{\mathrm{f}}\) (blue diamonds) should be identical to the reference solution \(\varepsilon ^{\mathrm{r}}\) (black stars). For the 11- and 22-components we observe that our extended version of the STP agrees very well with the reference solution, whereas the original STP has major deviations for the 11-component until after the bend of the sensor fiber at around \(s={90}\,{\mathrm{mm}}\) as well as for the 22-component beginning with the bend at around \(s={60}\,{\mathrm{mm}}\). This observation is independent of the matrix fiber volume fraction. For the 33-component the extended STP and original STP are identical, since the sensor fiber direction is always orthogonal to the z-direction. Also they both disagree with the reference solution by a significant amount for the \(0\%\) matrix fiber volume fraction case. The 23- and 13-components are one order of magnitude smaller than the other components and in theory they should be zero in view of the symmetry of the problem in z-direction (one would expect a sign change of displacements at the symmetry plane of the geometry \(z=0\) resulting in a constant zero displacement in z-direction in that plane and therefore resulting in zero strain 23- and 13-components). The reason that the 23- and 13-components are not zero originates from discretization error due to a non-symmetric mesh. Note that the 33-component is generally not expected to be zero in the symmetry plane \(z=0\), even when the z-displacement is zero. Finally, for the 12-component the original and extended STP agree everywhere except for the bend around the hole. Here again the extended STP outperforms the original STP in the case of a \(0\%\) matrix fiber volume fraction. The \(50\%\) case is indecisive.

7 Conclusion and outlook

In this paper, we proposed an improvement of the original strain transfer principle, which recovers the strain components inside a fiber embedded in a matrix material from simulations which do not require the fiber geometry to be resolved. The matrix material itself can be isotropic or fiber reinforced. The proposed modification to the classical STP consists of an additional term in the elastic energy of the total part, which takes into account the additional stretching energy using a simple one-dimensional fiber model. This modification is particularly relevant for fiber materials which are stiffer than the matrix, as it is often the case for fiber optical strain sensors made of glass. Our evaluation shows that the extended STP improves the classical STP in regions in which the presence of the sensor fiber restricts the displacement of the surrounding material in fiber direction.

It is a limitation that the one-dimensional fiber model (2.9) does not incorporate lateral strains, which results in moderate deviations in the 33-component. More deviations are expected due to the neglection of bending and twisting terms. While these have only a minor contribution to the overall energy of the fiber in the chosen example, they may become more relevant in other setups.

The inclusions of bending and twisting energies as well as the consideration of lateral strains are left to future research. We expect that these terms will be quite challenging to model, discretize and implement. Furthermore, as glass fiber sensors are usually coated with a protective layer made of a soft material, a model for coated fibers should be considered. Earlier investigations of [13, 14, 23] suggest that a low modulus coating layer absorbs some portion of the matrix strain in the transfer process, leading to a relative decrease in fiber strains compared to uncoated fibers. For this case there already exists an analytical STP proposed in the literature; see for instance [32]. Additionally, the integration of ideas from the shear-lag theory mentioned in the introduction might prove beneficial to address the usually vast differences in the material properties between fiber and coating materials.

Supplemental material for this publication can be found at https://github.com/fospald/strain-transfer-principle and is citeable via [22].