Modeling the non-trivial behavior of anisotropic beams: a simple Timoshenko beam with enhanced stress recovery and constitutive relations

TL;DR

This paper develops an enhanced Timoshenko beam model that accurately captures the complex effects of material anisotropy on structural behavior, including stress distribution and displacements, validated against finite element solutions.

Contribution

It introduces a stress recovery method and consistent constitutive relations for anisotropic beams, extending analytical capabilities beyond homogeneous cases.

Findings

Axial stress depends significantly on transversal internal force.

Anisotropy affects displacements more than shear deformation.

Model achieves errors less than 5% compared to finite element solutions.

Abstract

This paper analyzes the non-trivial influence of the material anisotropy on the structural behavior of an anisotropic multilayer planar beam. Indeed, analytical results available in literature are limited to homogeneous beams and several aspects has not been addressed yet, impeding an in-depth understanding of the mechanical response of anisotropic structural elements. This paper proposes an effective recovery of stress distribution and an energetically consistent evaluation of constitutive relations to be used within a planar Timoshenko beam model. The resulting structural-analysis tool highlights the following peculiarities of anisotropic beams: (i) the axial stress explicitly depends on transversal internal force, which can weigh up to 30% on the maximal magnitude of axial stress, and (ii) the anisotropy influences the beam displacements more than standard shear deformation and even for extremely slender beams. A rigorous comparison with analytical and accurate 2D Finite Element solutions confirms the accuracy of the proposed approach that leads to errors exceptionally greater than 5%.

Tables

Table 1. Table 1: Bi-layer anisotropic cantilever. Maximal axial displacement u ​ ( l ) 𝑢 𝑙 u\left(l\right) evaluated according to EB u E ​ B superscript 𝑢 𝐸 𝐵 u^{EB} and proposed u m ​ o ​ d superscript 𝑢 𝑚 𝑜 𝑑 u^{mod} beam models and relative errors.

$\lambda$	$\theta \left[\mathrm{deg}\right]$	$u^{ref}\left[\mathrm{mm}\right]$	$u^{EB}\left[\mathrm{mm}\right]$	$u^{mod}\left[\mathrm{mm}\right]$	$e_u^{EB}[\%]$	$e_u^{mod}[\%]$
5	+15	1.173e–1	0.000e+0	1.078e–1	100	8.10
5	–15	–1.125e–1	0.000e+0	–1.078e–1	100	4.18
10	+15	4.612e–1	0.000e+0	4.311e–1	100	6.53
10	–15	–4.504e–1	0.000e+0	–4.311e–1	100	4.29
20	+15	1.828e+0	0.000e+0	1.724e+0	100	5.69
20	–15	–1.803e+0	0.000e+0	–1.724e+0	100	4.38

Table 2. Table 2: Bi-layer anisotropic cantilever. Maximal rotations ϕ ​ ( l ) italic-ϕ 𝑙 \phi\left(l\right) evaluated according to EB ϕ E ​ B superscript italic-ϕ 𝐸 𝐵 \phi^{EB} and proposed ϕ m ​ o ​ d superscript italic-ϕ 𝑚 𝑜 𝑑 \phi^{mod} beam models and relative errors.

$\lambda$	$\theta \left[\mathrm{deg}\right]$	${\varphi }^{ref}\left[\text{rad}\right]$	${\varphi }^{EB}\left[\text{rad}\right]$	${\varphi }^{mod}\left[\text{rad}\right]$	$e_{\varphi }^{EB}[\%]$	$e_{\varphi }^{mod}[\%]$
5	+15	–3.569e–2	–3.189e–2	–3.513e–2	10.6	1.57
5	–15	–2.841e–2	–3.189e–2	–2.864e–2	12.2	0.81
10	+15	–2.706e–1	–2.551e–1	–2.681e–1	5.73	0.92
10	–15	–2.410e–1	–2.551e–1	–2.421e–1	5.85	0.46
20	+15	–2.103e+0	–2.041e+0	–2.093e+0	2.95	0.48
20	–15	–1.984e+0	–2.041e+0	–1.989e+0	2.87	0.25

Table 3. Table 3: Bi-layer anisotropic cantilever. Maximal transversal displacements v ​ ( l ) 𝑣 𝑙 v\left(l\right) evaluated according to EB v E ​ B superscript 𝑣 𝐸 𝐵 v^{EB} , Timoshenko v T superscript 𝑣 𝑇 v^{T} , and proposed v m ​ o ​ d superscript 𝑣 𝑚 𝑜 𝑑 v^{mod} beam models and relative errors.

$\lambda$	$\theta \left[\mathrm{deg}\right]$	$v^{ref}\left[\mathrm{mm}\right]$	$v^{EB}\left[\mathrm{mm}\right]$	$v^{EB}+v^T\left[\mathrm{mm}\right]$	$v^{mod}\left[\mathrm{mm}\right]$	$e_v^{EB}[\%]$	$e_v^{EB+T}[\%]$	$e_v^{mod}[\%]$
5	+15	–1.545e+1	–1.196e+1	–1.364e+1	–1.516e+1	22.6	11.7	1.88
5	–15	–1.177e+1	–1.196e+1	–1.364e+1	–1.191e+1	1.61	15.9	1.19
10	+15	–2.130e+2	–1.913e+2	–1.981e+2	–2.106e+2	10.2	7.00	1.13
10	–15	–1.834e+2	–1.913e+2	–1.981e+2	–1.847e+2	4.31	8.02	0.71
20	+15	–3.210e+3	–3.061e+3	–3.088e+3	–3.190e+3	4.64	3.80	0.62
20	–15	–2.972e+3	–3.061e+3	–3.088e+3	–2.982e+3	2.99	3.90	0.34

Table 4. Table 4: Doubly clamped bi-layer anisotropic beam. Axial constraint reactions N ​ ( 0 ) 𝑁 0 N\left(0\right) and N ​ ( l ) 𝑁 𝑙 N\left(l\right) evaluated according to EB N E ​ B superscript 𝑁 𝐸 𝐵 N^{EB} and proposed N m ​ o ​ d superscript 𝑁 𝑚 𝑜 𝑑 N^{mod} beam models and relative errors.

$\lambda$	$x$	$N^{ref}\left[N\right]$	$N^{EB}\left[N\right]$	$N^{mod}\left[N\right]$	$e_N^{EB}[\%]$	$e_N^{mod}[\%]$
5	0	7.356e+0	0.000e+0	8.742e+0	100	18.8
5	$l$	7.356e+0	0.000e+0	8.742e+0	100	18.8
10	0	8.959e+0	0.000e+0	1.092e+1	100	21.9
10	$l$	8.959e+0	0.000e+0	1.092e+1	100	21.9
20	0	9.472e+0	0.000e+0	1.165e+1	100	23.0
20	$l$	9.472e+0	0.000e+0	1.165e+1	100	23.0

