You do not have permission to edit this page, for the following reason:

You are not allowed to execute the action you have requested.


You can view and copy the source of this page.

x
 
1
<!-- metadata commented in wiki content
2
==SIMPLE AND ACCURATE TWO-NODED BEAM ELEMENT FOR COMPOSITE LAMINATED BEAMS USING A REFINED ZIGZAG THEORY==
3
4
'''E. Oñate<math>^{1,2}</math>, A. Eijo<math>^{1}</math> and S. Oller<math>^{1,2}</math>'''
5
6
{|
7
|-
8
|<math>^{1}</math> International Center for Numerical Methods in Engineering (CIMNE) 
9
|-
10
| <math>^{2}</math> Universitat Politecnica de Catalunya (UPC)
11
|-
12
| Campus Norte UPC, 08034 Barcelona, Spain 
13
|-
14
| [mailto:onate@cimne.upc.edu onate@cimne.upc.edu], [http://www.cimne.com/eo www.cimne.com/eo]
15
|}
16
-->
17
18
==Abstract==
19
20
In this work we present a new simple linear two-noded beam element adequate for the analysis of composite laminated and sandwich beams based on the combination of classical Timoshenko beam theory and the refined zigzag kinematics  proposed by Tessler ''et al.'' <span id='citeF-22'></span>[[#cite-22|[22]]]. The element has just four kinematic variables per node. Shear locking is eliminated by reduced integration. The accuracy of the new beam element is tested in a number  of applications to the analysis of composite laminated beams with simple supported and clamped ends under point loads and uniformly distributed loads. An example showing the capability of the new element for accurately reproducing delamination effects is also presented.
21
22
'''keywords''' Two-noded beam element, zigzag kinematics, Timoshenko theory, composite, sandwich beams
23
24
==1 INTRODUCTION==
25
26
It is well known that both the classical Euler-Bernouilli beam theory <span id='citeF-1'></span>[[#cite-1|[1]]] and the more advanced Timoshenko theory  <span id='citeF-2'></span>[[#cite-2|[2]]] produce inadequate predictions when applied to relatively thick composite laminated  beams with material layers that have highly dissimilar stiffness characteristics. Even with a judiciously chosen shear correction factor, Timoshenko theory tends to underestimate the axial stress at the top and bottom outer fibers of a beam. Also, along the layer interfaces of a laminated beam the transverse shear stresses predicted exhibit erroneous discontinuities. These difficulties are due to the higher complexity of the “true” variation of the axial displacement field across a highly heterogeneous beam cross-section.
27
28
Indeed to achieve accurate computational results, 3D finite element analyses are often preferred over beam models. For composite laminates with hundred of layers, however, 3D modelling becomes prohibitely expensive, specially for non linear and progressive failure analyses.
29
30
Improvements to the classical beam theories have been obtained by the so called equivalent single layer (ESL) theories that assume a priori the behaviour of the displacement and/or the stress through the laminate thickness <span id='citeF-3'></span>[[#cite-3|[3,4]]]. Despite of being computational efficient, ESL theories often produce innacurate distributions for the stresses and strains (in particular the transverse shear stress) across the thickness.
31
32
The need for composite laminated beam theories with better predictive capabilities has led to the development of the so-called ''higher order'' theories. In  these theories higher-order kinematic terms with respect to the beam depth  are added to the expression for the axial displacement and, in some cases, to the expressions for the deflection. A review of these theories can be found in <span id='citeF-3'></span>[[#cite-3|[3,4]]].
33
34
Accurate predictions of the correct shear and axial stresses for thick and highly heterogenous composite laminated and sandwich beams can be obtained by using ''layer-wise'' theory. In this theory the thickness coordinate is split into a number of ''analysis layers'' that may or not coincide with the number of laminate plies. The kinematics are independently described within each layer and certain physical continuity requirements are enforced <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[3,4]]].
35
36
A drawback of layer-wise theory is that the number of kinematic variables depends on the number of analysis layers. The layer displacements  can be condensed at each section in terms of the axial displacement for the top layer during the equation solution process <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[5,6]]]. The displacement condensation processes can be however expensive for problems involving many analysis layers.
37
38
Discrete layer theories in which the number of unknowns in the model does not depend on the number of layers in the laminate are described in <span id='citeF-7'></span>[[#cite-7|[7,8,9]]]. In this class of discrete layerwise theories (called zigzag theories) a piewise linear in-plane displacement function (the zigzag function) is superimposed over a linear displacement field <span id='citeF-7'></span>[[#cite-7|[7,8]]], a quadratic displacement field <span id='citeF-10'></span>[[#cite-7|[10,11]]] or a cubic displacement field <span id='citeF-12'></span>[[#cite-12|[12,13,14]]] through the thickness of the laminate.
39
40
Many zigzag theories  require <math display="inline">C^1</math> continuity for the deflection field, which is a drawback versus simpler <math display="inline">C^\circ </math> continuous FEM approximations. Also many zigzag theories run into theoretical difficulties to satisfy equilibrium of forces at a clamped support.
41
42
Averill et al. <span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-16'></span><span id='citeF-17'></span>[[#cite-9|[9,10,16,17]]] developed linear, quadratic and cubic zigzag beam theories that overcame the need for <math display="inline">C^1</math> continuity. The shear strain angle is introduced as a kinematic variable together with the deflection, the rotation and the zigzag function. A <math display="inline">C^\circ </math> interpolation can be used for all these variables. The relationship between the shear angle, the deflection and the rotation of each layer is introduced as a constraint via a penalty method. This also ensures the continuity of the transverse shear stress across the laminate depth and the satisfaction of the shear traction boundary conditions. However, Averill theories have difficulties to model correctly clamped boundary conditions. For this reason, analytical and numerical (FEM) studies based on Averill theory have mainly focused on simple supported beams <span id='citeF-16'></span><span id='citeF-17'></span>[[#cite-16|[16,17]]].
43
44
A 2-noded beam element based on Euler-Bernouilli beam theory and an extension of Averill's zigzag theory including a cubic in-plane displacement field within each layer has been recently proposed by Alam and Upadhyay <span id='citeF-18'></span>[[#cite-18|[18]]]. Good results are reported for cantilever and clamped composite and sandwich beams.
45
46
An assessment of different zigzag theories for beam is reported in <span id='citeF-19'></span><span id='citeF-20'></span>[[#cite-19|[19,20]]]. A review of zigzag theory for plate analysis can be found in <span id='citeF-21'></span>[[#cite-21|[21]]].
47
48
Tessler ''et al.'' <span id='citeF-22'></span>[[#cite-22|[22,23]]] have developed a refined zigzag (RZ) theory  starting from the standard Timoshenko kinematic assumptions. This allows one using <math display="inline">C^\circ </math> continuous interpolation for all the kinematic variables. Timoshenko beam theory also introduces naturally shear deformation effect for the homogeneous material case, which is advantageous for many problems. The zigzag functions chosen have the property of vanishing on the top and bottom surfaces of a laminate. A particular feature of this zigzag theory is that the transverse shear stresses are not required to be continuous at the layer interfaces. This results in simple piewice-constant functions that approximate the true shear stress distribution. An accurate continuous thickness distribution of the transverse shear stress can be obtained “a posteriori” in terms of the axial stress by integrating the equilibrium equations. This theory also provides good results for clamped supports.
49
50
Gherlone et al. <span id='citeF-24'></span>[[#cite-24|[24]]] have  developed two and three-noded <math display="inline">C^\circ </math> beam elements based on the RZ theory for analysis of multilayered composite and sandwich beams. Locking-free elements are obtained by using  ''anisoparametric'' interpolations that are adapted to approximate the four independent kinematic variables  that model the beam deformation. A family of beam elements is achieved by imposing different constraints  on the original displacement approximation. The constraint conditions requiring a constant variation of the transverse shear force provide an accurate 2-noded beam element <span id='citeF-24'></span>[[#cite-24|[24]]].
51
52
Quite simultaneously to the above work, Oñate et al. <span id='citeF-25'></span>[[#cite-25|[25]]] proposed a simple 2-noded beam element for composite laminated beams based on the RZ theory. A standard linear displacement field is used to model the four variables of the so called LRZ element. Shear locking is avoided by using reduced integration on selected terms of the shear stiffness matrix.
53
54
In this paper we present in detail the formulation of the LRZ beam element originally reported in  <span id='citeF-25'></span>[[#cite-25|[25]]] and explore the capabilities of the new element for multilayered beams and delamination analysis. A study of the locking-free behaviour of the LRZ element for slender beams is presented. The good performance of the element is demonstrated for simply supported and clamped composite laminated beams with different layers under point load and uniformly distributed loads.
55
56
Finally, an example showing the capability of the LRZ element to model delamination effects is presented.
57
58
==2 GENERAL CONCEPTS OF ZIGZAG BEAM THEORY==
59
60
The kinematic field in zigzag beam theory is generally written as
61
62
{| class="formulaSCP" style="width: 100%; text-align: left;" 
63
|-
64
| 
65
{| style="text-align: left; margin:auto;width: 100%;" 
66
|-
67
| style="text-align: center;" | <math>u^k (x,z)  = u_0 (x) - z\theta (x) + \bar u^k (x,z)</math>
68
|-
69
| style="text-align: center;" | <math> w(x,z)  =w_0(x) </math>
70
|}
71
| style="width: 5px;text-align: right;white-space: nowrap;" | (1.a)
72
|}
73
74
where
75
76
{| class="formulaSCP" style="width: 100%; text-align: left;" 
77
|-
78
| 
79
{| style="text-align: left; margin:auto;width: 100%;" 
80
|-
81
| style="text-align: center;" | <math>\bar u^k = \phi ^k (z)\Psi (x) \quad ,\quad k = 1, N </math>
82
|}
83
| style="width: 5px;text-align: right;white-space: nowrap;" | (1.b)
84
|}
85
86
is the ''zigzag displacement function'' (Figure [[#img-1|1]]).
87
88
In Eqs.(1) <math display="inline">N</math> is the number of layers, superscript <math display="inline">k</math> indicates quantities within the <math display="inline">k</math>th layer with <math display="inline">z_k \le z\le z_{k+1}</math> and <math display="inline">z_k</math> is the vertical coordinate of the <math display="inline">k</math>th interface. In Eq.(1a) the uniform axial displacement <math display="inline">u_0(x)</math>, the rotation <math display="inline">\theta (x)</math> and the transverse deflection <math display="inline">w_0(x)</math> are the primary kinematic  variables of the underlying equivalent single-layer Timoshenko beam theory. In Eq.(1b) function <math display="inline">\phi ^k (z)</math> denotes ''a piecewise linear zigzag function'', yet to be established, and <math display="inline">\Psi (x)</math> is a primary kinematic variable that defines the ''amplitude of the zigzag function'' along the beam. Collectively, the interfacial axial displacement field has a zigzag distribution, as shown in Figure [[#img-1|1]]c.
89
90
<div id='img-1'></div>
91
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
92
|-
93
|[[Image:Draft_Samper_624436055-test-Fig1.png|600px|Thickness distribution of the the zigzag function ϕ<sup>k</sup> (a), the zigzag displacement function ̄u<sup>k</sup> (b),  and the axial displacement (c) in refined zigzag beam theory]]
94
|- style="text-align: center; font-size: 75%;"
95
| colspan="1" | '''Figure 1:''' Thickness distribution of the the zigzag function <math>\phi ^k</math> (a), the zigzag displacement function <math>\bar u^k</math> (b),  and the axial displacement (c) in refined zigzag beam theory
96
|}
97
98
The  strain-displacement relations are derived by substituting Eq.(1a) into the expressions of classical Timoshenko beam theory, i.e.
99
100
{| class="formulaSCP" style="width: 100%; text-align: left;" 
101
|-
102
| 
103
{| style="text-align: left; margin:auto;width: 100%;" 
104
|-
105
| style="text-align: center;" | <math>\varepsilon _x^k = \frac{\partial u^k}{\partial x} = \frac{\partial u_0}{\partial x} - z \frac{\partial \theta }{\partial x}+ \phi ^k \frac{\partial \Psi }{\partial x} = [1,-z,\phi ^k]\left\{\begin{matrix}\displaystyle \frac{\partial u_0}{\partial x}\\[.3cm]\displaystyle \frac{\partial \theta }{\partial x} \\[.3cm]\displaystyle \frac{\partial \Psi }{\partial x} \end{matrix} \right\}= {\boldsymbol S}_p \hat {\boldsymbol \varepsilon }_p </math>
106
|}
107
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.a)
108
|}
109
110
{| class="formulaSCP" style="width: 100%; text-align: left;" 
111
|-
112
| 
113
{| style="text-align: left; margin:auto;width: 100%;" 
114
|-
115
| style="text-align: center;" | <math>\gamma _{xz}^k = \frac{\partial u^k}{\partial z}+ \frac{\partial w}{\partial x} = \frac{\partial w_0}{\partial x}-\theta + \frac{\partial \phi ^k}{\partial z}\Psi = \gamma +\beta ^k \Psi = [1,\beta ^k] \left\{\begin{matrix}\gamma \\ \Psi  \end{matrix} \right\}= {\boldsymbol S}_t \hat {\boldsymbol \varepsilon }_t </math>
116
|}
117
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.b)
118
|}
119
120
In Eq.s (2a) and (2b)
121
122
{| class="formulaSCP" style="width: 100%; text-align: left;" 
123
|-
124
| 
125
{| style="text-align: left; margin:auto;width: 100%;" 
126
|-
127
| style="text-align: center;" | <math>\begin{array}{ll}{\boldsymbol S}_p =[1,-z,\phi ^k]\quad ,  &  \hat {\boldsymbol \varepsilon }_p = \displaystyle \left[\frac{\partial u_0}{\partial x} ,  \frac{\partial \theta }{\partial x},\frac{\partial \Psi }{\partial x}\right]^T\\ {\boldsymbol S}_t =[1,\beta ^k]\quad ,  & \hat{\boldsymbol \varepsilon }_t = \left[\gamma ,\Psi \right]^T \end{array} </math>
128
|}
129
| style="width: 5px;text-align: right;white-space: nowrap;" | (2.c)
130
|}
131
132
where <math display="inline">\hat {\boldsymbol \varepsilon }_p</math> and <math display="inline">\hat{\boldsymbol \varepsilon }_t</math> are the generalized in-plane and transverse shear strain vectors, respectively. Vector <math display="inline">\hat {\boldsymbol \varepsilon }_p</math> contains the axial elongation <math display="inline">\left(\frac{\partial u_0}{\partial x}\right)</math>, the pseudo-curvature <math display="inline">\left(\frac{\partial \theta }{\partial x}  \right)</math> and the derivatives of the amplitude of the zigzag function <math display="inline">\left(\frac{\partial \Psi }{\partial x} \right)</math>. In <math display="inline">\hat{\boldsymbol \varepsilon }_t</math>, <math display="inline">\gamma = \frac{\partial w_0}{\partial x} -\theta </math> is the average transverse shear strain of Timoshenko beam theory and <math display="inline">\beta ^k = \frac{\partial \phi ^k}{\partial z}</math>. Note that since <math display="inline">\phi ^k (z)</math> is piecewise linear, <math display="inline">\beta ^k </math> is constant across each layer.
133
134
For major principal material axes that are coincident with the beam <math display="inline">x</math> axis, Hooke stress-strain relations for the <math display="inline">k</math>th orthotropic layer have the standard form
135
136
{| class="formulaSCP" style="width: 100%; text-align: left;" 
137
|-
138
| 
139
{| style="text-align: left; margin:auto;width: 100%;" 
140
|-
141
| style="text-align: center;" | <math>\sigma _x^k = E^k \varepsilon _x^k = E^k  {\boldsymbol S}_p\hat {\boldsymbol \varepsilon }_p</math>
142
| style="width: 5px;text-align: right;white-space: nowrap;" | (3.a)
143
|-
144
| style="text-align: center;" | <math> \tau _{xz}^k  = G^k \gamma _{xz}^k = G^k {\boldsymbol S}_t \hat{\boldsymbol \varepsilon }_t </math>
145
| style="width: 5px;text-align: right;white-space: nowrap;" | (3.b)
146
|}
147
|}
148
149
where <math display="inline">E^k</math> and <math display="inline">G^k</math> are the axial and shear moduli for the <math display="inline">k</math>th layer, respectively.
150
151
In the above equations we have distinguished all variables within a layer with superscript <math display="inline">k</math>.
152
153
==3 REFINED ZIGZAG THEORY==
154
155
===3.1 Zigzag kinematics===
156
157
The key attributes of the refined zigzag (RZ) theory  proposed by Tessler ''et al.'' [22] are, first, ''the zigzag function'' ''vanishes at the top and bottom surfaces'' ''of the beam'' section and does not require full shear-stress continuity across the laminated-beam depth. Second, all boundary conditions can be modelled adequately. And third, <math display="inline">C^\circ </math> continuity is only required for the FEM approximation of the kinematic variables.
158
159
Within each layer the zigzag function is expressed as
160
161
{| class="formulaSCP" style="width: 100%; text-align: left;" 
162
|-
163
| 
164
{| style="text-align: left; margin:auto;width: 100%;" 
165
|-
166
| style="text-align: center;" | <math>\phi ^k = \frac{1}{2} (1-\zeta ^k ) \bar \phi ^{k-1} + \frac{1}{2} (1+\zeta ^k ) \bar \phi ^k = \frac{\bar \phi ^k + \bar \phi ^{k-1}}{2}+ \frac{\bar \phi ^k - \bar \phi ^{k-1}}{2} \zeta ^k </math>
167
|}
168
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
169
|}
170
171
where <math display="inline">\bar \phi ^k</math> and <math display="inline">\bar \phi ^{k-1}</math> are the zigzag function values of the <math display="inline">k</math> and <math display="inline">k-1</math> interface, respectively with <math display="inline">\bar \phi ^0 = \bar \phi ^N =0</math> and <math display="inline">\zeta ^k = \frac{2(z - z^{k-1})}{h^k}-1</math>.
172
173
Collectively, function <math display="inline">\phi ^k</math> has the zigzag distribution shown in Figure [[#img-1|1]]a. Due to the dependence between the zigzag displacement function <math display="inline">\bar u^k</math> and <math display="inline">\bar \phi ^k</math> (see Eq.(1b)), <math display="inline">\bar u^k</math> also vanishes at the top and bottom layers. The axial displacement field is plotted in Figure [[#img-1|1]]c.
174
175
The above form of <math display="inline">\phi ^k</math> gives the constant value of <math display="inline">\beta ^k</math> for each layer as
176
177
{| class="formulaSCP" style="width: 100%; text-align: left;" 
178
|-
179
| 
180
{| style="text-align: left; margin:auto;width: 100%;" 
181
|-
182
| style="text-align: center;" | <math>\beta ^k = \frac{\partial \phi ^k}{\partial z} = \frac{\bar \phi ^k - \bar \phi ^{k-1}}{h^k} </math>
183
|}
184
| style="width: 5px;text-align: right;white-space: nowrap;" | (5.a)
185
|}
186
187
and
188
189
{| class="formulaSCP" style="width: 100%; text-align: left;" 
190
|-
191
| 
192
{| style="text-align: left; margin:auto;width: 100%;" 
193
|-
194
| style="text-align: center;" | <math>\iint _A \beta ^k dA =0 </math>
195
|}
196
| style="width: 5px;text-align: right;white-space: nowrap;" | (5.b)
197
|}
198
199
The <math display="inline">\beta ^k</math> parameter is useful for computing the zigzag function as explained in the next section.
200
201
===3.2 Computation of the zigzag function===
202
203
Integrating Eq.(2b) over the cross section and using Eq.(5b) and the fact that <math display="inline">\Psi </math> is independent of <math display="inline">z</math> yields
204
205
{| class="formulaSCP" style="width: 100%; text-align: left;" 
206
|-
207
| 
208
{| style="text-align: left; margin:auto;width: 100%;" 
209
|-
210
| style="text-align: center;" | <math>\gamma = \frac{1}{A} \iint _A \gamma _{xz}^k dA </math>
211
|}
212
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
213
|}
214
215
i.e. <math display="inline">\gamma </math> represents the average shear strain of the cross section, as expected.
216
217
The shear strain-shear stress relationship of Eq.(3b) is written as
218
219
{| class="formulaSCP" style="width: 100%; text-align: left;" 
220
|-
221
| 
222
{| style="text-align: left; margin:auto;width: 100%;" 
223
|-
224
| style="text-align: center;" | <math>\tau _{xz}^k =G^k \eta + G^k (1+\beta ^k) \Psi  </math>
225
|}
226
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
227
|}
228
229
where <math display="inline">\eta =\gamma - \Psi </math> is a difference function.
230
231
'''Remark 1'''. Function <math display="inline">\Psi </math>  can be interpreted as a weighted-average shear strain angle <span id='citeF-22'></span>[[#cite-22|[22]]]. The value of <math display="inline">\Psi </math> should be prescribed to zero at a clamped edge and left unprescribed  at free and simply supported edges.
232
233
Eq.(7) shows that the distribution of <math display="inline">\tau _{xz}^k</math> within each layer is constant, as <math display="inline">\eta </math> is independent of the zigzag function and <math display="inline">\beta ^k </math> is constant.
234
235
The distribution of <math display="inline">\tau _{xz}^k</math> is now ''enforced to be independent of the zigzag function''. This can be achieved by constraining the term multiplying <math display="inline">\Psi </math> in Eq.(7) to be constant, i.e.
236
237
{| class="formulaSCP" style="width: 100%; text-align: left;" 
238
|-
239
| 
240
{| style="text-align: left; margin:auto;width: 100%;" 
241
|-
242
| style="text-align: center;" | <math>G^k (1+\beta ^k) = G^{k+1} (1+\beta ^{k+1})=G, \quad \hbox{constant} </math>
243
|}
244
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
245
|}
246
247
This is equivalent to enforcing the interfacial continuity of the second term in the r.h.s. of Eq.(7).
248
249
'''Remark 2.''' We emphasize that this zigzag theory  does not enforce the continuity of the transverse shear stresses across the section. This is consistent with the kinematic freedom inherent in the lower order kinematic approximation of the underlying beam theory. An accurate continuous distribution of the transverse shear stress across the thickness of the laminate can be obtained “a posteriori” in terms of axial stresses by integrating the equilibrium equations as explained in Section 7.3.
250
251
From Eq.(8) we deduce
252
253
{| class="formulaSCP" style="width: 100%; text-align: left;" 
254
|-
255
| 
256
{| style="text-align: left; margin:auto;width: 100%;" 
257
|-
258
| style="text-align: center;" | <math>\beta ^k = \frac{G}{G^k} -1 </math>
259
|}
260
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
261
|}
262
263
Substituting <math display="inline">\beta ^k</math> in the integral of Eq.(5b) gives
264
265
{| class="formulaSCP" style="width: 100%; text-align: left;" 
266
|-
267
| 
268
{| style="text-align: left; margin:auto;width: 100%;" 
269
|-
270
| style="text-align: center;" | <math>G = \left[\frac{1}{A}\iint _A \frac{dA}{G^k} \right]^{-1}  =  \left[h \sum \limits _{k=1}^N \frac{h^k}{G^k}   \right]^{-1} </math>
271
|}
272
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
273
|}
274
275
where <math display="inline">h</math> is the section depth. Substituting Eq.(9) into Eq.(5a) gives the following recursion relation for the zigzag displacement function values at the layer interfaces
276
277
{| class="formulaSCP" style="width: 100%; text-align: left;" 
278
|-
279
| 
280
{| style="text-align: left; margin:auto;width: 100%;" 
281
|-
282
| style="text-align: center;" | <math>\bar u_k = \sum \limits _{i=1}^k h^i \beta ^i \quad \hbox{with}\quad u^0 =u^N=0 </math>
283
|}
284
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
285
|}
286
287
Introducing Eq.(11) into (4) gives the expression for the zigzag function as
288
289
{| class="formulaSCP" style="width: 100%; text-align: left;" 
290
|-
291
| 
292
{| style="text-align: left; margin:auto;width: 100%;" 
293
|-
294
| style="text-align: center;" | <math>\phi ^k = \frac{h^k\beta ^k}{2} (\zeta ^k -1)+  \sum \limits _{i=1}^k  h^i \beta ^i </math>
295
|}
296
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
297
|}
298
299
Recall that superindex <math display="inline">k</math> denotes the number of each material layer.
300
301
'''Remark 3.''' For homogeneous material <math display="inline">G^k= G</math> and <math display="inline">\beta ^k =0</math>. Hence, the zigzag function <math display="inline">\phi ^k</math> vanishes and we recover the kinematics and constitutive expressions of the standard Timoshenko composite laminated beam theory. This is a particular feature of this zigzag theory.
302
303
'''Remark 4.''' Note that differently from standard Timoshenko beam theory, a shear correction parameter is not needed in the RZ theory.
304
305
306
===3.3 Constitutive relationship===
307
308
The in-plane bending and transverse shear resultant stresses are defined as
309
310
{| class="formulaSCP" style="width: 100%; text-align: left;" 
311
|-
312
| 
313
{| style="text-align: left; margin:auto;width: 100%;" 
314
|-
315
| style="text-align: center;" | <math>\hat{\boldsymbol \sigma }_p = \left\{\begin{array}{c}N \\                             M\\ M_\phi                            \end{array} \right\}=\iint _A {\boldsymbol S}_p^T \sigma _x^k  dA = \left(\iint _A {\boldsymbol S}_p^T {\boldsymbol S}_p E^k dA\right)\hat{\boldsymbol \varepsilon }_p =\hat {\boldsymbol D}_p \hat{\boldsymbol \varepsilon }_p</math>
316
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
317
|-
318
| style="text-align: center;" | <math> \hat{\boldsymbol \sigma }_t = \left\{\begin{array}{c}Q \\ Q_\phi                            \end{array}  \right\}= \iint _A {\boldsymbol S}_t^T \tau _{xz}^k dA = \left(\iint _A  {\boldsymbol S}_t^T {\boldsymbol S}_t G^k dA\right)\hat{\boldsymbol \varepsilon }_t= \hat {\boldsymbol D}_t \hat{\boldsymbol \varepsilon }_t  </math>
319
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
320
|}
321
|}
322
323
In vectors <math display="inline">\hat{\boldsymbol \sigma }_p</math> and <math display="inline">\hat{\boldsymbol \sigma }_t</math>, <math display="inline">N,M</math> and <math display="inline">Q</math> are respectively the axial force, the bending moment and the transverse shear force of standard beam theory, whereas <math display="inline">M_\phi </math> and <math display="inline">Q_\phi </math> are an additional bending moment and an additional shear force which are conjugate to the new generalized strains <math display="inline">{\partial \Psi  \over \partial x}</math> and <math display="inline">\Psi </math>, respectively.
324
325
The generalized constitutive matrices <math display="inline">\hat {\boldsymbol D}_p</math> and <math display="inline">\hat {\boldsymbol D}_t</math> are
326
327
{| class="formulaSCP" style="width: 100%; text-align: left;" 
328
|-
329
| 
330
{| style="text-align: left; margin:auto;width: 100%;" 
331
|-
332
| style="text-align: center;" | <math>\hat{\boldsymbol D}_p =\iint _A E^k \left[\begin{matrix}1 & -z &  \phi ^k\\ -z &z^2&-z\phi ^k\\ \phi ^k & - z\phi ^k & (\phi ^k)^2 \end{matrix}\right]dA \quad ,\quad  \hat{\boldsymbol D}_t = \left[\begin{matrix}D_s &-g\\ -g & g \end{matrix}\right] </math>
333
|}
334
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.a)
335
|}
336
337
with
338
339
{| class="formulaSCP" style="width: 100%; text-align: left;" 
340
|-
341
| 
342
{| style="text-align: left; margin:auto;width: 100%;" 
343
|-
344
| style="text-align: center;" | <math>D_s =\iint _A G^k dA \quad ,\quad g = D_s - GA </math>
345
|}
346
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.b)
347
|}
348
349
In the derivation of the expression for <math display="inline">\hat {\boldsymbol D}_s</math> we have used the definition of <math display="inline">\beta ^k</math> of Eq.(9).
350
351
The generalized constitutive equation can be written as
352
353
{| class="formulaSCP" style="width: 100%; text-align: left;" 
354
|-
355
| 
356
{| style="text-align: left; margin:auto;width: 100%;" 
357
|-
358
| style="text-align: center;" | <math>\hat{\boldsymbol \sigma } = \left\{\begin{array}{c}\hat{\boldsymbol \sigma }_p \\                             \hat{\boldsymbol \sigma }_t                           \end{array}  \right\}= \hat {\boldsymbol D}  \hat{\boldsymbol \varepsilon } = \hat {\boldsymbol D} \left\{\begin{array}{c}\hat{\boldsymbol \varepsilon }_p \\                             \hat{\boldsymbol \varepsilon }_t                           \end{array} \right\}\quad \quad \hbox{with } \quad \hat {\boldsymbol D}=\left[\begin{matrix}\hat{\boldsymbol D}_p & {\boldsymbol 0}\\ {\boldsymbol 0}&\hat{\boldsymbol D}_t  \end{matrix}\right] </math>
359
|}
360
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
361
|}
362
363
===3.4 Virtual work expression===
364
365
The virtual work expression for a distributed load <math display="inline">q</math> is
366
367
{| class="formulaSCP" style="width: 100%; text-align: left;" 
368
|-
369
| 
370
{| style="text-align: left; margin:auto;width: 100%;" 
371
|-
372
| style="text-align: center;" | <math>\iiint _V (\delta \varepsilon _x^k \sigma _x^k + \delta \gamma _{xz}^k \tau _{xz}^k)dV - \int _l \delta w q dA =0 </math>
373
|}
374
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
375
|}
376
377
The l.h.s. of Eq.(17) contains the internal virtual work performed by the axial and tangential stresses and the r.h.s. is the external virtual work carried out by the distributed load. <math display="inline">V</math> and <math display="inline">l</math> are the volume and length of the beam, respectively.
378
379
Substituting Eqs.(3) into the expression for the virtual internal work  and using Eqs.(13) and (14) gives
380
381
{| class="formulaSCP" style="width: 100%; text-align: left;" 
382
|-
383
| 
384
{| style="text-align: left; margin:auto;width: 100%;" 
385
|-
386
| style="text-align: center;" | <math>\iiint _V \left(\delta \varepsilon _x^k \sigma _x^k + \delta \gamma _{xz}^k \tau _{xz}^k\right)dV = \!\!\iiint _V \left(\delta  \hat{\boldsymbol \varepsilon }_p^T {\boldsymbol S}_p^T \sigma _x^k + \delta \hat{\boldsymbol \varepsilon }_t^T {\boldsymbol S}_t^T \tau _{xz}^k \right)dV =</math>
387
|-
388
| style="text-align: center;" | <math> =\!\! \int _l \left(\delta \hat{\boldsymbol \varepsilon }_p^T \hat{\boldsymbol \sigma }_p + \delta \hat{\boldsymbol \varepsilon }_t^T \hat{\boldsymbol \sigma }_t\right)dx </math>
389
|}
390
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
391
|}
392
393
The virtual work is therefore written as
394
395
{| class="formulaSCP" style="width: 100%; text-align: left;" 
396
|-
397
| 
398
{| style="text-align: left; margin:auto;width: 100%;" 
399
|-
400
| style="text-align: center;" | <math>\int _l \left(\delta \hat{\boldsymbol \varepsilon }_p^T\hat{\boldsymbol \sigma }_p + \delta \hat{\boldsymbol \varepsilon }_t^T \hat{\boldsymbol \sigma }_t\right)dx - \int _l \delta w q dx=0 </math>
401
|}
402
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
403
|}
404
405
==4 TWO-NODED LRZ BEAM ELEMENT==
406
407
The  four kinematic variables are <math display="inline">u_0,w_0, \theta </math> and <math display="inline">\Psi </math>. They can be discretized using  2-noded linear <math display="inline">C^\circ </math> beam elements of length <math display="inline">l^e</math> in the standard form as
408
409
{| class="formulaSCP" style="width: 100%; text-align: left;" 
410
|-
411
| 
412
{| style="text-align: left; margin:auto;width: 100%;" 
413
|-
414
| style="text-align: center;" | <math>{\boldsymbol u} = \left\{                     \begin{array}{c}u_0 \\                       w_0 \\                       \theta \\                      \Psi \\                     \end{array}  \right\}= \sum \limits _{i=1}^2 N_i {\boldsymbol a}_i = {\boldsymbol N} {\boldsymbol a}^e  </math>
415
|}
416
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
417
|}
418
419
with
420
421
{| class="formulaSCP" style="width: 100%; text-align: left;" 
422
|-
423
| 
424
{| style="text-align: left; margin:auto;width: 100%;" 
425
|-
426
| style="text-align: center;" | <math>{\boldsymbol N}= [N_1{\boldsymbol I}_4,N_2{\boldsymbol I}_4]~~,~~{\boldsymbol a}^e = \left\{                     \begin{array}{c}{\boldsymbol a}_1\\{\boldsymbol a}_2\end{array}  \right\}~~,~~ {\boldsymbol a}_i= \left\{                     \begin{array}{c}u_{0_i} \\                       w_{0_i} \\                       \theta _i \\                      \Psi _i \\                     \end{array}  \right\} </math>
427
|}
428
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
429
|}
430
431
where <math display="inline">N_i = \frac{1}{2} (1+\xi \xi _i)</math> with <math display="inline">\xi = 1 - \frac{2x}{l^e}</math> are the standard one-dimensional linear shape functions, <math display="inline">{\boldsymbol a}_i</math> is the vector of nodal kinematic variables and <math display="inline">{\boldsymbol I}_4</math> is the <math display="inline">4\times 4</math> unit matrix.
432
433
Substituting Eq.(20) into the generalized strain vectors in Eq.(2c) gives
434
435
{| class="formulaSCP" style="width: 100%; text-align: left;" 
436
|-
437
| 
438
{| style="text-align: left; margin:auto;width: 100%;" 
439
|-
440
| style="text-align: center;" | <math>\hat {\boldsymbol \varepsilon }_p = {\boldsymbol B}_p {\boldsymbol a}^e \quad ,\quad \hat {\boldsymbol \varepsilon }_t = {\boldsymbol B}_t {\boldsymbol a}^e </math>
441
|}
442
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
443
|}
444
445
The  generalized strain matrices <math display="inline">{\boldsymbol B}_p</math> and <math display="inline">{\boldsymbol B}_t </math> are
446
447
{| class="formulaSCP" style="width: 100%; text-align: left;" 
448
|-
449
| 
450
{| style="text-align: left; margin:auto;width: 100%;" 
451
|-
452
| style="text-align: center;" | <math>{\boldsymbol B}_p =[{\boldsymbol B}_{p_1},{\boldsymbol B}_{p_2} ]\quad ,\quad {\boldsymbol B}_t =[{\boldsymbol B}_{t_1},{\boldsymbol B}_{t_2} ] </math>
453
|}
454
| style="width: 5px;text-align: right;white-space: nowrap;" | (23.a)
455
|}
456
457
with
458
459
{| class="formulaSCP" style="width: 100%; text-align: left;" 
460
|-
461
| 
462
{| style="text-align: left; margin:auto;width: 100%;" 
463
|-
464
| style="text-align: center;" | <math>{\boldsymbol B}_{p_i} = \begin{bmatrix}\displaystyle \frac{\partial N_i}{\partial x} & 0 & 0 & 0 \\ 0 & 0 & \displaystyle \frac{\partial N_i}{\partial x} & 0 \\ 0 & 0 & 0 & \displaystyle \frac{\partial N_i}{\partial x} \end{bmatrix} </math>
465
|}
466
| style="width: 5px;text-align: right;white-space: nowrap;" | (23.b)
467
|}
468
469
{| class="formulaSCP" style="width: 100%; text-align: left;" 
470
|-
471
| 
472
{| style="text-align: left; margin:auto;width: 100%;" 
473
|-
474
| style="text-align: center;" | <math>\displaystyle {\boldsymbol B}_{t_i} = \begin{bmatrix}0&\displaystyle \frac{\partial N_i}{\partial x}&-N_i & 0\\ --&--&--&--\\ 0&0&0&N_i \end{bmatrix}= \begin{bmatrix}{\boldsymbol B}_{s_i}\\ --\\ {\boldsymbol B}_{\psi _i}   \end{bmatrix} </math>
475
|}
476
| style="width: 5px;text-align: right;white-space: nowrap;" | (23.c)
477
|}
478
479
where <math display="inline">{\boldsymbol B}_{p_i}</math> and <math display="inline">{\boldsymbol B}_{t_i}</math> are the in-plane and transverse shear strain matrices for node <math display="inline">i</math>.
480
481
The virtual displacement and generalized strain fields are expressed in terms of the virtual nodal kinematic variables as
482
483
{| class="formulaSCP" style="width: 100%; text-align: left;" 
484
|-
485
| 
486
{| style="text-align: left; margin:auto;width: 100%;" 
487
|-
488
| style="text-align: center;" | <math>\delta {\boldsymbol u} = {\boldsymbol N} \delta {\boldsymbol a}^e~~,~~ \delta \hat{\boldsymbol \varepsilon }_p={\boldsymbol B}_p \delta {\boldsymbol a}^e ~~,~~ \delta \hat{\boldsymbol \varepsilon }_t ={\boldsymbol B}_t \delta {\boldsymbol a}^e </math>
489
|}
490
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
491
|}
492
493
The discretized equilibrium equations are obtained by substituting Eqs.(13), (14), (20), (22) and (24) into the virtual work expression (19). After simplification of the virtual nodal kinematic variables, the following standard matrix equation is obtained
494
495
{| class="formulaSCP" style="width: 100%; text-align: left;" 
496
|-
497
| 
498
{| style="text-align: left; margin:auto;width: 100%;" 
499
|-
500
| style="text-align: center;" | <math>{\boldsymbol K}{\boldsymbol a}-{\boldsymbol f}=0 </math>
501
|}
502
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
503
|}
504
505
where <math display="inline">{\boldsymbol a}</math> is the vector of nodal kinematic variables for the  whole mesh.
506
507
The stiffness matrix <math display="inline">{\boldsymbol K}</math> and the equivalent nodal force vector <math display="inline">{\boldsymbol f}</math> are obtained by assembling the element contributions <math display="inline">{\boldsymbol K}^e</math> and <math display="inline">{\boldsymbol f}^e</math> given by
508
509
{| class="formulaSCP" style="width: 100%; text-align: left;" 
510
|-
511
| 
512
{| style="text-align: left; margin:auto;width: 100%;" 
513
|-
514
| style="text-align: center;" | <math>{\boldsymbol K}^e = {\boldsymbol K}^e_p + {\boldsymbol K}^e_t </math>
515
|}
516
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
517
|}
518
519
with
520
521
{| class="formulaSCP" style="width: 100%; text-align: left;" 
522
|-
523
| 
524
{| style="text-align: left; margin:auto;width: 100%;" 
525
|-
526
| style="text-align: center;" | <math>{\boldsymbol K}^e_{p_{ij}} = \int _{l^e} {\boldsymbol B}^T_{p_i} \hat {\boldsymbol D}_p {\boldsymbol B}_{p_j}dx ~~,~~ {\boldsymbol K}^e_{t_{ij}} = \int _{l^e} {\boldsymbol B}^T_{t_i} \hat {\boldsymbol D}_t {\boldsymbol B}_{t_j}dx </math>
527
|}
528
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
529
|}
530
531
and
532
533
{| class="formulaSCP" style="width: 100%; text-align: left;" 
534
|-
535
| 
536
{| style="text-align: left; margin:auto;width: 100%;" 
537
|-
538
| style="text-align: center;" | <math>{\boldsymbol f}^e_i = \int _{l^e} N_i q [1,0,0,0]^T dx </math>
539
|}
540
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
541
|}
542
543
Matrix <math display="inline">{\boldsymbol K}^e_p</math> is integrated with a one-point numerical quadrature which is exact in this case. Full integration of matrix <math display="inline">{\boldsymbol K}^e_t</math> requires a two-point Gauss quadrature. This however leads to shear locking for slender composite laminated beams (Section 5).
544
545
In order to asses the influence of the reduced integration of matrix <math display="inline">{\boldsymbol K}^e_t</math> for overcoming the  shear locking problem we split <math display="inline">{\boldsymbol K}^e_t</math> as follows
546
547
{| class="formulaSCP" style="width: 100%; text-align: left;" 
548
|-
549
| 
550
{| style="text-align: left; margin:auto;width: 100%;" 
551
|-
552
| style="text-align: center;" | <math>{\boldsymbol K}^e_t = {\boldsymbol K}^e_s +{\boldsymbol K}^e_\psi + {\boldsymbol K}^e_{s\psi } + [{\boldsymbol K}^e_{s\psi }]^T </math>
553
|}
554
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.a)
555
|}
556
557
with
558
559
{| class="formulaSCP" style="width: 100%; text-align: left;" 
560
|-
561
| 
562
{| style="text-align: left; margin:auto;width: 100%;" 
563
|-
564
| style="text-align: center;" | <math>{\boldsymbol K}^e_{s_{ij}} = \int _{l^e}  D_s {\boldsymbol B}^T_{s_i} {\boldsymbol B}_{s_j}dx ~~,~~ {\boldsymbol K}^e_{\psi _{ij}} = \int _{l^e} g {\boldsymbol B}_{\psi _i}^T {\boldsymbol B}_{\psi _j}dx </math>
565
|}
566
| style="width: 5px;text-align: right;white-space: nowrap;" | (29.b)
567
|}
568
569
{| class="formulaSCP" style="width: 100%; text-align: left;" 
570
|-
571
| 
572
{| style="text-align: left; margin:auto;width: 100%;" 
573
|-
574
| style="text-align: center;" | <math> {\boldsymbol K}^e_{s\psi _{ij}}= \int _{l^e} (-g ) {\boldsymbol B}^T_{s_i}{\boldsymbol B}_{\psi _j}dx </math>
575
|}
576
|}
577
578
where <math display="inline">{\boldsymbol B}_{s_i}</math> and <math display="inline">{\boldsymbol B}_{\psi _i}</math> are defined in Eq.(23c) and <math display="inline">D_s</math> and <math display="inline">g</math> are given in Eq.(15b).
579
580
The new linear beam element based on the RZ theory is termed LRZ.
581
582
A study of the accuracy of the LRZ beam element for analysis of slender laminated beams using one and two-point quadratures for integrating matrices <math display="inline">{\boldsymbol K}^e_s</math>, <math display="inline">{\boldsymbol K}^e_\psi </math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>  is  presented in the next section.
583
584
==5 STUDY OF SHEAR LOCKING FOR THE LRZ BEAM ELEMENT==
585
586
We study the performance of the  LRZ beam element for the analysis of a cantilever beam of length <math display="inline">L</math> under an end point load of value <math display="inline">F=1</math> (Figure [[#img-2|2]]). The beam is formed by a symmetric three-layered material whose properties are listed on Table 1. The analysis is performed for four span-to-thickness ratios: <math display="inline">\lambda = 5,10,50,100</math> (<math display="inline">\lambda = L/h</math>) using a mesh of 100 LRZ beam elements. Results for the LRZ element are labelled “ZZ” in the figures.
587
588
The same beam was analized using a mesh of 27000 four-noded plane stress  rectangles for comparison purposes (Figure [[#img-3|3]]).  Results for the plane stress analysis are labeled “PS” in the figures.
589
590
<div id='img-2'></div>
591
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
592
|-
593
|[[Image:Draft_Samper_624436055-test-Fig2.png|600px|Cantilever beam under point  load ]]
594
|- style="text-align: center; font-size: 75%;"
595
| colspan="1" | '''Figure 2:''' Cantilever beam under point  load 
596
|}
597
598
<div id='img-3'></div>
599
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
600
|-
601
|[[Image:Draft_Samper_624436055-test-Fig3.png|600px|Mesh of 27000 4-noded plane stress rectangular elements for analysis of cantilever and simple supported beams]]
602
|- style="text-align: center; font-size: 75%;"
603
| colspan="1" | '''Figure 3:''' Mesh of 27000 4-noded plane stress rectangular elements for analysis of cantilever and simple supported beams
604
|}
605
606
607
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
608
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Symmetric 3-layered cantilever beam. Material properties for shear locking study
609
|- style="border-top: 2px solid;"
610
| style="text-align: left;" |   
611
| colspan='3' style="text-align: left;" | '''Composite material properties'''
612
|- style="border-top: 2px solid;"
613
| style="text-align: left;" |   
614
| Layer 1 
615
| Layer 2  
616
| Layer 3 
617
|-
618
| style="text-align: left;" | 
619
| (bottom)
620
| (core) 
621
| (top)
622
|- style="border-top: 2px solid;"
623
| style="text-align: left;" |  h [mm] 
624
| 6.6667 
625
| 6.6667 
626
| 6.6667 
627
|-
628
| style="text-align: left;" | E [MPa] 
629
| 2.19E5 
630
| 2.19E3 
631
| 2.19E5 
632
|- style="border-bottom: 2px solid;"
633
| style="text-align: left;" | G [MPa] 
634
| 0.876E5 
635
| 8.80E2 
636
| 0.876E5 
637
638
|}
639
640
Figure [[#img-4|4]] shows the ratio <math display="inline">r</math> between the end node deflection  obtained with the LRZ element (<math display="inline">w_{zz}</math>) and with the plane stress quadrilateral (<math display="inline">w_{ps}</math>) (i.e. <math display="inline">r=\frac{w_{zz}}{w_{ps}}</math>) versus the beam span-to-thickness ratio <math display="inline">d = \frac{L}{h}</math>. Results for the LRZ element have been obtained using ''exact'' two-point integration for all terms of  matrix <math display="inline">{\boldsymbol K}^e_t</math> (Eq.(27)) and  a one-point ''reduced'' integration for the following three groups of matrices: <math display="inline">{\boldsymbol K}^e_s</math>; <math display="inline">{\boldsymbol K}^e_s</math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>; and <math display="inline">{\boldsymbol K}_s</math>, <math display="inline">{\boldsymbol K}^e_{s\psi }</math> and <math display="inline">{\boldsymbol K}^e_{\psi }</math> (Eqs.(29b)).
641
642
Labels ``all'<nowiki/>''', ``S'<nowiki/>''', ``SPsi'<nowiki/>''' and `''`''Psi'''' in Figures [[#img-4|4]]&#8211;[[#img-7|7]] refer to matrices <math display="inline">{\boldsymbol K}^e_t</math>, <math display="inline">{\boldsymbol K}^e_s</math>, <math display="inline">{\boldsymbol K}^e_{s\psi }</math> and <math display="inline">{\boldsymbol K}^e_{\psi }</math> of Eq.(29a), respectively.
643
644
Results in Figure [[#img-4|4]] clearly show that the exact integration of <math display="inline">{\boldsymbol K}^e_t</math> leads to shear locking as expected. Good (locking-free) results are obtained by one-point reduced integration of the three groups of matrices.
645
646
The influence of reduced integration   in the distribution of the transverse shear stress was studied next for the three groups of matrices. Figures [[#img-5|5]]&#8211;[[#img-7|7]] show the thickness distribution of <math display="inline">\tau _{xz}</math> in sections located at distances <math display="inline">\frac{L}{20},\frac{L}{4},\frac{L}{2}</math> and <math display="inline">\frac{3}{4}L</math> from the clamped end for span-to-thickness ratios of <math display="inline">\lambda =5 , 10</math> and 100. Results are compared with the plane stress solution and also with results obtained with a mesh of 300 standard 2-noded elements based on laminated Timoshenko beam theory (labelled TBT in the figures). All TBT results presented in the paper have been used with a simple shear correction factor of <math display="inline">\frac{5}{6}</math>. Indeed a more accurate value of the shear correction factor in TBT can be used for laminated sections <span id='citeF-28'></span>[[#cite-28|[28]]].
647
648
<div id='img-4'></div>
649
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
650
|-
651
|[[Image:Draft_Samper_624436055-test-Fig4a.png|600px|r ratio \left(r = \fracw<sub>zz</sub>wₚₛ\right) versus L/h for cantilever under point load analyzed with the LRZ element. Labels ``''all'''', S, SPsi and Psi refer to matrices K<sup>e</sup>ₜ, Kₛ<sup>e</sup>,K<sub>sψ</sub><sup>e</sup> and K<sub>ψ</sub><sup>e</sup>, respectively]]
652
|- style="text-align: center; font-size: 75%;"
653
| colspan="1" | '''Figure 4:'<nowiki/>'' <math>r</math> ratio <math>\left(r = \frac{w_{zz}}{w_{ps}}\right)</math> versus <math>L/h</math> for cantilever under point load analyzed with the LRZ element. Labels ``''all'''', <math>S</math>, <math>SPsi</math> and <math>Psi</math> refer to matrices <math>{\boldsymbol K}^e_t</math>, <math>{\boldsymbol K}_s^e,{\boldsymbol K}_{s\psi }^e</math> and <math>{\boldsymbol K}_{\psi }^e</math>, respectively
654
|}
655
656
<div id='img-5a'></div>
657
<div id='img-5b'></div>
658
<div id='img-5c'></div>
659
<div id='img-5d'></div>
660
<div id='img-5'></div>
661
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
662
|-
663
|[[Image:Draft_Samper_624436055-test-Grafico2.png|360px|τ<sub>xz</sub>,λ= 5 , L/20]]
664
|[[Image:Draft_Samper_624436055-test-Grafico3.png|360px|τ<sub>xz</sub> ,λ= 5 , L/4]]
665
|- style="text-align: center; font-size: 75%;"
666
| (a) <math>\tau _{xz},\lambda = 5 , L/20</math>
667
| (b) <math>\tau _{xz} ,\lambda = 5 , L/4</math>
668
|-
669
|[[Image:Draft_Samper_624436055-test-Grafico4.png|360px|τ<sub>xz</sub> ,λ= 5 , L/2]]
670
|[[Image:Draft_Samper_624436055-test-Grafico5.png|360px|τ<sub>xz</sub> ,λ= 5 , 3L/4]]
671
|- style="text-align: center; font-size: 75%;"
672
| (c) <math>\tau _{xz} ,\lambda = 5 , L/2</math>
673
| (d) <math>\tau _{xz} ,\lambda = 5 , 3L/4</math>
674
|- style="text-align: center; font-size: 75%;"
675
| colspan="2" | '''Figure 5:''' Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of shear stress for <math>\lambda =5</math> at different sections
676
|}
677
678
<div id='img-6a'></div>
679
<div id='img-6b'></div>
680
<div id='img-6c'></div>
681
<div id='img-6d'></div>
682
<div id='img-6'></div>
683
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
684
|-
685
|[[Image:Draft_Samper_624436055-test-Grafico6.png|360px|τ<sub>xz</sub> ,λ= 10 , L/20]]
686
|[[Image:Draft_Samper_624436055-test-Grafico7.png|360px|τ<sub>xz</sub> ,λ= 10 , L/4]]
687
|- style="text-align: center; font-size: 75%;"
688
| (a) <math>\tau _{xz} ,\lambda = 10 , L/20</math>
689
| (b) <math>\tau _{xz} ,\lambda = 10 , L/4</math>
690
|-
691
|[[Image:Draft_Samper_624436055-test-Grafico8.png|360px|τ<sub>xz</sub> ,λ= 10 , L/2]]
692
|[[Image:Draft_Samper_624436055-test-Grafico9.png|360px|τ<sub>xz</sub> ,λ= 10 , 3L/4]]
693
|- style="text-align: center; font-size: 75%;"
694
| (c) <math>\tau _{xz} ,\lambda = 10 , L/2</math>
695
| (d) <math>\tau _{xz} ,\lambda = 10 , 3L/4</math>
696
|- style="text-align: center; font-size: 75%;"
697
| colspan="2" | '''Figure 6:''' Symmetric 3-layered cantilever thick beam under end point load.  Thickness distribution of shear stress for <math>\lambda =10</math> at different sections
698
|}
699
700
<div id='img-7a'></div>
701
<div id='img-7b'></div>
702
<div id='img-7c'></div>
703
<div id='img-7d'></div>
704
<div id='img-7'></div>
705
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
706
|-
707
|[[Image:Draft_Samper_624436055-test-Grafico10.png|360px|τ<sub>xz</sub> ,λ= 100 , L/20]]
708
|[[Image:Draft_Samper_624436055-test-Grafico11.png|360px|τ<sub>xz</sub> ,λ= 100 , L/4]]
709
|- style="text-align: center; font-size: 75%;"
710
| (a) <math>\tau _{xz} ,\lambda = 100 , L/20</math>
711
| (b) <math>\tau _{xz} ,\lambda = 100 , L/4</math>
712
|-
713
|[[Image:Draft_Samper_624436055-test-Grafico12.png|360px|τ<sub>xz</sub> ,λ= 100 , L/2]]
714
|[[Image:Draft_Samper_624436055-test-Grafico13.png|360px|τ<sub>xz</sub> ,λ= 100 , 3L/4]]
715
|- style="text-align: center; font-size: 75%;"
716
| (c) <math>\tau _{xz} ,\lambda = 100 , L/2</math>
717
| (d) <math>\tau _{xz} ,\lambda = 100 , 3L/4</math>
718
|- style="text-align: center; font-size: 75%;"
719
| colspan="2" | '''Figure 7:''' Symmetric 3-layered cantilever thick beam under end point load. Thickness distribution of transverse shear stress for <math>\lambda =100</math> at different sections
720
|}
721
722
The conclusion is that for small values of <math display="inline">\lambda </math> the reduced or exact reduced integration of matrix <math display="inline">{\boldsymbol K}^e_t</math> leads to similar results. For slender beams, however, results obtained using reduced integration for <math display="inline">{\boldsymbol K}^e_s</math>; <math display="inline">{\boldsymbol K}^e_s</math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>; and <math display="inline">{\boldsymbol K}^e_s</math>, <math display="inline">{\boldsymbol K}^e_{s\psi }</math> and <math display="inline">{\boldsymbol K}^e_{\psi }</math> are different. Slightly more accurate results are obtained with the second choice for the section at <math display="inline">x=L/4</math> and <math display="inline">\lambda =100</math> (Figure [[#img-7|7]]b).
723
724
In conclusion, we recommend using a reduced one-point integration for matrices <math display="inline">{\boldsymbol K}^e_s</math> and <math display="inline">{\boldsymbol K}^e_{s\psi }</math>, while matrix <math display="inline">{\boldsymbol K}^e_\psi </math> should be integrated with a 2-point quadrature.
725
726
==6 CONVERGENCE STUDY==
727
728
The same three-layered cantilever beam of Figure [[#img-2|2]] was studied next for three different set of thickness and material properties for the three layers as listed in Table [[#table-2|2]]. Material A is the more homogeneous one, while material C is clearly the more heterogeneous.
729
730
The problem was studied with  six meshes of LRZ elements ranging from 5 to 300 elements. Tables [[#table-3|3]]&#8211;[[#table-5|5]] show the convergence with the number of elements for the deflection and function <math display="inline">\Psi </math> at the beam end,  the maximum axial stress <math display="inline">\sigma _x</math> at the end section and the maximum shear stress <math display="inline">\tau _{xz}</math> at the mid section.
731
732
733
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
734
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Non symmetric 3-layered cantilever beams. Material properties for convergence analysis
735
736
|-
737
| 
738
| 
739
| colspan='3' | '''Material properties'''
740
741
|-
742
| 
743
| 
744
| '''Layer 1'''(bottom) 
745
| '''Layer 2''' (core) 
746
| '''Layer 3''' (top) 
747
748
|- style="border-top: 2px solid;"
749
| '''Composite A''' 
750
| h [mm] 
751
| 6.66 
752
| 6.66 
753
| 6.66 
754
|-
755
|  
756
| E [MPa] 
757
| 4.40E5 
758
| 2.19E4 
759
| 2.19E5 
760
|-
761
|  
762
| G [MPa] 
763
| 2.00E5 
764
| 8.80E3 
765
| 8.76E4 
766
767
|- style="border-top: 2px solid;"
768
| '''Composite B''' 
769
| h [mm] 
770
| 6.66 
771
| 6.66 
772
| 6.66 
773
|-
774
|  
775
| E [MPa] 
776
| 2.19E5 
777
| 2.19E3 
778
| 2.19E5 
779
|-
780
|  
781
| G [MPa] 
782
| 8.76E4 
783
| 8.80E2 
784
| 8.76E4 
785
786
|- style="border-top: 2px solid;"
787
| '''Composite C''' 
788
| h [mm] 
789
| 2 
790
| 16 
791
| 2 
792
|-
793
| 
794
| E [MPa] 
795
| 7.30E5 
796
| 7.30E2 
797
| 2.19E5 
798
|- style="border-bottom: 2px solid;"
799
| 
800
| G [MPa] 
801
| 2.92E5 
802
| 2.20E2 
803
| 8.76E4  
804
805
|}
806
807
808
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
809
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Non symmetric 3-layered cantilever thick beams under end point load (<math>\lambda =5</math>). Relative error for <math>w</math> at <math>x=L</math>
810
811
|- style="border-top: 2px solid;"
812
| colspan='4' | <math>e_r %  - w</math> at <math>x=L</math>
813
814
|- style="border-top: 2px solid;"
815
| '''Number of''' 
816
| colspan='3' | '''Composites'''
817
818
|-
819
| '''elements''' 
820
| '''A''' 
821
| '''B''' 
822
| '''C''' 
823
824
|- style="border-top: 2px solid;"
825
| 5 
826
| 1.800 
827
| 9.588 
828
| 42.289 
829
|-
830
| 10 
831
| 0.506 
832
| 2.901 
833
| 19.277 
834
|-
835
| 25 
836
| 0.0860 
837
| 0.499 
838
| 4.913 
839
|-
840
| 50 
841
| 0.0191 
842
| 0.123 
843
| 1.406 
844
|-
845
| 100 
846
| 0.0048 
847
| 0.031 
848
| 0.339 
849
|- style="border-bottom: 2px solid;"
850
| 300 
851
| 0.0000 
852
| 0.0000 
853
| 0.0000 
854
855
|}
856
857
858
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
859
|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 Non symmetric 3-layered cantilever thick beams under end point load (<math>\lambda =5</math>). Convergence study. Relative error for <math>\Psi </math> at <math>x=L</math>
860
861
|- style="border-top: 2px solid;"
862
| colspan='4' | <math>e_r %  - \Psi </math> at <math>x=L</math>
863
864
|- style="border-top: 2px solid;"
865
| '''Number of''' 
866
| colspan='3' | '''Composites'''
867
868
|-
869
| '''elements''' 
870
| '''A''' 
871
| '''B''' 
872
| '''C''' 
873
874
|- style="border-top: 2px solid;"
875
| 5 
876
| 0.040 
877
| 8.563 
878
| 36.113 
879
|-
880
| 10 
881
| 0.003
882
| 1.814
883
| 8.042 
884
|-
885
| 25 
886
| 0.000
887
| 0.259 
888
| 0.328 
889
|-
890
| 50 
891
| 0.000
892
| 0.063 
893
| 0.033 
894
|-
895
| 100 
896
| 0.000 
897
| 0.016 
898
| 0.007 
899
|- style="border-bottom: 2px solid;"
900
| 300 
901
| 0.000 
902
| 0.000 
903
| 0.000 
904
905
|}
906
907
908
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
909
|+ style="font-size: 75%;" |<span id='table-5'></span>Table. 5 Non symmetric 3-layered cantilever thick beams under end  point load (<math>\lambda =5</math>). Convergence study. (a) Relative error for the maximum value of <math>\sigma _x</math> at <math>x=L</math> and (b) idem for <math>\tau _{xz}</math> at <math>x=L/2</math> 
910
|-
911
| colspan='4' | (a)
912
913
|- style="border-top: 2px solid;"
914
| colspan='4' | <math>e_r %  - (\sigma _x)_{\max }</math> at <math>x=L</math>
915
916
|- style="border-top: 2px solid;"
917
| '''Number of''' 
918
| colspan='3' | '''Composites'''
919
920
|-
921
| '''elements''' 
922
| '''A''' 
923
| '''B''' 
924
| '''C''' 
925
926
|- style="border-top: 2px solid;"
927
| 5 
928
| 0.568 
929
| 6.923 
930
| 18.239 
931
|-
932
| 10 
933
| 0.076 
934
| 2.704 
935
| 12.437 
936
|-
937
| 25 
938
| 0.013 
939
| 0.568 
940
| 4.266 
941
|-
942
| 50 
943
| 0.003 
944
| 0.131 
945
| 1.095 
946
|-
947
| 100 
948
| 0.001 
949
| 0.029 
950
| 0.250 
951
|- style="border-bottom: 2px solid;"
952
| 300 
953
| 0.000 
954
| 0.000 
955
| 0.000 
956
957
|}
958
959
960
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
961
|-
962
| colspan='4' | (b)
963
964
|- style="border-top: 2px solid;"
965
| colspan='4' | <math>e_r %  - (\tau _{xz})_{\max }</math> at <math>\frac{L}{2}</math>
966
967
|- style="border-top: 2px solid;"
968
| '''Number of''' 
969
| colspan='3' | '''Composites'''
970
971
|-
972
| '''elements''' 
973
| '''A''' 
974
| '''B''' 
975
| '''C''' 
976
977
|- style="border-top: 2px solid;"
978
| 5 
979
| 7.020 
980
| 19.283 
981
| 50.938 
982
|-
983
| 10 
984
| 0.352 
985
| 5.176 
986
| 20.602 
987
|-
988
| 25 
989
| 0.052 
990
| 0.888 
991
| 3.408 
992
|-
993
| 50 
994
| 0.010 
995
| 0.210 
996
| 0.707 
997
|-
998
| 100 
999
| 0.003 
1000
| 0.049 
1001
| 0.147 
1002
|- style="border-bottom: 2px solid;"
1003
| 300 
1004
| 0.000 
1005
| 0.000 
1006
| 0.000 
1007
1008
|}
1009
1010
Convergence is measured by the relative error defined (in absolute value) as
1011
1012
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1013
|-
1014
| 
1015
{| style="text-align: left; margin:auto;width: 100%;" 
1016
|-
1017
| style="text-align: center;" | <math>e_r = \left|\frac{v_6-v_i}{v_6}\right| </math>
1018
|}
1019
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
1020
|}
1021
1022
where <math display="inline">v_6</math> and <math display="inline">v_i</math> are the values of the magnitude of interest obtained using the finest grid  (300 elements) and the <math display="inline">i</math>th mesh (<math display="inline">i=1,2,\cdots 5</math>), respectively.
1023
1024
Results clearly show that convergence is always slower for the heterogeneous material case, as expected.
1025
1026
For a mesh of 25 elements the errors for all the magnitudes considered are less than 1% for materials A and B. For material C the maximum error does not exceed 5% (Table [[#table-5|5]]). For the  50 element mesh errors of the order of 1% or less were obtained in all cases.
1027
1028
Results for a 10 element mesh are good for material A (errors less than 0.4%), relatively good for material B (errors less than around <math display="inline">5%</math>) and unacceptable for material C (errors ranging from around 8% to 20%
1029
1030
==7 EXAMPLES OF APPLICATION==
1031
1032
===7.1 Three-layered thick cantilever beam with non symmetric material properties===
1033
1034
We present results for  a laminated thick cantilever beam under an end point load. The material properties are those of Composite C in Table [[#table-2|2]]. The span-to-thickness ratio is <math display="inline">\lambda =5</math>.
1035
1036
For the laminated sandwich considered the core  is eight times thicker than the face sheets. In addition, the core is three orders of magnitude more compliant than the bottom face sheet. Moreover, the top face sheet has the same thickness as the bottom face sheet, but is about three times stiffer.   Note that this laminate does not possess material symmetry with respect to the mid-depth reference axis. The high heterogeneity of this stacking sequence is very challenging for the beam theories considered herein to model adequately.
1037
1038
As in previous section, the legend caption PS denotes  the ''reference'' solution obtained with the structured mesh  of 27000 four-noded plane stress quadrilaterals shown in  Figure [[#img-3|3]]. TBT denotes the solution obtained with a mesh of 300 2-noded beam elements based on standard laminated Timoshenko beam theory. LRZ-300, LRZ-50, LRZ-25, LRZ-10 refer to the solution obtained with the LRZ beam element using meshes of 300, 50, 25 and 10 elements, respectively.
1039
1040
Figure [[#img-8|8]] shows the deflection values along the beam length. Very good agreement with the plane stress solution is obtained already for the LRZ-50 mesh as expected from the conclusions of the previous section.
1041
1042
TBT results are considerable stiffer. The difference with the reference solution is about ''six times stiffer'' for the end deflection value.<br/>
1043
1044
Figure [[#img-9|9]] shows the distribution of the axial displacements at the upper and lower surfaces of layer 3 (top layer) along the beam length. Excellent results are again obtained with the 50 element mesh. The TBT results are far from the correct ones.
1045
1046
<div id='img-8'></div>
1047
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1048
|-
1049
|[[File:Draft_Samper_624436055_2824_Fig8.jpg|420px|]]
1050
|- style="text-align: center; font-size: 75%;"
1051
| colspan="1" | '''Figure 8:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Distribution of the vertical deflection <math>w</math> for different theories and meshes
1052
|}
1053
1054
<div id='img-9a'></div>
1055
<div id='img-9b'></div>
1056
<div id='img-9'></div>
1057
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1058
|-
1059
|[[Image:Draft_Samper_624436055-test-Fig9a.png|420px|]]
1060
|[[Image:Draft_Samper_624436055-test-Fig9b.png|420px|]]
1061
|- style="text-align: center; font-size: 75%;"
1062
| (a) 
1063
| (b) 
1064
|- style="text-align: center; font-size: 75%;"
1065
| colspan="2" | '''Figure 9:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Axial displacement <math>u</math> at the upper and lower surfaces of the top layer (layer 3)
1066
|}
1067
1068
<div id='img-10a'></div>
1069
<div id='img-10b'></div>
1070
<div id='img-10c'></div>
1071
<div id='img-10'></div>
1072
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1073
|-
1074
|[[Image:Draft_Samper_624436055-test-Fig10a.png|420px|]]
1075
|[[Image:Draft_Samper_624436055-test-Fig10b.png|420px|]]
1076
|- style="text-align: center; font-size: 75%;"
1077
| (a) 
1078
| (b) 
1079
|-
1080
| colspan="2"|[[Image:Draft_Samper_624436055-test-Fig10c.png|420px|]]
1081
|- style="text-align: center; font-size: 75%;"
1082
|  colspan="2" | (c) 
1083
|- style="text-align: center; font-size: 75%;"
1084
| colspan="2" | '''Figure 10:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Thickness distribution of the axial displacement <math>u</math> at <math>x=L/4</math> (a), <math>x=L/2</math> (b) and <math>x=L</math> (c)
1085
|}
1086
1087
Figure [[#img-10|10]] shows the thickness distribution for the axial displacement at sections located at distances <math display="inline">\frac{L}{4},\frac{L}{2}</math> and <math display="inline">\frac{3L}{4}</math> from the clamped end. Results for the LRZ element (LRZ-25, LRZ-50 and LRZ-300) are in good agreement  with the reference solution. The TBT results have the standard linear distribution which is far from the correct zigzag results.
1088
1089
Figure [[#img-11|11]] shows the distribution along the beam length of the axial stress <math display="inline">\sigma _x</math> at the top and bottom surfaces of the beam cross section. Very good agreement between the reference PS solution and the LRZ-50 and LRZ-300 results is obtained. Results for the LRZ-25 mesh compare reasonably well with the PS solution except in the vicinity of the clamped edge. This error is corrected for the LRZ-50 and LRZ-300 meshes. The TBT results yield a linear distribution of the axial stress along the beam, as expected. This introduces large errors in the axial stress values in the vicinity of the clamped edge, as clearly shown in  Figure [[#img-11|11]].
1090
1091
<div id='img-11a'></div>
1092
<div id='img-11b'></div>
1093
<div id='img-11'></div>
1094
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1095
|-
1096
|[[Image:Draft_Samper_624436055-test-Fig11a.png|420px|]]
1097
|[[Image:Draft_Samper_624436055-test-Fig11b.png|420px|]]
1098
|- style="text-align: center; font-size: 75%;"
1099
| (a) 
1100
| (b) 
1101
|- style="text-align: center; font-size: 75%;"
1102
| colspan="2" | '''Figure 11:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Axial stress <math>\sigma _x</math> at upper (a) and lower (b) surfaces of the cross section along the beam length
1103
|}
1104
1105
Figures [[#img-12|12]] and [[#img-13|13]] show the thickness distribution for the axial stress <math display="inline">\sigma _x</math> at the clamped section and at the center of the beam. LRZ results  agree quite well with those of the reference solution. TBT results have an erroneous stress distribution for the top and bottom layers at the clamped end. These differences are less important for the central section.
1106
1107
Figure [[#img-14|14]] shows the thickness distribution for the transverse shear stress <math display="inline">\tau _{xz}</math> at different sections (<math display="inline">\frac{L}{20},\frac{L}{4},\frac{L}{2}</math> and <math display="inline">\frac{3L}{4}</math>). LRZ results provide an accurate estimate of the average transverse shear stress value for each layer. The distribution of <math display="inline">\tau _{xz}</math> across the thickness can be substantially improved by using the equilibrium equations for computing <math display="inline">\tau _{xz}</math> “a posteriori” as explained in Section 7.3.
1108
1109
TBT results are acceptable for the central layer and clearly overestimate the transverse shear stress in sections far from the clamped end.
1110
1111
LRZ and TBT results for the distribution of the (constant) tangential shear stress <math display="inline">\tau _{xz}</math> for each of the three layers along the beam length are shown in Figure [[#img-15|15]]. TBT results are clearly inaccurate (except for the values at the clamped edge).
1112
1113
<div id='img-12'></div>
1114
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1115
|-
1116
|[[Image:Draft_Samper_624436055-test-Fig12.png|600px|Non symmetric 3-layered cantilever thick beam under end point load (λ=5). Thickness distribution of the axial stress σₓ at x=0]]
1117
|- style="text-align: center; font-size: 75%;"
1118
| colspan="1" | '''Figure 12:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Thickness distribution of the axial stress <math>\sigma _x</math> at <math>x=0</math>
1119
|}
1120
1121
<div id='img-13'></div>
1122
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1123
|-
1124
|[[Image:Draft_Samper_624436055-test-Fig13.png|600px|Non symmetric 3-layered cantilever thick beam under end point load (λ=5). Thickness distribution of the axial stress σₓ at x=L/2]]
1125
|- style="text-align: center; font-size: 75%;"
1126
| colspan="1" | '''Figure 13:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Thickness distribution of the axial stress <math>\sigma _x</math> at <math>x=L/2</math>
1127
|}
1128
1129
<div id='img-14a'></div>
1130
<div id='img-14b'></div>
1131
<div id='img-14c'></div>
1132
<div id='img-14d'></div>
1133
<div id='img-14'></div>
1134
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1135
|-
1136
|[[Image:Draft_Samper_624436055-test-Fig14a.png|420px|]]
1137
|[[Image:Draft_Samper_624436055-test-Fig14b.png|420px|]]
1138
|- style="text-align: center; font-size: 75%;"
1139
| (a) 
1140
| (b) 
1141
|-
1142
|[[Image:Draft_Samper_624436055-test-Fig14c.png|420px|]]
1143
|[[Image:Draft_Samper_624436055-test-Fig14d.png|420px|]]
1144
|- style="text-align: center; font-size: 75%;"
1145
| (c) 
1146
| (d) 
1147
|- style="text-align: center; font-size: 75%;"
1148
| colspan="2" | '''Figure 14:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). Thickness distribution of transverse shear stress <math>\tau _{xz}</math> at <math>L/20</math> (a), <math>L/4</math> (b), <math>L/2</math> (c) and <math>3L/4</math> (d) 
1149
|}
1150
1151
<div id='img-15a'></div>
1152
<div id='img-15b'></div>
1153
<div id='img-15c'></div>
1154
<div id='img-15'></div>
1155
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1156
|-
1157
|[[Image:Draft_Samper_624436055-test-Fig15a.png|420px|]]
1158
|[[Image:Draft_Samper_624436055-test-Fig15b.png|420px|]]
1159
|- style="text-align: center; font-size: 75%;"
1160
| (a) 
1161
| (b) 
1162
|-
1163
| colspan="2"|[[Image:Draft_Samper_624436055-test-Fig15c.png|420px|]]
1164
|- style="text-align: center; font-size: 75%;"
1165
|  colspan="2" | (c) 
1166
|- style="text-align: center; font-size: 75%;"
1167
| colspan="2" | '''Figure 15:''' Non symmetric 3-layered cantilever thick beam under end point load (<math>\lambda =5</math>). LRZ and TBT results for the transverse shear stress <math>\tau _{xz}</math> along the beam. Layer 1 (a), layer 2 (b) and layer 3(c)
1168
|}
1169
1170
===7.2 Three-layered simple supported (SS) thick beams under uniform load===
1171
1172
The next example is the analysis of a three-layered simple supported thick beam under a uniformly distributed load of unit value (<math display="inline">q=1</math>). The material properties and the thickness for the three layers are shown in Table [[#table-6|6]]. ''The material has a non symmetric distribution'' with respect to the beam axis. An unusually low value for the shear modulus of the core layer has been taken, thus reproducing the effect of a damaged material in this zone. The span-to-thickness ratio is <math display="inline">\lambda =5</math>. Results obtained with the LRZ element  are once more compared with those obtained with a mesh of 300 2-noded TBT elements and with the mesh of 27000 4-noded plane stress (PS) rectangles shown in Figure [[#img-3|3]]. The PS solution has been obtained by fixing the vertical displacement of all nodes at the end sections and the horizontal displacement of the mid-line node at <math display="inline">x=0</math> and <math display="inline">x=L</math> to a zero value. This way of approximating a simple support condition leads to some discrepancies between the PS results and those obtained with beam theory.
1173
1174
No advantage of the symmetry of the problem for the discretization has been taken.
1175
1176
Figure [[#img-16|16]] shows the distribution of the vertical deflection for the different methods. The error in the “best” maximum central deflection value versus the “exact” PS solution is <math display="inline">\simeq </math> 12% The discrepancy is due to the difference in the way the simple support condition is modelled in beam and PS theories, as well as to the limitations of beam theory to model accurately very thick beams. TBT results are inaccurate, as expected.
1177
1178
Figure [[#img-17|17]] shows the distribution of the axial stress <math display="inline">\sigma _x</math> along the beam for the top surface of the second and third layer.
1179
1180
The accuracy of the LRZ results is remarkable with a maximum error of 10% despite of the modeling limitations mentioned above. TBT results are incorrect.
1181
1182
Figure [[#img-18|18]] shows the thickness distribution of the axial displacement at the left end section and at <math display="inline">x = \frac{L}{4}</math>. The LRZ element captures very well the zigzag shape of the axial displacement field even for a coarse mesh of 10 elements. The TBT element yields an unrealistic linear distribution.
1183
1184
Figures [[#img-19|19]] and [[#img-20|20]] show the thickness distribution of the axial stress and the transverse shear stress at the left end and mid sections. The accuracy of the LRZ results is again noticeable (even for the coarse 10 element mesh). The TBT element fails to capture the zigzag distribution of the axial stress (Figure [[#img-19|19]]) and gives a wrong value of almost zero shear stress at the core layer for the two sections chosen (Figure [[#img-20|20]]).
1185
1186
Figure [[#img-21|21]] shows  the distribution of the shear stress <math display="inline">\tau _{xz}</math> along the beam for each of the three layers obtained with the LRZ and TBT elements. TBT results are accurate for the first and third layer but are wrong for the core layer.
1187
1188
Figure [[#img-22|22]] shows a similar set of results for a moderately thick SS beam (<math display="inline">\lambda =10</math>) and the same material properties. Results shown are the distribution along the beam of the deflection and the axial stress at the top surface of layer 2. The accuracy of the LRZ element is again noticeable.
1189
1190
1191
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1192
|+ style="font-size: 75%;" |<span id='table-6'></span>Table. 6 Thickness and material properties for 3-layered non-symmetric simple supported (SS)  beam
1193
1194
|-
1195
| 
1196
| colspan='3' | '''Thickness and material properties'''
1197
1198
|- style="border-top: 2px solid;"
1199
| 
1200
| Layer 1 (bottom) 
1201
| Layer 2 (core) 
1202
| Layer 3 (top)
1203
1204
|- style="border-top: 2px solid;"
1205
| h [mm] 
1206
| 6.6666 
1207
| 6.6666 
1208
| 6.6666 
1209
|-
1210
| E [MPa] 
1211
| 2.19E<math display="inline">^{5}</math> 
1212
| 5.30E<math display="inline">^{5}</math> 
1213
| 7.30E<math display="inline">^{5}</math> 
1214
|- style="border-bottom: 2px solid;"
1215
| G [MPa] 
1216
| 8.76E<math display="inline">^{4}</math> 
1217
| 2.90E<math display="inline">^{2}</math> 
1218
| 2.92E<math display="inline">^{5}</math> 
1219
1220
|}
1221
1222
<div id='img-16'></div>
1223
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1224
|-
1225
|[[Image:Draft_Samper_624436055-test-Fig16.png|600px|Non symmetric 3-layered SS thick beam under uniformly distributed load \left(λ=5\right). Distribution of vertical deflection w along the beam length]]
1226
|- style="text-align: center; font-size: 75%;"
1227
| colspan="1" | '''Figure 16:''' Non symmetric 3-layered SS thick beam under uniformly distributed load <math>\left(\lambda =5\right)</math>. Distribution of vertical deflection <math>w</math> along the beam length
1228
|}
1229
1230
<div id='img-17a'></div>
1231
<div id='img-17b'></div>
1232
<div id='img-17'></div>
1233
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1234
|-
1235
|[[Image:Draft_Samper_624436055-test-Fig17a.png|420px|]]
1236
|[[Image:Draft_Samper_624436055-test-Fig17b.png|420px|]]
1237
|- style="text-align: center; font-size: 75%;"
1238
| (a) 
1239
| (b) 
1240
|- style="text-align: center; font-size: 75%;"
1241
| colspan="2" | '''Figure 17:''' Non symmetric 3-layered SS thick beam under uniformly distributed load <math>(\lambda =5)</math>. Distribution of axial stress <math>\sigma _{x}</math> at upper surface of layer 2 (a) and  layer 3 (b)
1242
|}
1243
1244
<div id='img-18a'></div>
1245
<div id='img-18b'></div>
1246
<div id='img-18'></div>
1247
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1248
|-
1249
|[[Image:Draft_Samper_624436055-test-Fig18a.png|420px|x=0]]
1250
|[[Image:Draft_Samper_624436055-test-Fig18b.png|420px|x=\fracL2]]
1251
|- style="text-align: center; font-size: 75%;"
1252
| (a) <math>x=0</math>
1253
| (b) <math>x=\frac{L}{2}</math>
1254
|- style="text-align: center; font-size: 75%;"
1255
| colspan="2" | '''Figure 18:''' Non symmetric 3-layered SS thick beam under uniformly distributed load <math>(\lambda =5)</math>.  Thickness distribution of axial displacement at <math>x=0</math> (a) and at <math>x=L/2</math> (b)
1256
|}
1257
1258
<div id='img-19a'></div>
1259
<div id='img-19b'></div>
1260
<div id='img-19'></div>
1261
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1262
|-
1263
|[[Image:Draft_Samper_624436055-test-Fig19a.png|420px|]]
1264
|[[Image:Draft_Samper_624436055-test-Fig19b.png|420px|]]
1265
|- style="text-align: center; font-size: 75%;"
1266
| (a) 
1267
| (b) 
1268
|- style="text-align: center; font-size: 75%;"
1269
| colspan="2" | '''Figure 19:''' Non symmetric 3-layered SS thick beam under uniformly distributed load <math>(\lambda =5)</math>.  Thickness distribution of axial stress <math>\sigma _{x}</math> at <math>x=0</math> (a) and at <math>x=L/2</math> (b)
1270
|}
1271
1272
<div id='img-20a'></div>
1273
<div id='img-20b'></div>
1274
<div id='img-20'></div>
1275
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1276
|-
1277
|[[Image:Draft_Samper_624436055-test-Fig20a.png|420px|]]
1278
|[[Image:Draft_Samper_624436055-test-Fig20b.png|420px|]]
1279
|- style="text-align: center; font-size: 75%;"
1280
| (a) 
1281
| (b) 
1282
|- style="text-align: center; font-size: 75%;"
1283
| colspan="2" | '''Figure 20:''' Non symmetric 3-layered SS thick beam under uniformly distributed load <math>(\lambda =5)</math>.  Thickness distribution of the  shear stress at <math>x=L/20</math> (a) and at <math>x=L/4</math> (b)
1284
|}
1285
1286
<div id='img-21a'></div>
1287
<div id='img-21b'></div>
1288
<div id='img-21c'></div>
1289
<div id='img-21'></div>
1290
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1291
|-
1292
|[[Image:Draft_Samper_624436055-test-Fig21a.png|420px|]]
1293
|[[Image:Draft_Samper_624436055-test-Fig21b.png|420px|]]
1294
|- style="text-align: center; font-size: 75%;"
1295
| (a) 
1296
| (b) 
1297
|-
1298
| colspan="2"|[[Image:Draft_Samper_624436055-test-Fig21c.png|420px|]]
1299
|- style="text-align: center; font-size: 75%;"
1300
|  colspan="2" | (c) 
1301
|- style="text-align: center; font-size: 75%;"
1302
| colspan="2" | '''Figure 21:''' Non symmetric 3-layered SS thick beam under uniformly distributed load <math>(\lambda =5)</math>. LRZ and TBT results for the distribution of <math>\tau _{xz}</math> along the beam for layer 1 (a), layer 2 (b) and layer 3 (c)
1303
|}
1304
1305
<div id='img-22a'></div>
1306
<div id='img-22b'></div>
1307
<div id='img-22'></div>
1308
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1309
|-
1310
|[[Image:Draft_Samper_624436055-test-Fig22a.png|420px|]]
1311
|[[Image:Draft_Samper_624436055-test-Fig22b.png|420px|]]
1312
|- style="text-align: center; font-size: 75%;"
1313
| (a) 
1314
| (b) 
1315
|- style="text-align: center; font-size: 75%;"
1316
| colspan="2" | '''Figure 22:''' Non symmetric 3-layered SS moderately thick beam under uniformly distributed load <math>(\lambda =10)</math>. Distribution along the beam length of the vertical deflection <math>w</math> (a) and the axial stress <math>\sigma _x</math> at the upper of layer 2 (b)
1317
|}
1318
1319
===7.3 Non-symmetric ten-layered clamped slender beam under uniformly distributed loading===
1320
1321
We present results for  a ten-layered clamped slender rectangular beam (<math display="inline">L=100</math> mm, <math display="inline">h=5</math> mm, <math display="inline">b=1</math> mm, <math display="inline">\lambda =20</math>) under uniformly distributed loading (<math display="inline">q =1</math> KN/mm). The composite material has the non-symmetric distribution across the thickness shown in Table 7.
1322
1323
1324
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1325
|+ style="font-size: 75%;" |<span id='table-7'></span>Table. 7 10-layered clamped slender rectangular beam under uniformly distributed loading. (a) Thickness and material number for each of the 10 layers. (b) Properties of each material
1326
|-
1327
| colspan='3' | (a)
1328
1329
|- style="border-top: 2px solid;"
1330
| Layer 
1331
| <math>h_i</math>
1332
| Material
1333
1334
|- style="border-top: 2px solid;"
1335
| 1 
1336
| 0.5 
1337
| IV 
1338
|-
1339
| 2 
1340
| 0.6 
1341
| I 
1342
|-
1343
| 3 
1344
| 0.5 
1345
| V 
1346
|-
1347
| 4 
1348
| 0.4 
1349
| III 
1350
|-
1351
| 5 
1352
| 0.7 
1353
| IV 
1354
|-
1355
| 6 
1356
| 0.1 
1357
| III 
1358
|-
1359
| 7 
1360
| 0.4 
1361
| II 
1362
|-
1363
| 8 
1364
| 0.5 
1365
| V 
1366
|-
1367
| 9 
1368
| 0.3 
1369
| I 
1370
|- style="border-bottom: 2px solid;"
1371
| 10 
1372
| 1 
1373
| II 
1374
1375
|}
1376
1377
1378
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1379
|-
1380
| colspan='3' | (b)
1381
1382
|- style="border-top: 2px solid;"
1383
| Material 
1384
| E [MPa] 
1385
| G [MPa] 
1386
1387
|- style="border-top: 2px solid;"
1388
| I 
1389
| 2.19e5 
1390
| 0.876e5 
1391
|-
1392
| II 
1393
| 7.3e5 
1394
| 2.92e5 
1395
|-
1396
| III 
1397
| 0.0073e5 
1398
| 0.0029e5 
1399
|-
1400
| IV 
1401
| 5.3e5 
1402
| 2.12e5 
1403
|- style="border-bottom: 2px solid;"
1404
| V 
1405
| 0.82e5 
1406
| 0.328e5 
1407
1408
|}
1409
1410
Figure [[#img-23|23]] shows results for the deflection along the beam for LRZ meshes with 10 and 300 elements (LRZ-10 and LRZ-300). Results obtained with a mesh of 27.000 4-noded plane stress quadrilaterals  and with a mesh of 300 TBT elements are also shown for comparison. Note the accuracy of the coarse LRZ-10 mesh and the erroneous results of the TBT solution.
1411
1412
Figure [[#img-24|24]] shows the thickness distribution of the axial displacement and the axial stress (<math display="inline">\sigma _x</math>) for the section at <math display="inline">x=\frac{L}{4}</math>. The accuracy of the LRZ results is once more remarkable.
1413
1414
<div id='img-23'></div>
1415
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1416
|-
1417
|[[Image:Draft_Samper_624436055-test-Fig23.png|600px|10-layered clamped slender beam under uniform loading. Distribution of the deflection along the beam]]
1418
|- style="text-align: center; font-size: 75%;"
1419
| colspan="1" | '''Figure 23:''' 10-layered clamped slender beam under uniform loading. Distribution of the deflection along the beam
1420
|}
1421
1422
<div id='img-24a'></div>
1423
<div id='img-24b'></div>
1424
<div id='img-24'></div>
1425
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1426
|-
1427
|[[Image:Draft_Samper_624436055-test-Fig24a.png|420px|]]
1428
|[[Image:Draft_Samper_624436055-test-Fig24b.png|420px|]]
1429
|- style="text-align: center; font-size: 75%;"
1430
| (a) 
1431
| (b) 
1432
|- style="text-align: center; font-size: 75%;"
1433
| colspan="2" | '''Figure 24:''' 10-layered clamped slender beam under uniform loading. Thickness distribution of axial displacement (a) and axial stress <math>\sigma _x</math> (b) for <math>x=\frac{L}{4}</math>
1434
|}
1435
1436
Figure [[#img-25|25]] shows the thickness distribution of the transverse shear stress at <math display="inline">x=\frac{L}{4}</math>. Results in Figure [[#img-25|25]]a show the values directly obtained with the LRZ-10 and LRZ-300 meshes. These results are clearly better than those obtained with the TBT element but only coincide in an average sense with the plane stress FEM solution.
1437
1438
<div id='img-25a'></div>
1439
<div id='img-25b'></div>
1440
<div id='img-25'></div>
1441
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1442
|-
1443
|[[Image:Draft_Samper_624436055-test-Fig25a.png|420px|]]
1444
|[[Image:Draft_Samper_624436055-test-Fig25b.png|420px|]]
1445
|- style="text-align: center; font-size: 75%;"
1446
| (a) 
1447
| (b) 
1448
|- style="text-align: center; font-size: 75%;"
1449
| colspan="2" | '''Figure 25:''' 10-layered clamped slender beam under uniform loading. Thickness distribution of <math>\tau _{xz}</math> at <math>x=\frac{L}{4}</math>. (a) Comparison of LRZ-10 and LRZ-300  results with plane stress (PS) and TBT solutions. (b) PS solution and LRZ-10-<math>N_z</math> and LRZ-300-<math>N_z</math> results for <math>\tau _{xz}</math> obtained by thickness integration of the equilibrium equation  using the LRZ-10 and LRZ-300 results (Eq.(32))
1450
|}
1451
1452
LRZ results for the thickness distribution of <math display="inline">\tau _{xz}</math> can be much improved by computing <math display="inline">\tau _{xz}</math> “a posteriori” from the axial stress field using the equilibrium equation
1453
1454
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1455
|-
1456
| 
1457
{| style="text-align: left; margin:auto;width: 100%;" 
1458
|-
1459
| style="text-align: center;" | <math>\frac{\partial \sigma _x}{\partial x}+ \frac{\partial \tau _{xz}}{\partial z}=0  </math>
1460
|}
1461
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
1462
|}
1463
1464
The transverse shear stress at a point across the thickness with coordinate <math display="inline">z</math> is computed by integrating Eq.(31) as
1465
1466
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1467
|-
1468
| 
1469
{| style="text-align: left; margin:auto;width: 100%;" 
1470
|-
1471
| style="text-align: center;" | <math>\tau _{xz} (z) = - \int ^z_{-\frac{h}{2}} \frac{\partial \sigma _x}{\partial x} dz = - \frac{\partial N_z}{\partial x} </math>
1472
|}
1473
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
1474
|}
1475
1476
where
1477
1478
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1479
|-
1480
| 
1481
{| style="text-align: left; margin:auto;width: 100%;" 
1482
|-
1483
| style="text-align: center;" | <math>N_z = \int ^z_{-\frac{h}{2}} \sigma _x dz </math>
1484
|}
1485
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
1486
|}
1487
1488
is the axial force (per  unit width) resulting from the thickness integration of <math display="inline">\sigma _x</math> between the coordinates <math display="inline">-\frac{h}{2}</math> and <math display="inline">z</math>.
1489
1490
The space derivative of <math display="inline">N_z</math> in Eq.(32) is computed ''at a node'' <math display="inline">i</math> as
1491
1492
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1493
|-
1494
| 
1495
{| style="text-align: left; margin:auto;width: 100%;" 
1496
|-
1497
| style="text-align: center;" | <math>\frac{\partial N_z}{\partial x} = \frac{2}{l^e + l^{e-1}} (N_z^e- N_z^{e-1}) </math>
1498
|}
1499
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
1500
|}
1501
1502
where (<math display="inline">l^e,N_z^e</math>) and (<math display="inline">l^{e-1}, N_z^{e-1}</math>) are the element length and the value of <math display="inline">N_z</math> at elements <math display="inline">e</math> and <math display="inline">e-1</math> adjacent to node <math display="inline">i</math>, respectively. A value of <math display="inline">\tau _{xz} (-\frac{h}{2}) =0</math> is taken. It is remarkable that the method yields automatically <math display="inline">\tau _{xz} (\frac{h}{2})\simeq 0</math>.
1503
1504
Results for <math display="inline">\tau _{xz}</math> obtained with this procedure are termed LRZ-10-<math display="inline">N_z</math> and LRZ-300-<math display="inline">N_z</math> in Figure [[#img-25|25]]. We note the accuracy of the “recovered” thickness distribution for <math display="inline">\tau _{xz}</math>, even for the coarse mesh of 10 LRZ elements.
1505
1506
===7.4 Modeling of delamination with the LRZ element===
1507
1508
Prediction of delamination in composite laminated beams is a challenge for all beam models. A method for predicting delamination in beams using a Hermitian zigzag theory was presented in <span id='citeF-26'></span><span id='citeF-27'></span>[[#cite-26|[26,27]]]. A sub-laminate approach is used for which the number of kinematic unknowns  depends of the number of physical layers. This increases the number of variables but it yields the correct an accurate transverse shear stress distribution without integrating the equilibrium equations.
1509
1510
Delamination effects in composite laminated beams can be effectively reproduced with the LRZ element ''without introducing additional kinematic variables''. The delamination model simply implies ''introducing a very thin “interface layer” between adjacent material layers'' in the actual composite laminated section. Delamination is produced when the material properties of the interface layer are drastically reduced to almost a zero value in comparison with those of the adjacent layers due to interlamina failure.  This simple delamination model allows the LRZ element to take into account the reduction of the overall beam stiffness due to the failure of the interface layer leading to an increase in the deflection and rotation field. Moreover, the LRZ element can also accurately represent the jump in the axial displacement field across the interface layer and the change in the axial and tangential stress distributions over the beam sections as delamination progresses.
1511
1512
Figures [[#img-26|26]]&#8211;[[#img-30|30]] show an example of the capabilities of the LRZ beam element to model delamination. The problem represents the analysis of a cantilever thick rectangular beam (<math display="inline">\lambda = 5</math>) under an end point load. The beam section has three layers of composite material with properties  shown in Table [[#table-8|8]]. Delamination between the upper and core layers has been modelled by introducing a very thin interface layer (<math display="inline">h=0.01</math> mm) between these two layers (Figure [[#img-26|26]]). The initial properties of the interface layer coincide with those of the upper layer. Next, the  shear modulus value for the interface layer has been progressively reduced up to 11 orders of magnitude from <math display="inline">G2= 8.76 \times 10^4</math> MPa (Model 1) to <math display="inline">G2 =8.76 \times 10^{-7}</math> MPa (Model 12)  (Table [[#table-9|9]]).
1513
1514
<div id='img-26'></div>
1515
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1516
|-
1517
|[[Image:Draft_Samper_624436055-test-Fig26.png|600px|Modeling of interface layer for delamination study in 3-layered thick cantilever beam (λ= 5) under end point load]]
1518
|- style="text-align: center; font-size: 75%;"
1519
| colspan="1" | '''Figure 26:''' Modeling of interface layer for delamination study in 3-layered thick cantilever beam (<math>\lambda = 5</math>) under end point load
1520
|}
1521
1522
1523
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1524
|+ style="font-size: 75%;" |<span id='table-8'></span>Table. 8 Thickness and layer properties for delamination study in a 3-layered cantilever beam under end point load. Layer 2 is the interface layer. G2 values are given in Table 9
1525
|- style="border-top: 2px solid;"
1526
| style="text-align: left;" |  
1527
| colspan='4' | Composite material
1528
|- style="border-top: 2px solid;"
1529
| style="text-align: left;" |  
1530
| Layer 1
1531
| Layer 2
1532
| Layer 3
1533
| Layer 4
1534
1535
|- style="border-top: 2px solid;"
1536
| style="text-align: left;" | h [mm]
1537
| 2
1538
| 0.01
1539
| 16 
1540
| 2
1541
|-
1542
| style="text-align: left;" | E [MPa]
1543
| 2.19E5
1544
| 2.19E5
1545
| 0.0073E5
1546
| 7.30E5
1547
|- style="border-bottom: 2px solid;"
1548
| style="text-align: left;" | G [MPa] 
1549
| 0.876E5 
1550
| G2 
1551
| 0.0029E5 
1552
| 2.92E5
1553
1554
|}
1555
1556
1557
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1558
|+ style="font-size: 75%;" |<span id='table-9'></span>Table. 9 Shear modulus values for the interface layer for delamination study in a 3-layered cantilever beam. Values of G2 in MPa
1559
|- style="border-top: 2px solid;"
1560
|  Model 
1561
| G2  
1562
| Model 
1563
| G2  
1564
| Model
1565
| G2 
1566
|- style="border-top: 2px solid;"
1567
|  1
1568
| 8.76E+004
1569
| 5
1570
| 8.76E+000 
1571
| 9 
1572
| 8.76E-004
1573
|-
1574
| 2
1575
| 8.76E+003
1576
| 6
1577
| 8.76E-001 
1578
| 10 
1579
| 8.76E-005
1580
|-
1581
| 3
1582
| 8.76E+002
1583
| 7
1584
| 8.76E-002
1585
| 11 
1586
| 8.76E-006
1587
|- style="border-bottom: 2px solid;"
1588
| 4
1589
| 8.76E+001
1590
| 8
1591
| 8.76E-003
1592
| 12 
1593
| 8.76E-007
1594
1595
|}
1596
1597
We note that the reduction of the shear  modulus has been applied over the whole beam length in this case. However it can applied in selected beam regions as appropriate.
1598
1599
Figure [[#img-27|27]] shows results for the end deflection  in terms of the shear modulus value of the interface layer for the LRZ-100 mesh. Note that the deflection increases one order of magnitude versus the non-delaminated case. It is also interesting that the end deflection does not change after the shear modulus of the interface layer is reduced beyond eight orders of magnitude (results for Model 9 in Figure [[#img-27|27]]). Results agree reasonably well (error <math display="inline">\simeq 10%</math>) with those obtained with the plane stress model of Figure 3 introducing a similar reduction  in the shear modulus of an  ''ad hoc'' interface layer.
1600
1601
<div id='img-27'></div>
1602
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1603
|-
1604
|[[Image:Draft_Samper_624436055-test-Fig27.png|600px|Delamination study in 3-layered cantilever beam under end point load. Evolution of end deflection with the shear modulus value for the interface layer  LRZ-100 results and PS solution]]
1605
|- style="text-align: center; font-size: 75%;"
1606
| colspan="1" | '''Figure 27:''' Delamination study in 3-layered cantilever beam under end point load. Evolution of end deflection with the shear modulus value for the interface layer  LRZ-100 results and PS solution
1607
|}
1608
1609
Figure [[#img-28|28]] shows the thickness distribution for the axial displacement at the mid section for four decreasing values of the shear modulus at the interface layer: <math display="inline">G2=8.76, 8.76 \times 10^{-1}, 8.76\times 10^{-3}</math> and <math display="inline">8.76 \times 10^{-6}</math> MPa. The  jump of the axial displacement across the thickness at the interface layer during delamination is well captured. We again note that the displacement jump at the interface layer remains stationary after a reduction of the material properties in that layer of six orders of magnitude. Results agree well with the plane stress solution also shown in the figure.
1610
1611
<div id='img-28'></div>
1612
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1613
|-
1614
|[[Image:Draft_Samper_624436055-test-Fig28a.png|420px|]]
1615
|[[Image:Draft_Samper_624436055-test-Fig28b.png|420px|]]
1616
|-
1617
|[[Image:Draft_Samper_624436055-test-Fig28c.png|420px|]]
1618
|[[Image:Draft_Samper_624436055-test-Fig28d.png|420px|Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of axial displacement at x=\fracL2 for four decreasing values of the shear modulus at the interface layer  (Models 5, 6, 8 and 11, Table 9)]]
1619
|- style="text-align: center; font-size: 75%;"
1620
| colspan="2" | '''Figure 28:''' Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of axial displacement at <math>x=\frac{L}{2}</math> for four decreasing values of the shear modulus at the interface layer  (Models 5, 6, 8 and 11, Table 9)
1621
|}
1622
1623
Figure [[#img-29|29]]  shows  the thickness distribution of the axial stress (<math display="inline">\sigma _x</math>) for  the same four decreasing values of the shear modulus at the interface layer. The effect of delamination in the stress distribution is clearly visible. Once again the LRZ-100 results agree well  with the plane stress solution.
1624
1625
<div id='img-29'></div>
1626
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1627
|-
1628
|[[Image:Draft_Samper_624436055-test-Fig29a.png|420px|]]
1629
|[[Image:Draft_Samper_624436055-test-Fig29b.png|420px|]]
1630
|-
1631
|[[Image:Draft_Samper_624436055-test-Fig29c.png|420px|]]
1632
|[[Image:Draft_Samper_624436055-test-Fig29d.png|420px|Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of σₓ at x=\fracL2 for four decreasing values of the shear modulus at the interface layer  (Models 5, 6, 8 and 11, Table 9)]]
1633
|- style="text-align: center; font-size: 75%;"
1634
| colspan="2" | '''Figure 29:''' Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of <math>\sigma _x</math> at <math>x=\frac{L}{2}</math> for four decreasing values of the shear modulus at the interface layer  (Models 5, 6, 8 and 11, Table 9)
1635
|}
1636
1637
Figure [[#img-30|30]] finally shows the thickness distribution for the transverse shear stress at <math display="inline">x=\frac{L}{2}</math> for the same four values of the shear modulus at the interface layer. The three graphs show  the PS results,  the LRZ-100 results and the solution obtained by integrating the equilibrium equation (via Eqs.(31)&#8211;(34)) using the LRZ-100 results. Note the accuracy of the later solution versus the standard LRZ-100 results as  delamination develops and the transverse shear stress progressively vanishes at the interface layer.
1638
1639
Similar good results for predicting the delamination and the thickness distribution of the axial and transverse shear stresses are obtained over the entire beam length.
1640
1641
The example shows clearly the capability of the LRZ element to model a complex phenomenon such as delamination in composite laminated beams ''without introducing additional kinematic variables''. More evidence of the good behaviour of the LRZ beam element for predicting delamination in beams are reported in <span id='citeF-29'></span>[[#cite-29|[29]]].
1642
1643
<div id='img-30'></div>
1644
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1645
|-
1646
|[[Image:Draft_Samper_624436055-test-Fig30a.png|420px|]]
1647
|[[Image:Draft_Samper_624436055-test-Fig30b.png|420px|]]
1648
|-
1649
|[[Image:Draft_Samper_624436055-test-Fig30c.png|420px|]]
1650
|[[Image:Draft_Samper_624436055-test-Fig30d.png|420px|Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of τ<sub>xz</sub> at x=\fracL2 for four values of G at the interface layer (Models 5, 6, 8 and 11, Table 9). LRZ-100 results, plane stress (PS) solution and LRZ-100-Sx results obtained by integrating the equilibrium equation (Eq.(32)) using the LRZ-100 results]]
1651
|- style="text-align: center; font-size: 75%;"
1652
| colspan="2" | '''Figure 30:''' Delamination study in 3-layered cantilever beam under end point load. Thickness distribution of <math>\tau _{xz}</math> at <math>x=\frac{L}{2}</math> for four values of G at the interface layer (Models 5, 6, 8 and 11, Table 9). LRZ-100 results, plane stress (PS) solution and LRZ-100-Sx results obtained by integrating the equilibrium equation (Eq.(32)) using the LRZ-100 results
1653
|}
1654
1655
==8 CONCLUSIONS==
1656
1657
We have presented a simple and accurate 2-noded beam element based on the combination of the refined zigzag beam theory proposed by Tessler ''et al.'' <span id='citeF-22'></span>[[#cite-22|[22]]]. The element has four degrees of freedom per node (the axial displacement, the deflection, the rotation and the amplitude of the zigzag function). A standard <math display="inline">C^\circ </math> linear interpolation  is used for all variables. The resulting LRZ beam element is shear locking-free and has shown an excellent behaviour for analysis of thick and thin composite laminated beams with clamped and simple supported conditions. Numerical results  agree in practically all cases with  those obtained with a two-dimensional plane-stress FEM using a far larger number of degrees of freedom. It is remarkable that the zigzag distribution of the axial displacement and the axial stress across the thickness, typical of composite laminated beams, is very accurately captured with the basic approximation chosen. The possibilities of the new LRZ beam element for predicting delamination effects has been demonstrated in a simple but representative example of application.
1658
1659
==ACKNOWLEDGEMENTS==
1660
1661
This research was partially supported by project SEDUREC of the Consolider Programme of the Ministerio de Educación y Ciencia of Spain.
1662
1663
==REFERENCES==
1664
1665
<div id="cite-1"></div>
1666
[[#citeF-1|[1]]]  Timoshenko S.P. and Woinowsky-Krieger S., ''Theory of Plates and Shells''. McGraw-Hill, New York, 3rd Edition, 1959.
1667
1668
<div id="cite-2"></div>
1669
[[#citeF-2|[2]]]  Timoshenko S. On the correction for shear of differential equations for transverse vibrations of prismatic bars. Philosophical Magazine Series, Vol. 41, pp. 744-746, 1921.
1670
1671
<div id="cite-3"></div>
1672
[[#citeF-3|[3]]]  Liu D. and Li X. An overall view of laminate theories based on displacement hypothesis. Journal of Composite Materials, Vol. 30(14), pp. 1539-1561, 1996.
1673
1674
<div id="cite-4"></div>
1675
[[#citeF-4|[4]]]  Reddy, J.N., Mechanics of laminated composite plates and shells. Theory and analysis. 2nd Edition, CRC Press, Boca Raton, 2004.
1676
1677
<div id="cite-5"></div>
1678
[[#citeF-5|[5]]]  Owen D.R.J.  and  Li Z.H., A refined analysis of laminated plates by finite element displacement methods. Part I. Fundamentals and static analysis. II Vibration and stability.  Comp. Struct, Vol. 26, pp. 907&#8211;923, 1987.
1679
1680
<div id="cite-6"></div>
1681
[[#citeF-6|[6]]]  Botello, S. Oñate, E. and Canet, J.M., A layer-wise triangle for analysis of laminated composite plates and shells. Computers and Structures, Vol. 70, 635&#8211;646, 1999.
1682
1683
<div id="cite-7"></div>
1684
[[#citeF-7|[7]]]  Di Sciuva M. Bending, vibration and buckling of simply supported thick multilayered orthotropic plates: an evaluation of a new displacement model. Journal of Sound and Vibration, Vol. 105, pp. 425-442, 1986.
1685
1686
<div id="cite-8"></div>
1687
[[#citeF-8|[8]]]  Murakami H. Laminated composite plate theory with improved in-plane responses. ASME Journal of Applied Mechanics, Vol. 53, pp. 661-666, 1986.
1688
1689
<div id="cite-9"></div>
1690
[[#citeF-9|[9]]]  Aitharaju V.R. and Averill R.C., An assessment of zig-zag kinematic displacement models for the anlaysis of laminated composites. Mech. Composite Mat. and Struct., 6, 1&#8211;26, 1999a.
1691
1692
<div id="cite-10"></div>
1693
[[#citeF-10|[10]]]  Aitharaju V.R. and Averill R.C., C<math display="inline">^\circ </math> zig-zag finite element for analysis of laminated composite beams. J. Engineering Mechanics, ASCE, 323&#8211;330, 1999b.
1694
1695
<div id="cite-11"></div>
1696
[[#citeF-11|[11]]]  Li X. and Liu D. An interlaminar shear stress continuity theory for both thin and thick composite laminates. J. Appl. Mech., 59, 502&#8211;509, 1992.
1697
1698
<div id="cite-12"></div>
1699
[[#citeF-12|[12]]]  Di Sciuva M. An improved shear-deformation theory for moderately thick multilayered anisotropic shells and plates. J. Appl. Mech. 54, 589&#8211;594, 1987.
1700
1701
<div id="cite-13"></div>
1702
[[#citeF-13|[13]]]  Toledano A and Murakami H. A higher-order laminate plate theory with improved in-plane response. Int. J. Solids and Struct., 23, 111&#8211;131, 1987b.
1703
1704
<div id="cite-14"></div>
1705
[[#citeF-14|[14]]]  Cho M. and Parmerter R.R. Efficient higher order composite plate theory for general laminations configuration. AIAA J., 31, 1299&#8211;1306, 1993.
1706
1707
<div id="cite-15"></div>
1708
[[#citeF-15|[15]]]  Kapuria, S., Dumir, P.C., Ahmed, A. and Alam, N. Finite element model of efficient zigzag theory for static analysis of hybrid piezoelectric beams, Computational Mechanics,  Vol. 34 (6), pp. 475&#8211;483, 2004.
1709
1710
<div id="cite-16"></div>
1711
[[#citeF-16|[16]]]  Averill R.C. Static and dynamic response of moderately thick laminated beams with damage. Composites Engineering, Vol. 4(4), pp. 381-395, 1994.
1712
1713
<div id="cite-17"></div>
1714
[[#citeF-17|[17]]]  Averill R.C. and Yuen Cheong Yip. Development of simple, robust finite elements based on refined theories for thick laminated beams. Computers & Structures, Vol. 59(3), pp. 529-546, 1996.
1715
1716
<div id="cite-18"></div>
1717
[[#citeF-18|[18]]]  Alam N.M. and Upadhyay N.Kr. Finite element analysis of laminated composite beams for zigzag theory using MATLAB. Int. J. of Mechanics and Solids, Vol. 5(1), pp. 1&#8211;14, 2010.
1718
1719
<div id="cite-19"></div>
1720
[[#citeF-19|[19]]]  Savoia M. On the accuracy of one-dimensional models for multilayered composite beams. Int. J. Solids Struct.,  Vol. 33, pp. 521–-44, 1996.
1721
1722
<div id="cite-20"></div>
1723
[[#citeF-20|[20]]]   Kapuria S, Dumir P.C. and Jain N.K. Assessment of zigzag theory for static loading, buckling, free and forced response of composite and sandwich beams. Composite Structures,  Vol. 64, pp. 317–-27, 2004.
1724
1725
<div id="cite-21"></div>
1726
[[#citeF-21|[21]]]  Carrera, E. Historical review of zigzag theories for multilayered plate and shell. Applied Mechanics Review, 56(3), 287&#8211;308, 2003.
1727
1728
<div id="cite-22"></div>
1729
[[#citeF-22|[22]]]  Tessler, A., Di Sciuva M. and Gherlone, M. A refined zigzag beam theory for composite and sandwich beams. J. of Composite Materials, Vol. 43, 1051&#8211;1081, 2009.
1730
1731
<div id="cite-23"></div>
1732
[[#citeF-23|[23]]]  Tessler, A., Di Sciuva M. and Gherlone, M. A consistent refinement of first-order shear-deformation theory for laminated composite and sandwich plates using improved zigzag kinematics. J. of Mechanics of Materials and Structures, 5(2), 341&#8211;367, 2010.
1733
1734
<div id="cite-24"></div>
1735
[[#citeF-24|[24]]]  Gherlone, M., Tessler, A. and Di Sciuva M. C<math display="inline">^\circ </math> beam element based on the refined zigzag theory for multilayered composite and sandwich laminates. Composite Structures. Accepted for publication, May 2011.
1736
1737
<div id="cite-25"></div>
1738
[[#citeF-25|[25]]]  Oñate, E., Eijo, A. and Oller, S. Two-noded beam element for composite and sandwich beams using Timoshenko theory and refined zigzag kinematics. Publication CIMNE PI346, CIMNE, Barcelona, October 2010.
1739
1740
<div id="cite-26"></div>
1741
[[#citeF-26|[26]]]  Di Sciuva M. and Gherlone, M. A global/local third-order Hermitian displacement field with damaged interfaces and transverse extensibility: FEM formulation. Composite Structures, 59(4), 433-444, 2003.
1742
1743
<div id="cite-27"></div>
1744
[[#citeF-27|[27]]]  Di Sciuva M. and Gherlone, M. Quasi-3D static and dynamic analysis of undamaged and damaged sandwich beams. Journal of Sandwich Structures & Materials, 7(1), 31&#8211;52, 2005.
1745
1746
<div id="cite-28"></div>
1747
[[#citeF-28|[28]]]  Oñate, E. ''Structural Analysis with the Finite Element Method''. Vol. 2 Beams, Plates and Shells. CIMNE-Springer, 2011.
1748
1749
<div id="cite-29"></div>
1750
[[#citeF-29|[29]]]  Oñate, E., Eijo, A. and Oller, S. Modeling of delamination in composite laminated beams using a two-noded beam element based in refined zigzag theory. Publication CIMNE, PI367, November 2011.
1751

Return to Onate et al 2012a.

Back to Top