Table 5. Table 5: Doubly clamped bi-layer anisotropic beam. Bending moment constraint reactions M ​ ( 0 ) 𝑀 0 M\left(0\right) and M ​ ( l ) 𝑀 𝑙 M\left(l\right) evaluated according to EB M E ​ B superscript 𝑀 𝐸 𝐵 M^{EB} and proposed M m ​ o ​ d superscript 𝑀 𝑚 𝑜 𝑑 M^{mod} beam models and relative errors.

$\lambda$	$x$	$M^{ref}\left[\mathrm{Nm}\right]$	$M^{EB}\left[\mathrm{Nm}\right]$	$M^{mod}\left[\mathrm{Nm}\right]$	$e_M^{EB}[\%]$	$e_M^{mod}[\%]$
5	0	–1.732e+4	–2.083e+4	–1.794e+4	20.3	3.58
5	$l$	–2.482e+4	–2.083e+4	–2.415e+4	16.1	2.70
10	0	–7.458e+4	–8.333e+4	–7.583e+4	11.7	1.68
10	$l$	–9.244e+4	–8.333e+4	–9.137e+4	9.86	1.16
20	0	–3.138e+5	–3.333e+5	–3.170e+5	6.21	1.02
20	$l$	–3.526e+5	–3.333e+5	–3.502e+5	5.47	0.68

Table 6. Table 6: Doubly clamped bi-layer anisotropic beam. Shear constraint reactions V ​ ( 0 ) 𝑉 0 V\left(0\right) and V ​ ( l ) 𝑉 𝑙 V\left(l\right) evaluated according to EB V E ​ B superscript 𝑉 𝐸 𝐵 V^{EB} and proposed V m ​ o ​ d superscript 𝑉 𝑚 𝑜 𝑑 V^{mod} beam models and relative errors.

$\lambda$	$x$	$V^{ref}\left[N\right]$	$V^{EB}\left[N\right]$	$V^{mod}\left[N\right]$	$e_V^{EB}[\%]$	$e_V^{mod}[\%]$
5	0	–2.333e+2	–2.500e+2	–2.376e+2	7.16	1.84
5	$l$	–2.645e+2	–2.500e+2	–2.624e+2	5.48	0.79
10	0	–4.817e+2	–5.000e+2	–4.845e+2	3.80	0.58
10	$l$	–5.184e+2	–5.000e+2	–5.155e+2	3.55	0.56
20	0	–9.921e+2	–1.000e+3	–9.834e+2	0.80	0.88
20	$l$	–1.026e+2	–1.000e+3	–1.017e+3	2.53	0.88

Equations

$$L:=\left\{x\in \left[0,l\right]\right\};H:=\left\{y\in \left[-\beta h,\left(1-\beta \right)h\right]\right\}$$

$$\Omega :=L\times H$$

$${ϵ}_x\left(x,y\right)=$$

$${ϵ}_y\left(x,y\right)=$$

$${\gamma }_{xy}\left(x,y\right)=$$

$${\sigma }_{x,x}\left(x,y\right)+{\tau }_{,y}\left(x,y\right)=$$

$${\tau }_{,x}\left(x,y\right)+{\sigma }_{y,y}\left(x,y\right)=$$

$${ϵ}_x\left(x,y\right)=$$

$${ϵ}_y\left(x,y\right)=$$

$${\gamma }_{xy}\left(x,y\right)=$$

$$\boldsymbol{D}\left(y\right)=\left[\begin{array}[]{ccc}\frac{1}{E_{xx}\left(y\right)}&\frac{1}{E_{xy}\left(y\right)}&\frac{1}{G_{x}\left(y\right)}\\ \frac{1}{E_{xy}\left(y\right)}&\frac{1}{E_{yy}\left(y\right)}&\frac{1}{G_{y}\left(y\right)}\\ \frac{1}{G_{x}\left(y\right)}&\frac{1}{G_{y}\left(y\right)}&\frac{1}{G\left(y\right)}\end{array}\right]=\boldsymbol{R}^{T}\left(y\right)\left[\begin{array}[]{ccc}\frac{1}{{E}_{11}\left(y\right)}&-\frac{\nu\left(y\right)}{{E}_{11}\left(y\right)}&0\\ -\frac{\nu\left(y\right)}{{E}_{11}\left(y\right)}&\frac{1}{E_{22}\left(y\right)}&0\\ 0&0&\frac{1}{G_{12}\left(y\right)}\end{array}\right]\boldsymbol{R}\left(y\right)$$

$$\boldsymbol{R}\left(y\right)=\left[\begin{array}[]{ccc}\cos^{2}\left(\theta\left(y\right)\right)&\sin^{2}\left(\theta\left(y\right)\right)&2\sin\left(\theta\left(y\right)\right)\cos\left(\theta\left(y\right)\right)\\ \sin^{2}\left(\theta\left(y\right)\right)&\cos^{2}\left(\theta\left(y\right)\right)&-2\sin\left(\theta\left(y\right)\right)\cos\left(\theta\left(y\right)\right)\\ -\sin\left(\theta\left(y\right)\right)\cos\left(\theta\left(y\right)\right)&\sin\left(\theta\left(y\right)\right)\cos\left(\theta\left(y\right)\right)&\cos^{2}\left(\theta\left(y\right)\right)-\sin^{2}\left(\theta\left(y\right)\right)\end{array}\right]$$

$$A^{\ast }={\int }_H{E}_{xx}\left(y\right)dy;\beta =\frac{1}{h{A}^{\ast }}{\int }_H{E}_{xx}\left(y\right)ydy;I^{\ast }={\int }_H{E}_{xx}\left(y\right)y^{2}dy$$

$$s_x\left(x,y\right)\approx$$

$$s_y\left(x,y\right)\approx$$

$$ϵ\left(x\right)$$

$$\chi \left(x\right)$$

$$\gamma \left(x\right)$$

$$N\left(x\right)={\int }_H{\sigma }_x\left(x,y\right)dy;M\left(x\right)={\int }_H{\sigma }_x\left(x,y\right)\left(-y\right)dy;V\left(x\right)={\int }_H\tau \left(x,y\right)dy$$

$$q={\int }_H{f}_{y}dy=h{f}_y$$

$$N^{\prime }\left(x\right)$$

$$M^{\prime }\left(x\right)$$

$$V^{\prime }\left(x\right)$$

$${\sigma }_x\left(x,y\right)=E_{xx}\left(y\right){ϵ}_x\left(x,y\right)-\frac{E_{xx}\left(y\right)}{G_x\left(y\right)}\tau \left(x,y\right)$$

$$\tau \left(x,y\right)=-{\int }_{-\beta h}^y{\sigma }_{x,x}\left(x,\hat{y}\right)d\hat{y}$$

$${\sigma }_x\left(x,y\right)=d_{{\sigma }_x}^N\left(y\right)N\left(x\right)+d_{{\sigma }_x}^M\left(y\right)M\left(x\right)+d_{{\sigma }_x}^V\left(y\right)V\left(x\right)+d_{{\sigma }_x}^q\left(y\right)q$$

$$\tau \left(x,y\right)=d_{\tau }^V\left(y\right)V\left(x\right)+d_{\tau }^q\left(y\right)q$$

$$d_{{\sigma }_x}^N\left(y\right)$$

$$d_{{\sigma }_x}^M\left(y\right)$$

$$d_{{\sigma }_x}^V\left(y\right)$$

$$d_{{\sigma }_x}^q\left(y\right)$$

$$d_{\tau }^V\left(y\right)$$

$$d_{\tau }^q\left(y\right)$$

$${\Psi }^{\ast }\left(x,y\right)=\frac{1}{2}{\mathbf{\sigma }}^T\left(x,y\right)\cdot \mathbf{D}\left(y\right)\cdot \mathbf{\sigma }\left(x,y\right)=\frac{1}{2}\left(\frac{{\sigma }_x^2\left(x,y\right)}{E_{xx}\left(y\right)}+\frac{{\tau }^2\left(x,y\right)}{G\left(y\right)}+2\frac{{\sigma }_x\left(x,y\right)\tau \left(x,y\right)}{G_x\left(y\right)}\right)$$

$$ϵ\left(x\right)=$$

$$\chi \left(x\right)=$$

$$\gamma \left(x\right)=$$

$${ϵ}_N$$

$${ϵ}_M={\chi }_N$$

$${ϵ}_V={\gamma }_N$$

Full text

1 Introduction

Effective analysis and design of timber and composite structures unavoidably require beam and plate models capable to handle the material anisotropy. Nowadays the need of accurate analysis tools is even more urgent due to the fast development of novel technologies. As an example, additive manufacturing allows to create structural elements with variable orientation of fibers [14]. Likewise, laser scanners detect grain orientation on the surfaces of timber boards, allowing for accurate analyses of glued laminated timber beams [21] and advanced optimization of structural elements [30].

Engineering research has spent great effort in terms of beam and plate modeling within the last decades. Nevertheless, most of models were derived under the hypothesis of isotropic or, at most, orthotropic material [36, 8]. As a consequence, several features of anisotropic structural elements are not yet well addressed, in particular when principal directions of the material are not aligned with the beam axis or the plate reference surface.

Limiting the discussion to planar problems, the generalized 2D Hook’s law of an isotropic material is represented by a block diagonal matrix. Consistently, beam constitutive relations are represented by a diagonal matrix i.e., axial deformation, curvature, and shear deformation uniquely depend on axial internal force, bending moment, and transversal internal force, respectively. Conversely, anisotropy leads the generalized 2D Hook’s law to be represented by a full matrix [24]. As a consequence, also constitutive relations for anisotropic structural elements may be represented by a full matrix i.e., all generalized deformations may depend on all internal forces.

Despite its importance, the above-mentioned problematic was only partially addressed in literature. Murakami et al. [28] proposed a First-order Shear Deformation Theory (FSDT) for a planar homogeneous anisotropic beam where a coupling term (mentioned also as coefficient of mutual influence [27, 26]) relates axial deformation with shear force (and vice-versa). In the successive years, several researchers [29, 20, 31, 35, 27] used different approaches for the estimation of the coupling term, reaching slightly different solutions. As will be discussed in the following, the introduction of a single coupling term allows to define models that are effective only for an extremely limited set of cross-section geometries, while simple and effective models capable to handle more general cases has not been proposed jet.

More recently, Karttunen and Von Hertzen [22] have proposed an accurate analysis of the structural behavior of an homogeneous anisotropic planar beam. The analytical expression for the stress distribution is calculated using the Airy’s stress function, analytical expression for deformations is computed using 2D constitutive relations, and the 2D displacement field is recovered using the compatibility Partial Differential Equations. Simple calculations allow to reformulate the obtained analytical expression for 2D displacements and stresses in terms of 1D functions coinciding with internal forces and FSDT kinematic parameters. Taking advantage of this simplification, the authors provide the analytical expression of the Finite Element (FE) stiffness matrix of the beam. Analogous stress distributions was obtained also by Hashin [17], nevertheless the simplification proposed by Karttunen and Von Hertzen [22] highlights that the axial stress explicitly depends on transversal internal force due to the non-trivial constitutive relations of the material. The above discussed analytical results represent a milestone for the development of effective anisotropic beams. Nevertheless, the derivation procedure can not be easily generalized to multilayer structures, resulting therefore of limited interest for practitioners.

An other significant aspect that has to be carefully handled is the length of zones where boundary effects extinguish according to the Saint-Venant principle. While for a isotropic beam boundary effects are negligible at a distance greater than the maximal size of the cross-section, for an anisotropic beams such a distance depends on the ratio between axial and shear modulus and may be greater than six or seven times the maximal cross-section dimension [9, 6, 18]. On the one hand, this reduces the effectiveness of beam models and, on the other hand, it introduces further phenomena to be considered in the analysis of structural elements, impeding a straightforward interpretation of both numerical and experimental results.

Nowadays, effective and accurate cross-section analysis tools (e.g., [15], Variational Asymptotic Beam Sectional Analysis [37, 38, 32, 16], Semi-Analytical Finite Element [12, 23, 25, 11], Generalized Beam Thory [33]) that may accurately handle the so far introduced problems are available. Nevertheless, all the cross-section analysis tools are based on auxiliary PDEs and functionals, impeding an immediate physical understanding of the analysis results. As a consequence, engineers use the above-mentioned analysis tools as black-boxes. Furthermore, the scarce awareness about the effects of anisotropy on the structural behavior leads engineers to erroneously believe that coarse adaptations of isotropic beam models are effective [4].

This paper proposes a simple planar beam model that effectively describes the linear elastic behavior of anisotropic multi-layer structural elements accounting for the previously introduced issues. Specifically, the beam model will assume that the beam cross-section behaves rigidly, in analogy to the Timoshenko beam. On the one hand, this choice limits the accuracy and the applicability of the proposed beam model. On the other hand, it leads to Ordinary Differential Equations for which an analytical solution can be computed and easily interpreted by simple physical considerations, allowing for a deep understanding of the structural behavior of anisotropic beams.

The main novelty of the developed model is an enhanced and effective stress recovery based on a two-steps iterative procedure. The former step uses the first 2D constitutive relation for the recovery of axial stress and allows to handle the effects of the anisotropy on the stress distribution. The latter step uses the horizontal equilibrium PDE for the recovery of shear stress, in analogy with standard Jourawsky approach [19, 7]. Such a procedure allows to identify the non-trivial and explicit dependence of axial stress on transversal internal force and to manage also multi-layer anisotropic beams. Beam constitutive relations are derived from the stress potential and the outcomes of stress recovery procedure. Such an approach properly embeds anisotropy effects within the beam model and the effectiveness of the proposed constitutive relation derivation path was already demonstrated for non-prismatic and functionally graded material beams [3, 1, 2, 5].

Numerical results will demonstrate that the proposed beam model describes the behavior of anisotropic structural element with a good accuracy, leading to an extremely convenient cost-benefit ratio. In particular, the proposed beam model effectively predicts: (i) the highly non-linear distribution of axial stresses, obtained despite deformations have a linear distribution and the material is linear-elastic, (ii) the explicit dependency of horizontal stress on transversal internal force and load, (iii) the fact that the anisotropy influences the beam displacements more than shear deformation.

The outline of the paper is as follows: Section 2 defines the problem and illustrates the beam model ODEs, Section 3 derives the ODEs analytical solution, Sections 4 and 5 discuss some meaningful examples, and Section 6 resumes main properties, advantages, and limitations of the proposed method and delineates future research.

2 Beam model

This section introduces the 2D problem as well as notations (Section 2.1), it discusses the beam compatibility and equilibrium ODEs (Section 2.2), the recovery of cross-section stress distribution (Section 2.3), and the beam constitutive relations (Section 2.4).

2.1 2D problem definition

The beam longitudinal axis $L$ and the beam cross-section $H$ are closed and bounded subsets of $x$- and $y$- axes defined as

[TABLE]

where $l$, $h$, and $0<\beta<1$ are the beam length, the beam thickness, and a dimensionless parameter defining the distance between the $x$-axis and the lower boundary of the cross-section, respectively. The beam depth $b$ denotes the cross-section size along the $z$ coordinate and, in the following, we assume that $b=1$. As illustrated in Figure 1, the 2D beam body $\Omega$ is defined as

[TABLE]

Finally, we assume that the body is slender (i.e., $l\gg h$) and behaves under the hypothesis of plane stress and small displacements.

We introduce the displacement vector field $\boldsymbol{s}\left(x,y\right)=\left[s_{x}\left(x,y\right),s_{y}\left(x,y\right)\right]$, the stress tensor field $\boldsymbol{\sigma}\left(x,y\right)=\left[\sigma_{x}\left(x,y\right),\sigma_{y}\left(x,y\right),\tau\left(x,y\right)\right]^{T}$, and the strain tensor field $\boldsymbol{\varepsilon}\left(x,y\right)=\left[\epsilon_{x}\left(x,y\right),\epsilon_{y}\left(x,y\right),\gamma_{xy}\left(x,y\right)\right]^{T}$ using the engineering notation. Furthermore, a distributed load $\boldsymbol{f}=\left[f_{x},f_{y}\right]$ is applied within the domain and suitable Boundary Conditions are assigned on the boundary of domain $\Omega$. The 2D compatibility PDEs read

[TABLE]

where the notation $\left(\cdot\right),i$ for $i=x,y$ represents partial derivatives. The 2D equilibrium PDEs read

[TABLE]

The beam is made of a linear-elastic and anisotropic material. As represented in Figure 1, the material properties do not depend on the beam axis $x$ coordinate and are piecewise constant within the beam thickness. Following the notation introduced by Murakami et al. [28], the 2D anisotropic constitutive relations can be represented as

[TABLE]

The coefficients of the material constitutive relations can be collected in a matrix $\boldsymbol{D}$ that is defined as

[TABLE]

where $\boldsymbol{R}\left(y\right)$ reads

[TABLE]

${E}_{11}\left(y\right),E_{22}\left(y\right),G_{12}\left(y\right)$, and $\nu\left(y\right)$ are the parameters defining the mechanical properties of the material with respect to the principal directions. The quantity $\theta\left(y\right)$, with $-\pi/2<\theta\left(y\right)<\pi/2$, is the rotation of the principal direction of the material with respect to the $x$-axis.

Remark 2.1

Due to the definition (6), $E_{xx}\left(y\right)$, $E_{xy}\left(y\right)$, $E_{yy}\left(y\right)$, and $G\left(y\right)$ are even functions of the material principal direction rotation $\theta\left(y\right)$ whereas the material coupling terms $G_{x}\left(y\right)$ and $G_{y}\left(y\right)$ are odd.

2.2 Compatibility and equilibrium ODEs

For convenience, we define the axial stiffness $A^{*}$, the dimensionless parameter $\beta$ introduced in Equation (1), and the bending stiffness $I^{*}$

[TABLE]

Remark 2.2

Due to Definition (8), the origin of the adopted Cartesian coordinate system $O$ coincides with the so-called stiffness centroid that is equal to the cross-sectional geometric centroid only if the cross-section is symmetric and, in particular, when the beam is homogeneous, as discussed also by Kosmatka et al. [23].

As usual for standard Timoshenko beam models, the 2D displacement field $\boldsymbol{s}\left(x,y\right)=\left[s_{x}\left(x,y\right),s_{y}\left(x,y\right)\right]^{T}$ is represented in terms of three 1D functions, indicated as axial displacement $u\left(x\right)$, cross-section rotation $\phi\left(x\right)$, and transversal displacement $v\left(x\right)$. Therefore, the displacement field components are approximated as follows

[TABLE]

Introducing the generalized strains defined as the axial strain $\epsilon\left(x\right)$, the curvature $\chi\left(x\right)$, and the shear strain $\gamma\left(x\right)$, the beam compatibility is expressed through the following ODEs

[TABLE]

where the notation $\left(\cdot\right)^{\prime}$ denotes derivatives with respect to $x$.

Remark 2.3

In light of Remark 2.2, kinematic approximation (9) differs from standard Timoshenko one. In fact, $u\left(x\right)$ represents the axial displacement of the stiffness centroids and, in general, it does not coincide with the mean value of the cross-section axial displacements (i.e., $u\left(x\right)=s_{x}\left(x,0\right)\neq 1/h\int_{H}s_{x}\left(x,y\right)dy$). Similarly, $\epsilon\left(x\right)$ is not the mean value of the axial strain evaluated within the cross-section, but just represents the axial elongation evaluated at $y=0$ (i.e., $\epsilon\left(x\right)\neq 1/h\int_{H}\partial s_{x}\left(x,y\right)/\partial x\,dy$).

We introduce the axial internal force $N\left(x\right)$, the bending moment $M\left(x\right)$, and the transversal internal force $V\left(x\right)$ defined as

[TABLE]

Furthermore, we assume that $f_{x}=0$ and we introduce the transversal load $q$

[TABLE]

Considering the axial, rotational, and transversal equilibrium of a infinitesimally long beam-segment, the equilibrium ODEs read

[TABLE]

2.3 Stress recovery

The stress recovery is based on a recursive procedure that leads to define a distribution of stresses that satisfies the first constitutive relation (5a) and the first equilibrium PDE (4a). Conversely, we assume that transversal stress vanishes i.e., $\sigma_{y}\left(x,y\right)=0$, aiming at the maximal simplicity of the model.

In order to set the recursive procedure up, it is convenient to isolate Equations (5a) and (4a) for the variables $\sigma_{x}\left(x,y\right)$ and $\tau\left(x,y\right)$, respectively

[TABLE]

[TABLE]

The iterative procedure is resumed in Figure 2 and leads to the following distribution of stresses

[TABLE]

where

[TABLE]

Remark 2.4

In Definitions (17c) and (17d), the second and the third addends satisfy Equation (11), maintaining standard physical meaning of beam model variables.

Equations (16a) highlights that axial stress $\sigma_{x}$ explicitly depends on transversal internal force $V\left(x\right)$ and load $q$. Similarly, also the shear stress $\tau$ explicitly depends on the transversal load $q$. To the authors’ knowledge, only Karttunen and Von Hertzen [22] and Hashin [17] obtained similar dependencies, but their analysis was limited to homogeneous beams. Furthermore, $d_{\sigma_{x}}^{V}$ and $d_{\sigma_{x}}^{q}$ depend on $E_{xx}/G_{x}$ and $E_{xx}^{2}/G_{x}^{2}$, respectively (see Equation (17c) and (17d)), and similar coefficients was reported also by Karttunen and Von Hertzen [22].

On the one hand, similarity of Equations (16a) and (16b) with analytical solutions reported in [22] and [17] indicates that the procedure summarized in Figure 2 may be effective. On the other hand, non-trivial dependency of $\sigma_{x}$ on transversal internal force $V\left(x\right)$ indicates that stress recovery procedures developed for isotropic or orthotropic structural elements available in literature [13, 34] and implemented in most of structural analysis commercial softwares can lead to coarse results.

2.4 Beam constitutive relations

To complete the Timoshenko-like beam model, simplified constitutive relations have to be defined. To this aim, we introduce the stress potential

[TABLE]

Substituting the stress recovery relations (16) into Equation (18), the generalized strains result as the cross-section integral of the derivatives of the stress potential with respect to the corresponding internal forces, reading

[TABLE]

with

[TABLE]

Introducing the definitions of stress distributions (16a) into Definition (20b), we obtain that $\epsilon_{M}=\chi_{N}=0$ due to the choice of the origin of the Cartesian coordinate system introduced in Equation (8). Conversely, the transversal internal force $V\left(x\right)$ produces not only shear deformation $\gamma\left(x\right)$, but also axial strain $\epsilon\left(x\right)$ and curvature $\chi\left(x\right)$ since $\epsilon_{V}\neq 0$ and $\chi_{V}\neq 0$. Equations (20a) and (20d) lead to a definition of axial and bending stiffness analogous to the one obtained for isotropic beams. Conversely, Equation (20f) highlights that the shear stiffness of anisotropic beams $\gamma_{V}$ depends not only on shear modulus $G\left(y\right)$, but also on both the axial modulus of elasticity $E_{xx}\left(y\right)$ and the coupling term $G_{x}\left(y\right)$. Finally, all the deformations explicitly depend on the transversal load $q$ (see Equation (19)).

Such an deep influence of the material anisotropy on beam constitutive relations is ignored by most of the literature. To the authors’ knowledge, the coefficients $\epsilon_{V}=\gamma_{N}$ was analyzed only in [28, 29, 20, 31, 35]. Conversely, the existence of the coefficients $\chi_{V}=\gamma_{M}$ was mentioned by [23] in the framework of the derivation of an enhanced 3D beam model, but their influence on the beam structural response was never analyzed.

Remark 2.5

The extremely simple assumptions on kinematics (9) do not allow to tackle any higher order and boundary effects, as usual for all FSDT. Therefore, the proposed beam model has not the capability to describe deformation of the cross-section and the phenomena that occur in the neighborhood of constraints and concentrated loads.

3 ODEs analytical solution

This section discusses the analytical solution of beam model ODEs (4), (10), and (19) for the two-layer cantilever depicted in Figure 3.

The two layers are made of the same anisotropic material and their thicknesses are $h_{1}=\alpha h$ and $h_{2}=\left(1-\alpha\right)h$ with $0\leq\alpha\leq 1$. In the bottom layer, material principal direction is aligned with the beam axis, therefore $E_{xx}=E_{11}$, $G=G_{12}$, and $1/G_{x}=0$. In the top layer material principal direction is rotated with respect to the beam axis of an angle $\theta$, therefore $E_{xx}=E_{11}/\mu$ and $G=G_{12}/\kappa$. The dimensionless parameters $\mu$ and $\kappa$ account for the reduction of axial and shear modulus due to the rotation $\theta$ of the material principal direction (6). They are defined in Appendix A, together with the material coupling term $G_{x}$. For $\alpha=1$, the beam reduces to an homogeneous orthotropic beam with material orientation aligned with beam axis. Conversely, for $\alpha=0$, the beam reduces to an homogeneous anisotropic beams, similar to the one analyzed by [28, 29, 20, 31, 35, 17, 22]. Finally, aiming at the maximal simplicity of the analytical solution, we are going to neglect the influence of transversal load on the beam deformation (i.e., we assume $\epsilon^{q}=\chi^{q}=\gamma^{q}=0$).

The coefficients of beam constitutive relations introduced in Definition (19) read

[TABLE]

where the dimensionless coefficients $P_{i}$ for $i=1\dots 6$ (reported in Appendix B) account for the influence of geometry and material on the structural element stiffness.

The solution of ODEs (10), (4), and (19) leads to the following analytical expressions for beam model variables

[TABLE]

where $C_{i}$ for $i=1\dots 6$ depend on BCs. Notations $\left(\cdot\right)_{EB}$, $\left(\cdot\right)_{T}$, $\left(\cdot\right)_{c}$, and $\left(\cdot\right)_{r}$ highlight dependency of addends on the axial stiffness $E_{11}$, the shear stiffness $G_{12}$, the coupling term $G_{x}$, and the ratio $G_{x}^{2}/E_{11}$, respectively. Furthermore, few calculations allow to conclude that the addends denoted as $\left(\cdot\right)_{EB}$ coincide with the solution of the Euler Bernoulli (EB) beam theory whereas the addend denoted as $\left(\cdot\right)_{T}$ coincides with the shear deformation considered by Timoshenko beam theory.

Considering a cantilever (see Figure 3), the following BCs have to be enforced

[TABLE]

Requiring ODEs solution (22) to satisfy BCs (23) leads to determine the following value of $C_{i}$ for $i=1\dots 6$

[TABLE]

Finally, introducing the dimensionless parameter $\lambda=l/h$, the maximal transversal displacement of the beam reads

[TABLE]

where the dimensionless coefficients $Q_{i}$ for $i=1\dots 4$ are reported in Appendix C.

Equation (25) shows that the maximal transversal displacement $v\left(l\right)$ is the sum of four terms. The first addend $v_{EB}\left(l\right)$ depends on the Young’s modulus along the principal direction $E_{11}$ and, for $\alpha=1$, it corresponds to the classical EB solution. The second addend $v_{T}\left(l\right)$ depends on the shear modulus $G_{12}$ and, for $\alpha=1$, it coincides with the contribution due to shear deformation handled by the Timoshenko beam model. The third term $v_{c}\left(l\right)$ depends on the material coupling term $G_{x}$ and its existence is just a consequence of the fact that material principal directions are not aligned with the beam axis. The fourth term $v_{r}\left(l\right)$ depends on the material coupling term $G_{x}$ and on the ratio $E_{11}/G_{x}$ that appears in the definition of axial stress (see Equation (17c)).

Looking at Equation (25) from a different perspective, the first addend $v_{EB}\left(l\right)$ depends on $\lambda^{3}$, the second $v_{T}\left(l\right)$ and the fourth $v_{r}\left(l\right)$ ones depend on $\lambda$, and the third one $v_{c}\left(l\right)$ depends on $\lambda^{2}$. On the one hand, the so far highlighted result is conformal to what stated in standard literature. Shear deformation $v_{T}\left(l\right)$ has a negligible influence on the total displacement of the beam for slender beams (i.e., for $\lambda\gg 1$). On the other hand, the third term $v_{c}\left(l\right)$ can weigh on the total displacement more than the shear deformation whereas the forth term $v_{r}\left(l\right)$ can have an influence similar to the shear deformation. To the authors’ knowledge, the existence of terms $v_{c}\left(l\right)$ and $v_{r}\left(l\right)$ was never mentioned in the literature and their role will be analyzed in the following section.

Finally, it is worth mentioning that

stress distributions (16) reduce to linear and quadratic functions for $\alpha=1$ and $1/G_{x}=0$, as usual in homogeneous prismatic beams,

• 2.

$\epsilon_{N}=1/\left(E_{11}h\right)$, $\epsilon_{V}=\gamma_{N}=\chi_{V}=\gamma_{M}=0$, $\chi_{M}=12/\left(h^{3}E_{11}\right)$, and $\gamma_{V}=6/\left(5G_{12}h\right)$ for $\alpha=1$, analogously to homogeneous prismatic isotropic beams,

• 3.

$\epsilon_{N}=\mu/\left(E_{11}h\right)$, $\epsilon_{V}=\gamma_{N}=1/\left(G_{x}h\right)$, $\chi_{M}=12\mu/\left(h^{3}E_{11}\right)$, $\chi_{V}=\gamma_{M}=0$, and $\gamma_{V}=6\kappa/\left(5G_{12}h\right)$ for $\alpha=0$, similarly to anisotropic beam model proposed by Murakami et al. [28].

confirming that the presented beam model can recover analytical solutions already available in literature.

4 Comparison with analytical solution, simply-supported homogeneous beam,

This section compares the solution of the beam model discussed in Section 2 with the analytical solution derived in [22] for a simply supported homogeneous beam. Numerical results are obtained assuming the following parameters

[TABLE]

It is worth mentioning that the material is highly anisotropic: $E_{11}/E_{22}=20$ and $E_{11}/G=10$. Such a choice aims at magnifying the effects of both shear deformation and coupling on the behavior of the structural element, allowing for a more accurate discussion of the beam model effectiveness.

Being $\psi\left(x\right)$ a beam model variable, the solution computed by means of ODEs (10), (4), and (19) is denoted in the following as $\psi^{mod}$. Conversely, the reference solution $\psi^{ref}$ is computed using the analytical expressions reported in [22].

Assuming $\theta=45\,\deg$ the rotation of the constitutive relation (6) leads to set

[TABLE]

Figure 4 compares the cross-section distribution of stresses evaluated according to reference [22] $\psi^{ref}$, the proposed beam model $\psi^{mod}$, and standard Timoshenko beam $\psi^{T}$.

Numerical results demonstrate that the proposed beam model provides results substantially identical to the reference solution. In particular, Figure 4(b) highlights that shear does not vanish in the beam mid-span $\tau\left(l/2,y\right)\neq 0$, despite the vertical internal force vanishes $V\left(l/2\right)=0$. Further comments about this peculiarity of simply supported beams can be found in [22]. More interestingly, Figure 4(b) highlights that axial stress does not vanish at the bearing $\sigma_{x}\left(l,y\right)\neq 0$, despite both bending moment and axial internal force vanish $M\left(l\right)=N\left(l\right)=0$. Considering also Figures 4(c) and 4(d), it is possible to conclude that anisotropy influences the distribution of both axial and shear stresses. In particular, stress-recovery procedures developed for isotropic structural elements can underestimate the maximal magnitude of axial stress with errors greater than $10\,\%$.

5 Comparison with 2D FE, bi-layer beam

This section reports numerical results for two examples: Subsection 5.1 considers the cantilever already introduced in Section 3 and Subsection 5.2 analyzes a doubly-clamped beam (see Figure 3). In both cases, numerical results are obtained assuming the following parameters

[TABLE]

In this section, the reference solution $\psi^{ref}$ is computed using the commercial software Abaqus [10], in which the 2D problem domain $\Omega$ was discretized with a structured mesh of square bilinear elements CPS4. As discussed in Section 1, boundary effects could significantly affect the structural element behavior. Aiming at limiting their influence in reference solution, the BCs are imposed requiring only vanishing mean value of cross-section displacements and rotation. In this manner constrained cross-sections can warp and deform, but stress concentrations are limited in magnitude. Aiming at guaranteeing negligible numerical errors in the reference results, a sequence of analysis has been performed considering the bilayer cantilever and defining the element size $\delta$ according to the series $1/2^{n}$ for $n=0,1,2,\dots$. The procedure has been interrupted when the relative increase of the maximal displacement magnitude was smaller than $10^{-4}$, leading to set $\delta=0.25$

Transversal displacement $v^{ref}\left(x\right)$ and shear strain $\gamma^{ref}\left(x\right)$ have been obtained computing the mean value over the cross-section of the 2D transversal displacements $s_{y}^{ref}\left(x,y\right)$ and shear strains $\gamma_{xy}^{ref}\left(x,y\right)$, respectively. Conversely, axial $N^{ref}\left(x\right)$ and shear $V^{ref}\left(x\right)$ internal forces have been obtained as the integral over the cross-section of stress components $\sigma_{x}^{ref}\left(x,y\right)$ and $\tau^{ref}\left(x,y\right)$, respectively. The bending moment $M^{ref}\left(x\right)$ has been obtained as the integral over the cross section of axial stress $\sigma_{x}^{ref}\left(x,y\right)$ times the $y$ coordinate. Finally, according to Remark 2.3, the axial displacement $u^{ref}\left(x\right)$ and the rotation $\phi^{ref}\left(x\right)$ have been computed as the coefficients of the linear least squares with respect to $y$ of the axial displacements $s_{x}^{ref}\left(x,y\right)$. Similarly, the axial strain $\epsilon^{ref}\left(x\right)$ and the curvature $\chi^{ref}\left(x\right)$ have been computed as the coefficients of the linear least squares with respect to $y$ of the strains $\epsilon_{x}^{ref}\left(x,y\right)$.

5.1 Cantilever

In the following we set $l=500\,\mathrm{m}\mathrm{m}$ i.e., we choose $\lambda=5$. This assumption leads to consider a beam geometry that is close to the well known limit of validity of the FSDTs and, therefore, it will allow to identify every potential critical issue of the proposed model. Assuming $\theta=15\,\deg$ the rotation of constitutive relation (6) leads to set

[TABLE]

Figure 5 reports numerical results concerning the displacements.

Figures 5(b) and 5(d) analyze the model solution, highlighting the deep influence of the material coupling term $G_{x}$ on global structural response. In particular, Figure 5(d) highlights that transversal displacement component $v_{c}\left(x\right)$ (see Equation (22)) has a magnitude similar to the shear deformation component $v_{T}\left(x\right)$ and it contributes to total transversal displacement more than $10\,\%$. Conversely, $v_{r}\left(x\right)$ influences the total displacement less than $1\,\%$. Figure 5(b) highlights that $\phi_{c}\left(x\right)$ contributes to total cross-section rotation up to $10\,\%$. Figures 5(c) and 5(e) solution reveal a good accuracy of the proposed model. Indeed, relative errors are smaller than $2\,\%$ for rotation and transversal displacements.

Due to considered loads $N^{mod}\left(x\right)=0$ and, due to BCs, also $u_{EB}\left(x\right)=0$ (see Equation (22)). As a consequence, $u^{mod}\left(x\right)=u_{c}\left(x\right)$ i.e., the axial displacement is uniquely controlled by the material coupling term $G_{x}$ and the transversal internal force $V\left(x\right)$. Figure 5(a) shows that beam model correctly predict a non-vanishing distribution of axial displacement, but the error is near to $8\,\%$. Nevertheless, the axial displacement is two order of magnitude smaller than the transversal one, not affecting the errors evaluated on the total displacement $\boldsymbol{s}\left(x,y\right)$.

Figure 6 reports numerical results concerning generalized strains.

Figure 6(a) shows that axial strain is uniquely attributed to the transversal internal force $V\left(x\right)$ by means of the coefficient $\epsilon_{V}$ as already discussed above (see also Figure 5(a)). The comparison with reference solution (Figure 6(b)) confirms the goodness of the estimation provided by the beam model. Anyway, the reference solution reveals the presence of some higher order effects close to the clamp that can not be detected by the proposed beam model (see Remark 2.5) and may be also responsible of the errors on axial displacements.

Figure 6(c) shows that transversal internal force $V\left(x\right)$ produces non-negligible curvature, up to $10\,\%$ of the total. Similarly, Figure 6(e) shows that bending moment $M\left(x\right)$ deeply influences the shear strain which has a non-linear distribution despite transversal internal force $V\left(x\right)$ is linear. In particular, bending moment $M\left(x\right)$ produces non-negligible shear strain, up to $30\,\%$ of the total. For both shear deformation and curvature, Figures 6(d) and 6(f) demonstrate that generalized strains predicted by the beam model are in extremely good agreement with reference solution. Only near to the clamp, reference solution reveals the presence of some higher order effects that are not handled by the beam model.

Figures 7 and 8 report cross-section stress distributions at $1/2l=250\,\mathrm{m}\mathrm{m}$ and $3/4l=375\,\mathrm{m}\mathrm{m}$.

Figures 7(a) and 7(c) highlight that the axial stress depending on the transversal internal force $d_{\sigma_{x}}^{V}\left(y\right)V\left(x\right)$ is not negligible at all, but can increase the magnitude of maximal stress up to $30\,\%$. Conversely, the effect of the axial stress depending on the transversal load $d_{\sigma_{x}}^{q}\left(y\right)q$ is less significant. The comparison with reference solution (Figures 7(b) and 7(d)) reveals that the stress recovery developed in Section 2.3 provides accurate estimations of the stress magnitude, with relative errors rarely bigger than $10\,\%$. In particular, the stress recovery correctly predicts the jump of axial stress at the interlayer surface. Conversely, reference solution reveals the presence of some higher order effects near to the free end of the cantilever that the beam model is not able to catch (see Remark 2.5), and may locally lead to an increase of the relative errors up to $40\,\%$.

Figures 8(a) and 8(c) highlight that the shear stress depending on transversal load $d_{\tau}^{q}\left(y\right)q$ is not negligible, but it can lead to the creation of a local minimum on the interlayer surface. The comparison with reference solution (Figures 8(a) and 8(c)) reveals that the proposed stress recovery is in good agreement with reference solution, leading to relative errors rarely bigger than $5\,\%$.

In order to complete the discussion of the proposed model capabilities, Tables 1, 2, and 3 compare the solutions obtained using different models. Maximal displacements $\psi\left(l\right)$, with $\psi=u,\phi,v$, are evaluated using 2D FE $\psi^{ref}\left(l\right)$, standard EB beam model $\psi^{EB}\left(l\right)$, and the proposed beam model $\psi^{mod}\left(l\right)$. Only for the transversal displacement, also the Timoshenko beam model is considered since its solution differs form EB ($v^{EB}\left(l\right)+v^{T}\left(l\right)$. Relative errors are computed as

[TABLE]

Finally, numerical results and relative errors are provided for $\lambda=5,\,10,\mbox{ and }20$ and $\theta=\pm 15\,\deg$.

On the one hand, relative errors decrease increasing the slenderness for all the considered beam models, consistently with standard beam model assumptions. On the other hand, it is worth highlighting that EB and Timoshenko beam models lead to errors that are often greater than $10\,\%$ and are therefore not acceptable for most of engineering applications. Conversely, the proposed beam model leads to errors that are usually smaller than $5\,\%$, reaching an accuracy adequate for most of engineering applications.

Furthermore, Table 3 highlights that Timoshenko beam model not always performs better than EB, even considering tick beams. In particular, for $\lambda=5$ and $\theta=-15\,\deg$ maximal transversal displacement predicted by the proposed beam model qualitatively coincides with the EB solution i.e., $v^{mod}\left(l\right)\approx v^{EB}\left(l\right)$. On the one hand, such a result highlights (i) the deep influence of fiber direction on the structural element stiffness (see Remarck 2.1) and (ii) the inappropriateness of beam models developed for isotropic structural element. On the other hand, the extremely low errors obtained for both positive and negative $\theta$ confirm the effectiveness of the proposed model in handling all peculiar aspects of anisotropic beams.

Figure 9 reports the weight of the four components $v_{EB}\left(l\right)$, $v_{T}\left(l\right)$, $v_{c}\left(l\right)$, and $v_{r}\left(l\right)$ on the total transversal displacement as a function of $\lambda$.

The analysis is limited to geometry and material mechanical properties introduced at the beginning of Section 5, but it highlights several effects of the anisotropy on the structural response of beams. The component $v_{c}\left(l\right)$, depending on the material coupling term $G_{x}$, is always bigger than component depending of shear $v_{T}\left(l\right)$. As an example, considering $\lambda=20$ $v_{c}\left(l\right)/v\left(l\right)>3\,\%$ whereas $v_{T}\left(l\right)/v\left(l\right)<1\%$. Furthermore, for $\lambda=80$ $v_{c}\left(l\right)/v\left(l\right)\approx 1\,\%$ whereas $v_{T}\left(l\right)/v\left(l\right)\approx 0.05\,\%$. As a consequence, it is possible to conclude that the material coupling contribution $v_{c}\left(x\right)$ could be significant also for slender structural elements for which, instead, shear contribution is negligible. Conversely, the transversal displacement $v_{r}\left(l\right)$ always contributes to the total displacement less than $1\,\%$.

5.2 Doubly-clamped beam

This section considers the statically indeterminate multilayer beam depicted in Figure 3, aiming at confirming the capabilities of the proposed beam model in effectively estimating anisotropic beam stiffness. Analytical solution reported in Equation (22) is still valid. Conversely the following BCs has to be considered:

[TABLE]

Analytical expression for the coefficients $C_{i}$ for $i=1\dots 6$ turns out to be extremely complex and, for brevity, they will not be reported. Anyway, it has to be noticed that all the coefficients depends on all mechanical properties $E_{11}$, $G_{12}$ , and $G_{x}$. As a consequence, the subdivision of displacements in components $\left(\cdot\right)_{EB}$, $\left(\cdot\right)_{T}$, $\left(\cdot\right)_{c}$, and $\left(\cdot\right)_{r}$ introduced in Equation (22) is no longer meaningful and it will not be considered in the following. We set $l=1000\,\mathrm{m}\mathrm{m}$ and we use the geometrical and mechanical properties reported in Equations (28) and (29).

Figure 10 reports numerical results concerning distribution of internal forces $N\left(x\right)$, $M\left(x\right)$, and $V\left(x\right)$.

Numerical results highlight a non-trivial effect of the material anisotropy. The distribution of internal forces is non-symmetric and reactions on the right hand side clamp are greater than the ones in the left hand side clamp despite BCs and load are symmetric with respect to the beam mid-span. Once more, comparison with reference solution reveals the high accuracy of the proposed beam model that predicts both bending moment and transversal internal force with negligible errors. Finally, the proposed beam model correctly predicts a non-vanishing, constant distribution of axial internal force, which magnitude is anyway negligible if compared with transversal internal force and bending moment.

Tables 4, 5, and 6 report constraint reactions (i.e., $N\left(x\right)$, $M\left(x\right)$, and $V\left(x\right)$ for $x=0,l$) evaluated using 2D FE $\psi^{ref}\left(l\right)$, standard EB beam model $\psi^{EB}\left(l\right)$, and the proposed beam model $\psi^{mod}\left(l\right)$ for $\lambda=5,\,10,\mbox{ and }20$. Relative errors are computed according to Equation (30).

As already remarked in Section 5.1, EB beam model leads to errors that are often bigger than $10\,\%$, leading to estimations that are too coarse for most of engineering applications. Conversely, the proposed beam model leads to errors that are generally three-six times smaller and always below $5\,\%$. Only Table 4 highlights that model estimates axial internal force $N\left(x\right)$ with errors over $20\,\%$. Anyway, since the magnitude of axial internal force is approximatively 50 times smaller than transversal one, the relative error on axial internal force might not have a deep influence on the global response of structural element.

6 Conclusions

This paper has proposed a simple beam model that effectively handles the influence of anisotropy on the beam constitutive relations and the stress distribution. The independent variables of the model are the internal forces and the standard Timoshenko kinematic parameters. Despite its simplicity, the beam model has allowed to highlight the following peculiarities of anisotropic beams.

Material anisotropy leads transversal internal force to contribute up to $30\,\%$ of the magnitude of axial stress, deeply affecting also the beam strength, not explicitly considered in this paper. 2. 2.

In beam constitutive relations, non-vanishing out-of-diagonal terms that relate transversal internal force with curvature (and bending moment with shear strain) exist and deeply influence the response of the structural element. 3. 3.

In addition to the standard bending contribution (proportional to cube beam-slenderness) and the shear one (proportional to beam-slenderness), a third term, depending on material coupling term and proportional to square beam-slenderness, contributes to transversal displacement. 4. 4.

The contribution depending on material coupling terms can be bigger than the contribution given by shear deformation and it may be non-negligible for length vs thickness ratios greater than fifty.

A systematic comparison with analytical results and 2D FE solutions, obtained using highly refined meshes, demonstrates the effectiveness of the proposed modeling approach. In general, the proposed beam model has a computational cost similar to simplest beam models used in engineering practice and it estimates significant displacements and internal forces with relative errors usually smaller than $5\,\%$. Conversely, coarse adaptations of beam models developed for isotropic structural elements may lead to errors greater than $20\,\%$ in the prediction of both internal forces and displacements. Furthermore, analysis of stress distributions demonstrates that stress recovery tools developed for isotropic structural elements are no longer effective for anisotropic ones, but ad-hoc routines has to be developed. The main limitations of the proposed model are the assumptions on kinematics that do not allow to describe higher order effects like cross-section warping and distortion as well as phenomena that occur in the neighborhood of constraints and concentrated loads.

Future research will include the application of the proposed modeling strategy to higher order planar beams and its generalization to 3D beams and plates.

7 Acknowledgments

This work was funded by the Austrian Science Found (FWF) [M 2009-N32]. F. Auricchio and S. Morganti would like to acknowledge the strategic theme of the University of Pavia ”Virtual Modeling and Additive Manufacturing for Advanced Materials (3D@UniPV)”

Appendix A Mechanical properties coefficients

[TABLE]

[TABLE]

[TABLE]

Appendix B Dimensionless coefficients $P_{i}$ (for $i=1\dots 6$)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Maximal displacement coefficients $Q_{i}$ (for $i=1\dots 3$)

[TABLE]

[TABLE]

[TABLE]

[TABLE]