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
Published in ''Computer Methods in Applied Mechanics and Engineering, vol. 194, 2406-2443, 2005''
2
3
4
5
== SUMMARY ==
6
7
A family of rotation-free three node triangular shell elements is presented. The simplest element of the family is based on an assumed constant curvature field expressed in terms of the nodal deflections of a patch of four elements and a constant membrane field computed from the standard linear interpolation of the displacements within each triangle. An enhanced version of the element is obtained by using a quadratic interpolation of the geometry in terms of the six patch nodes. This allows to compute an assumed linear membrane strain field which improves the in-plane behaviour of the original element. A simple and economic version of the element using a single integration point is presented. The efficiency of the different rotation-free shell triangles is demonstrated in many examples of application including linear and non linear analysis of shells under static and dynamic loads, the inflation and de-inflation of membranes and a sheet stamping problem.
8
9
==1 INTRODUCTION==
10
11
Triangular shell elements are very useful for the solution of large scale shell problems such as those occurring in many practical engineering situations. Typical examples are the analysis of shell roofs under static and dynamic loads, sheet stamping processes, vehicle dynamics and crash-worthiness situations. Many of these problems involve high geometrical and material non linearities and time changing frictional contact conditions. These difficulties are usually increased by the need of discretizing complex geometrical shapes. Here the use of shell triangles and non-structured meshes becomes a critical necessity. Despite recent advances in the field <span id='citeF-1'></span>[[#cite-1|[1]]]&#8211;<span id='citeF-6'></span>[[#cite-6|[6]]] there are not so many simple shell triangles which are capable of accurately modelling the deformation of a shell structure under arbitrary loading conditions.
12
13
A promising line to derive simple shell triangles is to use the nodal displacements as the only unknown for describing the shell kinematics. This idea goes back to the original attempts to solve thin plate bending problems using finite difference schemes with the deflection as the only nodal variable <span id='citeF-7'></span>[[#cite-7|[7]]]&#8211;<span id='citeF-9'></span>[[#cite-9|[9]]].
14
15
In past years some authors have derived a number of thin plate and shell triangular elements free of rotational degrees of freedom (d.o.f.) based on Kirchhoff's theory [10]&#8211;<span id='citeF-27'></span>[[#cite-27|[27]]]. In essence all methods attempt to express the curvatures field over an element in terms of the displacements of a collection of nodes belonging to a patch of adjacent elements. Oñate and Cervera [14] proposed a general procedure of this kind combining finite element and finite volume concepts for deriving thin plate triangles and quadrilaterals with the deflection as the only nodal variable and presented a simple and competitive rotation-free three d.o.f. triangular element termed BPT (for Basic Plate Triangle). These ideas were extended and formalized in <span id='citeF-29'></span>[[#cite-29|29]] to derive a number of rotation-free thin plate and shell triangles. The basic ingredients of the method are a mixed Hu-Washizu formulation, a standard discretization into three node triangles, a linear finite element interpolation of the displacement field within each triangle and a finite volume type approach for computing constant curvature and bending moment fields within appropriate non-overlapping control domains. The so called cell-centered and cell-vertex triangular domains yield different families of rotation-free plate and shell triangles. Both the BPT plate element and its extension to shell analysis (termed BST for Basic Shell Triangle) can be derived from the cell-centered formulation. Here the control domain is an individual triangle. The constant curvatures field within a triangle is computed in terms of the displacements of the six nodes belonging to the four elements patch formed by the chosen triangle and the three adjacent triangles. The cell-vertex approach yields a different family of rotation-free plate and shell triangles. Details of the derivation of both rotation-free triangular shell element families can be found in <span id='citeF-29'></span>[[#cite-29|[21]]].
16
17
An extension of the BST element to the non linear analysis of shells was implemented in an explicit dynamic code by Oñate ''et al.'' [26] using an updated lagrangian formulation and a hypo-elastic constitutive model. Excellent numerical results were obtained for non linear dynamics of shells involving frictional contact situations and sheet stamping problems [18,19,20,26].
18
19
A large strain formulation for the BST element using a total lagrangian description was presented by Flores and Oñate [24]. A recent extension of this formulation is based on a quadratic interpolation of the geometry of the patch formed by the BST element and the three adjacent triangles [27]. This yields a linear displacement gradient field over the element from which linear membrane strains and  constant curvatures  can be computed within the BST element.
20
21
In this paper the formulation of the BST element is revisited using an assumed strain approach. The content of the paper is the following. First some basic concepts of the formulation of the original BST element using an assumed constant curvature field are given. Next, the basic equations of the non linear thin shell theory chosen based on a total lagrangian description are presented. Then the non linear formulation of the BST element is presented. This is based on an assumed constant membrane field derived from the linear displacement interpolation and an assumed constant curvature field expressed in terms of the displacements of the nodes of the four element patch using a finite volume type approach. An enhanced version of the BST element is derived using an assumed linear field for the membrane strains and an assumed constant curvature field. Both assumed fields are obtained from the quadratic interpolation of the patch geometry following the ideas presented in [27]. Details of the derivation of the tangent stiffness matrix needed  for a quasi-static implicit solution are given for both the BST and EBST elements. An efficient version of the  EBST element using one single quadrature point for integration of the tangent matrix is  presented. An explicit  scheme adequate for dynamic analysis is  briefly described.
22
23
The efficiency and accuracy of the standard and enhanced versions of the BST element is validated in a number of examples of application including linear and non linear analysis of shells under static and dynamic loads, the inflation and de-inflation of membranes and a sheet stamping problem.
24
25
==2 FORMULATION OF THE BASIC PLATE TRIANGULAR USING AN ASSUMED CONSTANT CURVATURE FIELD==
26
27
Let us consider a patch of four plate three node triangles (Figure [[#img-1|1]]). The nodes 1, 2, and 3 in the main central triangle (M) are marked with circles while the external nodes in the patch (nodes 4, 5 and 6) are marked with squares. Mid side points in the central triangle are also marked with smaller squares. Kirchhoff's thin plate theory will be assumed to hold. The deflection is linearly interpolated within each three node triangle in the standard finite element manner as
28
29
{| class="formulaSCP" style="width: 100%; text-align: left;" 
30
|-
31
| 
32
{| style="text-align: left; margin:auto;" 
33
|-
34
| style="text-align: center;" | <math>w=\sum _{i=1}^{3}L_{i}^{e}w_{i}^{e} </math>
35
|}
36
| style="width: 5px;text-align: right;" | (1)
37
|}
38
39
where <math display="inline">L_{i}^{e}</math> are the linear shape functions of the three node triangle, <math display="inline">w_{i}^{e}</math> are nodal deflections and superindex <math display="inline">e</math> denotes element values.
40
41
<div id='img-1'></div>
42
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
43
|-
44
|[[Image:draft_Samper_226033773-fig1.png|400px|Patch of three node triangular elements including the central triangle (M) and three adjacent triangles (1, 2 and 3)]]
45
|- style="text-align: center; font-size: 75%;"
46
| colspan="1" | '''Figure 1:''' Patch of three node triangular elements including the central triangle (M) and three adjacent triangles (1, 2 and 3)
47
|}
48
49
The curvature field within the central triangle can be expressed in terms of a constant assumed curvature field as
50
51
<span id="eq-2"></span>
52
{| class="formulaSCP" style="width: 100%; text-align: left;" 
53
|-
54
| 
55
{| style="text-align: left; margin:auto;" 
56
|-
57
| style="text-align: center;" | <math>{\boldsymbol \kappa }=\left\{ \begin{array}{c} \kappa _{xx}\\ \kappa _{yy}\\ \kappa _{xy} \end{array} \right\} =\hat{\boldsymbol \kappa } </math>
58
|}
59
| style="width: 5px;text-align: right;" | (2)
60
|}
61
62
where <math display="inline">{\boldsymbol \kappa }</math> is the curvature vector and <math display="inline">\hat{\boldsymbol \kappa }</math> is the assumed constant curvature field defined as
63
64
<span id="eq-3"></span>
65
{| class="formulaSCP" style="width: 100%; text-align: left;" 
66
|-
67
| 
68
{| style="text-align: left; margin:auto;" 
69
|-
70
| style="text-align: center;" | <math>\hat{\boldsymbol \kappa }={\frac{1}{A_{M}}}\int \int _{A_{M}}\left[ -{\frac{\partial ^{2} w}{\partial x^{2}}}-,{\frac{\partial ^{2}w}{\partial y^{2}}},-2{\frac{\partial ^{2}w}{\partial x\partial y}}\right] ^{T}dA </math>
71
|}
72
| style="width: 5px;text-align: right;" | (3)
73
|}
74
75
where <math display="inline">A_{M}</math> is the area of the central triangle in Figure [[#img-1|1]].
76
77
Integrating by parts Eq.([[#eq-3|3]]) and substituting the resulting expression for <math display="inline">\hat{\boldsymbol \kappa }</math> into Eq.([[#eq-2|2]]) gives the constant curvature field within the element as
78
79
<span id="eq-4"></span>
80
{| class="formulaSCP" style="width: 100%; text-align: left;" 
81
|-
82
| 
83
{| style="text-align: left; margin:auto;" 
84
|-
85
| style="text-align: center;" | <math>{\boldsymbol \kappa }={\frac{1}{A_{M}}}   {\displaystyle \oint _{\Gamma _{M}}}  \left[ \begin{array}{cc} -n_{x} & 0\\ 0 & -n_{y}\\ -n_{y} & -n_{x} \end{array} \right] \left\{ \begin{array}{c} \dfrac{\partial w}{\partial x}\\ \dfrac{\partial w}{\partial y} \end{array} \right\} d\Gamma </math>
86
|}
87
| style="width: 5px;text-align: right;" | (4)
88
|}
89
90
where <math display="inline">\Gamma _{M}</math> is the boundary of the central triangle.
91
92
Eq.([[#eq-4|4]]) defines the assumed constant curvature field within the central triangle in terms of the deflection gradient along the sides of the triangle. Equation (4) can be found to be equivalent to the standard conservation laws used in finite volume procedures as described in [28,29].
93
94
The computation of the line integral in Eq.([[#eq-4|4]]) poses a difficulty as the deflection gradient is discontinuous along the element sides. A simple method to overcome this problem is to compute the deflection gradient at the element sides as the average values of the gradient contributed by the two triangles sharing the side <span id='citeF-29'></span>[[#cite-29|[21,29]]]. Following this idea the constant curvature field with the element is computed as
95
96
<span id="eq-5"></span>
97
{| class="formulaSCP" style="width: 100%; text-align: left;" 
98
|-
99
| 
100
{| style="text-align: left; margin:auto;" 
101
|-
102
| style="text-align: center;" | <math>{\boldsymbol \kappa }={\frac{1}{A_{M}}}\sum _{j=1}^{3}{\frac{l_{j}^{M}}{2}}\left[ \begin{array}{cc} -n_{x}^{j} & 0\\ 0 & -n_{y}^{j}\\ -n_{y}^{j} & -n_{x}^{j} \end{array} \right] ^{M}\left[ {\nabla }L_{i}^{M}w_{i}^{M}+{\nabla }L_{i}^{j}w_{i} ^{j}\right] =\mathbf{B}_{b}\mathbf{w}^{p} </math>
103
|}
104
| style="width: 5px;text-align: right;" | (5)
105
|}
106
107
where <math display="inline">\mathbf{w}^{p}=[w_{1},w_{2},w_{3},w_{4},w_{5},w_{6}]^{T}</math> is the deflection vector of the six  nodes in the patch. In Eq.([[#eq-5|5]]) the sum extends over the three sides of the central element <math display="inline">M</math>, <math display="inline">l_{j}^{M}</math> are the lengths of the element sides and superindexes <math display="inline">M</math> and <math display="inline">j</math> refer to the central triangle and to each of the adjacent elements, respectively. The standard sum convention for repeated indexes is used.
108
109
Note that the constant curvature field is expressed in terms of the six nodes of the four element patch linked to the element <math display="inline">M</math>. The expression of the <math display="inline">3\times{6}</math> <math display="inline">\mathbf{B}_{b}</math> matrix can be found in [14,21].
110
111
The virtual work expression is written as
112
113
<span id="eq-6"></span>
114
{| class="formulaSCP" style="width: 100%; text-align: left;" 
115
|-
116
| 
117
{| style="text-align: left; margin:auto;" 
118
|-
119
| style="text-align: center;" | <math>\int \int _{A}\delta{\boldsymbol \kappa }^{T}\mathbf{m}\,dA=\int \int _{A}\delta w\,q\,dA </math>
120
|}
121
| style="width: 5px;text-align: right;" | (6)
122
|}
123
124
where '''m''' is the bending moment field related to the curvature by the standard constitutive equations
125
126
{| class="formulaSCP" style="width: 100%; text-align: left;" 
127
|-
128
| 
129
{| style="text-align: left; margin:auto;" 
130
|-
131
| style="text-align: center;" | <math>\mathbf{m}=[M_{xx},M_{yy},M_{xy}]^T = \mathbf{D}_{b}{\boldsymbol \kappa }\quad ,\quad \mathbf{D}_{b}={\frac{Eh^{3} }{(1-\nu ^{2})}}\left[ \begin{array}{ccc} 1 & \nu & 0\\ \nu & 1 & 0\\ 0 & 0 & \frac{1-\nu }{2} \end{array} \right] </math>
132
|}
133
| style="width: 5px;text-align: right;" | (7)
134
|}
135
136
In Eqs.(6) and (7) <math display="inline">h</math> is the plate thickness, <math display="inline">E</math> is the Young's modulus, <math display="inline">\nu </math> is the Poisson's ratio, <math display="inline">\delta{\boldsymbol \kappa }</math> and <math display="inline">\delta w</math> are the virtual curvatures and the virtual deflection, respectively, and <math display="inline">q</math> is a distributed vertical load.
137
138
Substituting the approximation for the vertical deflection and the assumed constant curvature field into ([[#eq-6|6]]) leads to the standard linear system of equations
139
140
{| class="formulaSCP" style="width: 100%; text-align: left;" 
141
|-
142
| 
143
{| style="text-align: left; margin:auto;" 
144
|-
145
| style="text-align: center;" | <math>\mathbf{K}\mathbf{w}=\mathbf{f}</math>
146
|}
147
| style="width: 5px;text-align: right;" | (8)
148
|}
149
150
where the stiffness matrix <math display="inline">\mathbf{K}</math> and the equivalent nodal force <math display="inline">\mathbf{f}</math> can be found by assembly of the element contributions given by
151
152
{| class="formulaSCP" style="width: 100%; text-align: left;" 
153
|-
154
| 
155
{| style="text-align: left; margin:auto;" 
156
|-
157
| style="text-align: center;" | <math>\mathbf{K}^{e}=\int \int _{A^{e}}\mathbf{B}_{b}^{T}\mathbf{D}_{b}\mathbf{B}_{b}dA </math>
158
|}
159
| style="width: 5px;text-align: right;" | (9)
160
|}
161
162
{| class="formulaSCP" style="width: 100%; text-align: left;" 
163
|-
164
| 
165
{| style="text-align: left; margin:auto;" 
166
|-
167
| style="text-align: center;" | <math>\mathbf{f}^{e}=\int \int _{A^{e}}q\left\{ \begin{array}{c} L_{1}^e\\ L_{2}^e\\ L_{3}^e \end{array} \right\} dA </math>
168
|}
169
| style="width: 5px;text-align: right;" | (10)
170
|}
171
172
Note that <math display="inline">\mathbf{K}^{e}</math> is  a <math display="inline">6\times{6}</math> matrix, whereas <math display="inline">\mathbf{f}^{e}</math> has the same structure than for the standard linear triangle.
173
174
The explicit form of <math display="inline">\mathbf{K}^{e}</math> and <math display="inline">\mathbf{f}^{e}</math> can be found in [14].
175
176
The resulting Basic Plate Triangle (BPT) has one degree of freedom per node and a wider bandwidth than the standard three node triangles as each triangular element is linked to its three neighbours through Eq.([[#eq-5|5]]).
177
178
Examples of the good performance of the BPT element for analysis of thin plates can be found in [14,21]. The extension of the BPT element to the analysis of shells yields the Basic Shell Triangle (BST) <span id='citeF-29'></span>[[#cite-29|29]]. Different applications of the BST element to linear and non linear analysis of shells are reported in 14,18&#8211;21,24,26,27].
179
180
The ideas used to derive the BPT element will now be extended to derive two families of Basic Shell Triangles using a total lagrangian description.
181
182
==3 BASIC THIN SHELL EQUATIONS USING A TOTAL LAGRANGIAN FORMULATION==
183
184
===3.1 Shell kinematics===
185
186
A summary of the most relevant hypothesis related to the kinematic behaviour of a thin shell are presented. Further details may be found in the wide literature dedicated to this field [8,9].
187
188
Consider a shell with undeformed middle surface occupying the domain <math display="inline">\Omega ^{0}</math> in <math display="inline">R^{3}</math> with a boundary <math display="inline">\Gamma ^{0}</math>. At each point of the middle surface a thickness <math display="inline">h^{0}</math> is defined. The positions <math display="inline">\mathbf{x}^{0}</math> and <math display="inline">\mathbf{x}</math> of a point in the undeformed and the deformed configurations can be respectively written as a function of the coordinates of the middle surface <math display="inline">{\boldsymbol \varphi }</math> and the normal <math display="inline">\mathbf{t}_{3}</math> at the point as
189
190
{| class="formulaSCP" style="width: 100%; text-align: left;" 
191
|-
192
| 
193
{| style="text-align: left; margin:auto;" 
194
|-
195
| style="text-align: right;" | <math>\mathbf{x}^{0}\left( \xi _{1},\xi _{2},\zeta \right)  </math>
196
| <math>  ={\boldsymbol \varphi }^{0}\left( \xi _{1},\xi _{2}\right) +\lambda \mathbf{t}_{3}^{0}</math>
197
| style="width: 5px;text-align: right;" | (11)
198
|-
199
| style="text-align: right;" | <math> \mathbf{x}\left( \xi _{1},\xi _{2},\zeta \right)  </math>
200
| <math>  ={\boldsymbol \varphi }\left( \xi  _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3} </math>
201
| style="width: 5px;text-align: right;" | (12)
202
|}
203
|}
204
205
where <math display="inline">\xi _{1},\xi _{2}</math> are curvilinear local coordinates defined over the middle surface of the shell, and <math display="inline">\zeta </math> is the distance in the undeformed configuration of the point to the middle surface. The product <math display="inline">\zeta \lambda </math> is the distance of the point to the middle surface measured on the deformed configuration. This implies a constant strain in the normal direction associated to the parameter <math display="inline">\lambda </math> relating the thickness at the present and initial configurations, i.e.
206
207
{| class="formulaSCP" style="width: 100%; text-align: left;" 
208
|-
209
| 
210
{| style="text-align: left; margin:auto;" 
211
|-
212
| style="text-align: center;" | <math>\lambda =\frac{h}{h^{0}} </math>
213
|}
214
| style="width: 5px;text-align: right;" | (13)
215
|}
216
217
A convective coordinate system is defined at each point as
218
219
{| class="formulaSCP" style="width: 100%; text-align: left;" 
220
|-
221
| 
222
{| style="text-align: left; margin:auto;" 
223
|-
224
| style="text-align: center;" | <math>\mathbf{g}_{i}\left( \mathbf{\xi }\right) =\frac{\partial \mathbf{x}}{} {\partial \xi _{i}}\qquad i=1,2,3 </math>
225
|}
226
| style="width: 5px;text-align: right;" | (14)
227
|}
228
229
{| class="formulaSCP" style="width: 100%; text-align: left;" 
230
|-
231
| 
232
{| style="text-align: left; margin:auto;" 
233
|-
234
| style="text-align: left;" | <math>\mathbf{g}_{\alpha }\left( \mathbf{\xi }\right)</math> <math> =\frac{\partial \left( \mathbf{\boldsymbol \varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \xi _{\alpha }}={\boldsymbol \varphi }_{^{\prime }\alpha }+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }\alpha }\quad \alpha=1,2</math>
235
|}
236
| style="width: 5px;text-align: right;" | (15)
237
|}
238
239
{| class="formulaSCP" style="width: 100%; text-align: left;" 
240
|-
241
| 
242
{| style="text-align: left; margin:auto;" 
243
|-
244
| style="text-align: left;" | <math> \mathbf{g}_{3}\left( \mathbf{\xi }\right)   =\frac{\partial \left( \mathbf{\boldsymbol \varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \zeta }=\lambda \mathbf{t}_{3} </math>
245
|}
246
| style="width: 5px;text-align: right;" | (16)
247
|}
248
249
This can be particularized for the points on the middle surface as
250
251
{| class="formulaSCP" style="width: 100%; text-align: left;" 
252
|-
253
| 
254
{| style="text-align: left; margin:auto;" 
255
|-
256
| style="text-align: left;" | <math>\mathbf{a}_{\alpha } =\mathbf{g}_{\alpha }\left( \zeta=0\right) ={\boldsymbol \varphi  }_{^{\prime }\alpha }</math>
257
|}
258
| style="width: 5px;text-align: right;" | (17)
259
|}
260
261
{| class="formulaSCP" style="width: 100%; text-align: left;" 
262
|-
263
| 
264
{| style="text-align: left; margin:auto;" 
265
|-
266
| style="text-align: left;" | <math> \mathbf{a}_{3}   =\mathbf{g}_{3}\left( \zeta=0\right) =\lambda  \mathbf{t}_{3} </math>
267
|}
268
| style="width: 5px;text-align: right;" | (18)
269
|}
270
271
272
The covariant (first fundamental form) and contravariant metric tensors of the middle surface are
273
274
{| class="formulaSCP" style="width: 100%; text-align: left;" 
275
|-
276
| 
277
{| style="text-align: left; margin:auto;" 
278
|-
279
| style="text-align: center;" | <math>a_{\alpha \beta }=\mathbf{a}_{\alpha }\cdot \mathbf{a}_{\beta } </math>
280
|}
281
| style="width: 5px;text-align: right;" | (19)
282
|}
283
284
{| class="formulaSCP" style="width: 100%; text-align: left;" 
285
|-
286
| 
287
{| style="text-align: left; margin:auto;" 
288
|-
289
| style="text-align: center;" | <math>a^{\alpha \beta }=\mathbf{a}^{\alpha }\cdot \mathbf{a}^{\beta }={\tilde{\boldsymbol \varphi }}_{^{\prime }\alpha }\cdot{\tilde{\boldsymbol \varphi }}_{^{\prime }\beta } </math>
290
|}
291
| style="width: 5px;text-align: right;" | (20)
292
|}
293
294
The Green-Lagrange strain vector of the middle surface points (membrane strains) is defined as
295
296
{| class="formulaSCP" style="width: 100%; text-align: left;" 
297
|-
298
| 
299
{| style="text-align: left; margin:auto;" 
300
|-
301
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_{m}=[\varepsilon _{m_{11}},\varepsilon _{m_{12}},\varepsilon _{m_{12}}]^{T} </math>
302
|}
303
| style="width: 5px;text-align: right;" | (21)
304
|}
305
306
with
307
308
<span id="eq-22"></span>
309
{| class="formulaSCP" style="width: 100%; text-align: left;" 
310
|-
311
| 
312
{| style="text-align: left; margin:auto;" 
313
|-
314
| style="text-align: center;" | <math>\varepsilon _{m_{ij}}=\frac{1}{2}(a_{ij}-\delta _{ij}) </math>
315
|}
316
| style="width: 5px;text-align: right;" | (22)
317
|}
318
319
The curvatures (second fundamental form) of the middle surface are obtained by
320
321
{| class="formulaSCP" style="width: 100%; text-align: left;" 
322
|-
323
| 
324
{| style="text-align: left; margin:auto;" 
325
|-
326
| style="text-align: center;" | <math>\kappa _{\alpha \beta }=\frac{1}{2}\left( {\boldsymbol \varphi }_{^{\prime }\alpha } \cdot \mathbf{t}_{3^{\prime }\beta }+{\boldsymbol \varphi }_{^{\prime }\beta }\cdot  \mathbf{t}_{3^{\prime }\alpha }\right) =- \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{{\prime }\alpha \beta }\quad , \quad \alpha ,\beta=1,2 </math>
327
|}
328
| style="width: 5px;text-align: right;" | (23)
329
|}
330
331
The deformation gradient tensor is
332
333
{| class="formulaSCP" style="width: 100%; text-align: left;" 
334
|-
335
| 
336
{| style="text-align: left; margin:auto;" 
337
|-
338
| style="text-align: center;" | <math>\mathbf{F=} [{x}_{{\prime }1},{x}_{{\prime }2},{x}_{{\prime }3}]=\left[ \begin{array}{ccc} {\boldsymbol \varphi }_{^{\prime }1}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime  }1} & {\boldsymbol \varphi }_{^{\prime }2}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }2} & \lambda \mathbf{t}_{3} \end{array} \right] </math>
339
|}
340
| style="width: 5px;text-align: right;" | (24)
341
|}
342
343
The product <math display="inline">\mathbf{F}^{T}\mathbf{F=U}^{2}=\mathbf{C}</math> (where <math display="inline">\mathbf{U}</math> is the right stretch tensor, and <math display="inline">\mathbf{C}</math> the right Cauchy-Green deformation tensor) can be written as
344
345
<span id="eq-25"></span>
346
{| class="formulaSCP" style="width: 100%; text-align: left;" 
347
|-
348
| 
349
{| style="text-align: left; margin:auto;" 
350
|-
351
| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{ccc} a_{11}+2\kappa _{11}\zeta \lambda & a_{12}+2\kappa _{12}\zeta \lambda & 0\\ a_{12}+2\kappa _{12}\zeta \lambda & a_{22}+2\kappa _{22}\zeta \lambda & 0\\ 0 & 0 & \lambda ^{2} \end{array} \right] </math>
352
|}
353
| style="width: 5px;text-align: right;" | (25)
354
|}
355
356
In the derivation of expression ([[#eq-25|25]]) the derivatives of the thickness ratio <math display="inline">\lambda _{^{\prime }a}</math> and the terms associated to <math display="inline">\zeta ^{2}</math> have been neglected.
357
358
Eq.([[#eq-25|25]]) shows that <math display="inline">\mathbf{U}^{2}</math> is not a unit tensor at the original configuration for curved surfaces (<math display="inline">\kappa _{ij}^{0}\neq{0}</math>). The changes of curvature of the middle surface are computed by
359
360
{| class="formulaSCP" style="width: 100%; text-align: left;" 
361
|-
362
| 
363
{| style="text-align: left; margin:auto;" 
364
|-
365
| style="text-align: center;" | <math>\chi _{ij}=\kappa _{ij}-\kappa _{ij}^{0} </math>
366
|}
367
| style="width: 5px;text-align: right;" | (26)
368
|}
369
370
Note that <math display="inline">\delta \chi _{ij}=\delta \kappa _{ij}</math>.
371
372
For computational convenience the following approximate expression (which is exact for initially flat surfaces) will be adopted
373
374
{| class="formulaSCP" style="width: 100%; text-align: left;" 
375
|-
376
| 
377
{| style="text-align: left; margin:auto;" 
378
|-
379
| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{ccc} a_{11}+2\chi _{11}\zeta \lambda & a_{12}+2\chi _{12}\zeta \lambda & 0\\ a_{12}+2\chi _{12}\zeta \lambda & a_{22}+2\chi _{22}\zeta \lambda & 0\\ 0 & 0 & \lambda ^{2} \end{array} \right] </math>
380
|}
381
| style="width: 5px;text-align: right;" | (27)
382
|}
383
384
This expression is useful to compute different lagrangian strain measures. An advantage of these measures is that they are associated to material fibres, what makes it easy to take into account material anisotropy. It is also useful to compute the eigen decomposition of <math display="inline">\mathbf{U}</math> as
385
386
{| class="formulaSCP" style="width: 100%; text-align: left;" 
387
|-
388
| 
389
{| style="text-align: left; margin:auto;" 
390
|-
391
| style="text-align: center;" | <math>\mathbf{U=}\sum _{\alpha=1}^{3}\lambda _{\alpha } \mathbf{r}_{\alpha } \otimes \mathbf{r}_{\alpha } </math>
392
|}
393
| style="width: 5px;text-align: right;" | (28)
394
|}
395
396
where <math display="inline">\lambda _{\alpha }</math> and <math display="inline">\mathbf{r}_{\alpha }</math> are the eigenvalues and eigenvectors of <math display="inline">\mathbf{U}</math>.
397
398
The resultant stresses  (axial forces and moments) are obtained by integrating across the original thickness the second Piola-Kirchhoff stress vector <math display="inline">{ \boldsymbol \sigma }</math> using the actual distance to the middle surface for  evaluating the bending moments, i.e.
399
400
<span id="eq-29"></span>
401
{| class="formulaSCP" style="width: 100%; text-align: left;" 
402
|-
403
| 
404
{| style="text-align: left; margin:auto;" 
405
|-
406
| style="text-align: center;" | <math>{\boldsymbol \sigma }_{m}\equiv \lbrack N_{11},N_{22},N_{12}]^{T}=\int _{h^{0}}{\boldsymbol \sigma } d\zeta </math>
407
|}
408
| style="width: 5px;text-align: right;" | (29)
409
|}
410
411
<span id="eq-30"></span>
412
{| class="formulaSCP" style="width: 100%; text-align: left;" 
413
|-
414
| 
415
{| style="text-align: left; margin:auto;" 
416
|-
417
| style="text-align: center;" | <math>{\boldsymbol \sigma }_{b}\equiv \lbrack M_{11},M_{22},M_{12}]^{T}=\int _{h^{0}}{\boldsymbol \sigma  }\lambda \zeta  d\zeta </math>
418
|}
419
| style="width: 5px;text-align: right;" | (30)
420
|}
421
422
With these values the virtual work can be written as
423
424
<span id="eq-31"></span>
425
{| class="formulaSCP" style="width: 100%; text-align: left;" 
426
|-
427
| 
428
{| style="text-align: left; margin:auto;" 
429
|-
430
| style="text-align: center;" | <math>\int \int _{A^{0}}\left[ \delta{\boldsymbol \varepsilon }_{m}^{T}{\boldsymbol \sigma }_{m}+\delta{\boldsymbol \kappa  }^{T}{\boldsymbol \sigma }_{b}\right] dA=\int \int _{A^{0}}\delta \mathbf{u}^{T}\mathbf{t}dA </math>
431
|}
432
| style="width: 5px;text-align: right;" | (31)
433
|}
434
435
where <math display="inline">\delta \mathbf{u}</math> are virtual displacements, <math display="inline">\delta{\boldsymbol \varepsilon }_{m}</math> is the virtual Green-Lagrange membrane strain vector, <math display="inline">\delta{\boldsymbol \kappa }</math> are the virtual curvatures and <math display="inline">\mathbf{t}</math> are the surface loads. Other load types can be easily included into ([[#eq-31|31]]).
436
437
===3.2 Constitutive models===
438
439
In order to treat plasticity at finite strains an adequate stress-strain pair must be used. The Hencky measures will be adopted here. The (logarithmic) strains are defined as
440
441
{| class="formulaSCP" style="width: 100%; text-align: left;" 
442
|-
443
| 
444
{| style="text-align: left; margin:auto;" 
445
|-
446
| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=}\left[ \begin{array}{ccc} \varepsilon _{11} & \varepsilon _{21} & 0\\ \varepsilon _{12} & \varepsilon _{22} & 0\\ 0 & 0 & \varepsilon _{33} \end{array} \right] =\sum _{\alpha=1}^{3}\ln \left( \lambda _{\alpha }\right) \mathbf{r}_{\alpha }\otimes \mathbf{r}_{\alpha } </math>
447
|}
448
| style="width: 5px;text-align: right;" | (32)
449
|}
450
451
Two types of material models are considered here: an elastic-plastic material associated to thin rolled metal sheets and a hyper-elastic material for rubbers.
452
453
In the case of metals, where the elastic strains are small, the use of a logarithmic strain measure reasonably allows to adopt an additive decomposition of elastic and plastic components as
454
455
{| class="formulaSCP" style="width: 100%; text-align: left;" 
456
|-
457
| 
458
{| style="text-align: left; margin:auto;" 
459
|-
460
| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=E}_{\ln }^{e}+\mathbf{E}_{\ln }^{p} </math>
461
|}
462
| style="width: 5px;text-align: right;" | (33)
463
|}
464
465
A constant linear relationship between the (plane) Hencky stresses and the logarithmic elastic strains is  adopted giving
466
467
{| class="formulaSCP" style="width: 100%; text-align: left;" 
468
|-
469
| 
470
{| style="text-align: left; margin:auto;" 
471
|-
472
| style="text-align: center;" | <math>\mathbf{T}=\mathbf{CE}_{\ln }^{e} </math>
473
|}
474
| style="width: 5px;text-align: right;" | (34)
475
|}
476
477
These constitutive equations are integrated using a standard return algorithm. The following Mises-Hill [30] yield function with non-linear isotropic hardening is chosen here
478
479
{| class="formulaSCP" style="width: 100%; text-align: left;" 
480
|-
481
| 
482
{| style="text-align: left; margin:auto;" 
483
|-
484
| style="text-align: center;" | <math>\left( G+H\right) \;T_{11}^{2}+\left( F+H\right) \;T_{22}^{2} -2H\;T_{11}T_{22}+2N\;T_{12}^{2}=\sigma _0\left(e_{0}+e^{p}\right) ^{n} </math>
485
|}
486
| style="width: 5px;text-align: right;" | (35)
487
|}
488
489
where <math display="inline">F, G, H</math> and <math display="inline">N</math> define the non-isotropic shape of the yield surface and the parameters <math display="inline">\sigma _{0}</math>, <math display="inline">e_{0}</math> and <math display="inline">n</math> define its size as a function of the effective plastic strain <math display="inline">e^{p}</math>.
490
491
The simple Mises-Hill yield function  allows, as a first approximation, to treat rolled thin metal sheets with planar and transversal anisotropy.
492
493
For the case of rubbers, the Ogden [31] model extended to the compressible range is considered. The material behaviour is characterized by the strain energy density per unit undeformed volume defined as
494
495
{| class="formulaSCP" style="width: 100%; text-align: left;" 
496
|-
497
| 
498
{| style="text-align: left; margin:auto;" 
499
|-
500
| style="text-align: center;" | <math>\psi =\frac{K}{2}\left( \ln J\right) ^{2}+\sum _{p=1}^{N}\frac{\mu _{p}}{} {\alpha _{p}}\left[ J^{-\frac{\alpha _{p}}{3}}\left( \sum _{i=1}^{3}\lambda  _{i}^{\alpha _{p}-1}\right) -3\right] </math>
501
|}
502
| style="width: 5px;text-align: right;" | (36)
503
|}
504
505
where <math display="inline">K</math> is the bulk modulus of the material, <math display="inline">J</math> is the determinant of <math display="inline">\mathbf{U}</math>, <math display="inline">N</math>, <math display="inline">\mu _{i}</math> and <math display="inline">\alpha _{i}</math> are material parameters, <math display="inline">\mu _{i}\,,\,\alpha _{i}</math> are real numbers such that <math display="inline">\mu _{i}\alpha _{i}>0</math> <math display="inline"> (\forall i=1,N)</math> and <math display="inline">N</math> is a positive integer.
506
507
The stress measures associated to the principal logarithmic strains are denoted by <math display="inline">\beta _{i}</math>. They can be computed noting that
508
509
{| class="formulaSCP" style="width: 100%; text-align: left;" 
510
|-
511
| 
512
{| style="text-align: left; margin:auto;" 
513
|-
514
| style="text-align: center;" | <math>\beta _{i}=\frac{\partial \psi \left( \Lambda \right) }{\partial \left( \ln \lambda _{i}\right) }=K\left( \ln J\right) +\lambda _{i}\sum _{p=1}^{N} \mu _{p}J^{-\frac{\alpha _{p}}{3}}\left( \lambda _{i}^{\alpha _{p}-1}-\frac{1}{} {3}\frac{1}{\lambda _{i}}\sum _{j=1}^{3}\lambda _{j}^{\alpha _{p}}\right) </math>
515
|}
516
| style="width: 5px;text-align: right;" | (37)
517
|}
518
519
we define now
520
521
{| class="formulaSCP" style="width: 100%; text-align: left;" 
522
|-
523
| 
524
{| style="text-align: left; margin:auto;" 
525
|-
526
| style="text-align: center;" | <math>a^{p}=\sum _{j=1}^{3}\lambda _{j}^{\alpha _{p}} </math>
527
|}
528
| style="width: 5px;text-align: right;" | (38)
529
|}
530
531
which gives
532
533
{| class="formulaSCP" style="width: 100%; text-align: left;" 
534
|-
535
| 
536
{| style="text-align: left; margin:auto;" 
537
|-
538
| style="text-align: center;" | <math>\beta _{i}=K\left( \ln J\right) +\sum _{p=1}^{N}\mu _{p}J^{-\frac{\alpha _{p} }{3}}\left( \lambda _{i}^{\alpha _{p}}-\frac{1}{3}a_{p}\right) </math>
539
|}
540
| style="width: 5px;text-align: right;" | (39)
541
|}
542
543
The values of <math display="inline">\beta _{i}</math>, expressed in the principal strains directions, allow to evaluate the Hencky stresses in the convective coordinate system as
544
545
{| class="formulaSCP" style="width: 100%; text-align: left;" 
546
|-
547
| 
548
{| style="text-align: left; margin:auto;" 
549
|-
550
| style="text-align: center;" | <math>\mathbf{T}=\sum _{i=1}^{3}\beta _{i}\;\mathbf{r}_{i}\otimes \mathbf{r}_{i} </math>
551
|}
552
| style="width: 5px;text-align: right;" | (40)
553
|}
554
555
The Hencky stress tensor <math display="inline">\mathbf{T}</math> can be easily particularized for the plane stress case.
556
557
We define the rotated Hencky and second Piola-Kirchhoff stress tensors as
558
559
{| class="formulaSCP" style="width: 100%; text-align: left;" 
560
|-
561
| 
562
{| style="text-align: left; margin:auto;" 
563
|-
564
| style="text-align: left;" | <math>\mathbf{T}_{L}  </math>
565
| <math>  =\mathbf{R}_{L}^{T}\;\mathbf{T\;R}_{L}</math>
566
|}
567
| style="width: 5px;text-align: right;" | (41)
568
|-
569
| 
570
{| style="text-align: left; margin:auto;" 
571
|-
572
| style="text-align: left;" | <math> \mathbf{S}_{L}  </math>
573
| <math>  =\mathbf{R}_{L}^{T}\;\mathbf{S\;R}_{L} </math>
574
|}
575
| style="width: 5px;text-align: right;" | (42)
576
|}
577
578
579
where <math display="inline">\mathbf{R}_{L}</math> is the rotation tensor obtained from the eigenvectors of <math display="inline">\mathbf{U}</math> given by
580
581
{| class="formulaSCP" style="width: 100%; text-align: left;" 
582
|-
583
| 
584
{| style="text-align: left; margin:auto;" 
585
|-
586
| style="text-align: center;" | <math>\mathbf{R}_{L}=\left[ \begin{array}{ccc} \mathbf{r}_{1} & \mathbf{r}_{2} & \mathbf{r}_{3} \end{array} \right] </math>
587
|}
588
| style="width: 5px;text-align: right;" | (43)
589
|}
590
591
The relationship between the rotated Hencky and Piola-Kirchhoff stresses is
592
593
{| class="formulaSCP" style="width: 100%; text-align: left;" 
594
|-
595
| 
596
{| style="text-align: left; margin:auto;" 
597
|-
598
| style="text-align: right;" | <math>\left[ S_{L}\right] _{\alpha \alpha }  </math>
599
| <math>  =\frac{1}{\lambda _{\alpha }^{2} }\left[ T_{L}\right] _{\alpha \alpha }</math>
600
|-
601
| style="text-align: right;" | <math> \left[ S_{L}\right] _{\alpha \beta }  </math>
602
| <math>  =\frac{\ln \left( \lambda _{\alpha  }/\lambda _{\beta }\right) }{\frac{1}{2}\left( \lambda _{\alpha }^{2} -\lambda _{\beta }^{2}\right) }\left[ T_{L}\right] _{\alpha \beta } </math>
603
|}
604
| style="width: 5px;text-align: right;" | (44)
605
|}
606
607
The second Piola-Kirchhoff stress tensor can be computed by
608
609
{| class="formulaSCP" style="width: 100%; text-align: left;" 
610
|-
611
| 
612
{| style="text-align: left; margin:auto;" 
613
|-
614
| style="text-align: center;" | <math>\mathbf{S}=\mathbf{R}_{L}\;\mathbf{S}_{L}\mathbf{\;R}_{L}^{T} </math>
615
|}
616
| style="width: 5px;text-align: right;" | (45)
617
|}
618
619
The second Piola-Kirchhoff stress vector <math display="inline">{\boldsymbol \sigma }</math> of Eqs.([[#eq-29|29]]&#8211;[[#eq-30|30]]) can be readily extracted from the <math display="inline">\mathbf{S}</math> tensor.
620
621
==4 TOTAL LAGRANGIAN FORMULATION OF THE BASIC SHELL TRIANGLE==
622
623
===4.1 Definition of the element geometry and discretization of the displacement field===
624
625
The rotation-free BST element has three nodes with three displacement degrees of freedom at each node. As before an element patch is defined by the central triangle  and the three adjacent elements (Figure [[#img-1|1]]). This patch helps to define the curvature field within the central triangle (the BST element) in terms of the displacement of the six patch nodes.
626
627
The node-ordering in the patch is the following (see Figure [[#img-1|1]])
628
629
* The nodes in the main element (M) are numbered locally as 1, 2 and 3. They are defined counter-clockwise around the positive normal.
630
631
* The sides in the main element are numbered locally as 1, 2, and 3. They are defined by the local node opposite to the side.
632
633
* The adjacent elements (which are part of the cell) are numbered with the number associated to the common side.
634
635
* The extra nodes of the cell are numbered locally as 4, 5 and 6, corresponding to nodes on adjacent elements opposite to sides 1, 2  and 3 respectively.
636
637
* The connectivities in the adjacent elements are defined beginning with the extra node as shown in Table 1.
638
639
640
<div class="center" style="font-size: 75%;">
641
'''Table 1'''. Element numbering and nodal connectivities of the four elements patch of Figure 1.</div>
642
643
{| class="wikitable" style="text-align: left; margin: 1em auto;"
644
|- style="border-top: 2px solid;"
645
| '''Element''' 
646
| N1 
647
| N2 
648
| N3
649
|- style="border-top: 2px solid;"
650
|  '''M''' 
651
| 1 
652
| 2 
653
| 3
654
|- style="border-top: 2px solid;"
655
|  '''1''' 
656
| 4 
657
| 3 
658
| 2
659
|- style="border-top: 2px solid;"
660
|  '''2''' 
661
| 5 
662
| 1 
663
| 3
664
|- style="border-top: 2px solid;border-bottom: 2px solid;"
665
|  '''3''' 
666
| 6 
667
| 2 
668
| 1
669
670
|}
671
672
The following local cartesian coordinate system can be defined for the patch. In the main element the unit vector <math display="inline">\mathbf{t}_{1}</math>(associated to the local coordinate <math display="inline">\xi _{1}</math>) is directed along side 3 (from node 1 to node 2), <math display="inline">\mathbf{t}_{3}</math> (associated to the coordinate <math display="inline">\zeta </math>) is the unit normal to the plane, and finally <math display="inline">\mathbf{t}_{2}=\mathbf{t}_{3}\times \mathbf{t}_{1}</math> (associated to the coordinate <math display="inline">\xi _{2}</math>).
673
674
The coordinates and the displacements are linearly interpolated within each three node triangle in the mesh in the standard manner, i.e.
675
676
{| class="formulaSCP" style="width: 100%; text-align: left;" 
677
|-
678
| 
679
{| style="text-align: left; margin:auto;" 
680
|-
681
| style="text-align: center;" | <math>{\boldsymbol \varphi } = \sum \limits _{i=1}^{3} L_{i}^e {\boldsymbol \varphi }_{i} = \sum \limits _{i=1}^{3} L_{i}^e ({\boldsymbol \varphi }^{0}_{i} + \mathbf{u}_{i}) </math>
682
|}
683
| style="width: 5px;text-align: right;" | (46)
684
|}
685
686
{| class="formulaSCP" style="width: 100%; text-align: left;" 
687
|-
688
| 
689
{| style="text-align: left; margin:auto;" 
690
|-
691
| style="text-align: center;" | <math>\mathbf{u}=\left\{ \begin{array}{c} u_{1}\\ u_{2}\\ u_{3} \end{array} \right\} =\sum \limits _{i=1}^{3}L_{i}^e\mathbf{u}_{i}\quad ,\quad \mathbf{u}_{i}=\left\{ \begin{array}{c} u_{1}\\ u_{2}\\ u_{3} \end{array} \right\} _{i} </math>
692
|}
693
| style="width: 5px;text-align: right;" | (47)
694
|}
695
696
In above <math display="inline">{\boldsymbol \varphi }_{i}</math> and <math display="inline">\mathbf{u}_{i}</math> contain respectively the three coordinates and the three displacements of node <math display="inline">i</math>.
697
698
===4.2 Computation of the membrane strains===
699
700
The Green-Lagrange membrane strains are expressed by substituting the linear displacement interpolation into Eq.([[#eq-22|22]]). This gives
701
702
{| class="formulaSCP" style="width: 100%; text-align: left;" 
703
|-
704
| 
705
{| style="text-align: left; margin:auto;" 
706
|-
707
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_{m}=\frac{1}{2}\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }1}\cdot \boldsymbol \varphi _{^{\prime  }1}-1 \\ \boldsymbol \varphi _{^{\prime }2}\cdot \boldsymbol \varphi _{^{\prime  }2}-1 \\ 2\boldsymbol \varphi _{^{\prime }1}\cdot \boldsymbol \varphi _{^{\prime }2} \end{array}\right] </math>
708
|}
709
| style="width: 5px;text-align: right;" | (48)
710
|}
711
712
The membrane strain field is constant within each triangle similarly as in the standard CST element. The variation of the membrane strains is simply obtained by
713
714
{| class="formulaSCP" style="width: 100%; text-align: left;" 
715
|-
716
| 
717
{| style="text-align: left; margin:auto;" 
718
|-
719
| style="text-align: center;" | <math>\delta{\boldsymbol \varepsilon }_{m}=\mathbf{B}_{m}\delta \mathbf{a}^{e} </math>
720
|}
721
| style="width: 5px;text-align: right;" | (49)
722
|}
723
724
with
725
726
<span id="eq-50"></span>
727
{| class="formulaSCP" style="width: 100%; text-align: left;" 
728
|-
729
| 
730
{| style="text-align: left; margin:auto;" 
731
|-
732
| style="text-align: center;" | <math>\mathbf{B}_{m}=[\mathbf{B}_{m_{1}},\mathbf{B}_{m_{2}},\mathbf{B}_{m_{3}} ]\quad ,\quad \mathbf{a}^{e}=\left\{ \begin{array}{c} \mathbf{u}_{1}\\ \mathbf{u}_{2}\\ \mathbf{u}_{3} \end{array} \right\} </math>
733
|}
734
| style="width: 5px;text-align: right;" | (50)
735
|}
736
737
and
738
739
<span id="eq-51"></span>
740
{| class="formulaSCP" style="width: 100%; text-align: left;" 
741
|-
742
| 
743
{| style="text-align: left; margin:auto;" 
744
|-
745
| style="text-align: center;" | <math>\begin{array}{c} \\ \mathbf{B}_{m_{i}}\\ 3\times{3} \end{array} =\left[ \begin{array}{c} L_{i,1}^M\boldsymbol \varphi _{^{\prime }1}^{T}\\ L_{i,2}^M\boldsymbol \varphi _{^{\prime }2}^{T}\\ L_{i,1}^M\boldsymbol \varphi _{^{\prime }2}^{T}+L_{i,2}^M\boldsymbol \varphi _{^{\prime }1}^{T} \end{array} \right]  </math>
746
|}
747
| style="width: 5px;text-align: right;" | (51)
748
|}
749
750
===4.3 Computation of bending strains (curvatures)===
751
752
We will assume the following constant curvature field within each element
753
754
{| class="formulaSCP" style="width: 100%; text-align: left;" 
755
|-
756
| 
757
{| style="text-align: left; margin:auto;" 
758
|-
759
| style="text-align: center;" | <math>\kappa _{\alpha \beta }=\hat{\kappa }_{\alpha \beta } </math>
760
|}
761
| style="width: 5px;text-align: right;" | (52)
762
|}
763
764
where <math display="inline">\hat{\kappa }_{\alpha \beta }</math> is the assumed constant curvature field defined by
765
766
<span id="eq-53"></span>
767
{| class="formulaSCP" style="width: 100%; text-align: left;" 
768
|-
769
| 
770
{| style="text-align: left; margin:auto;" 
771
|-
772
| style="text-align: center;" | <math>\hat{\kappa }_{\alpha \beta }=-\frac{1}{A_{M}^{0}}\int _{A_{M}^{0}}\mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }\beta \alpha }\;dA^{0} </math>
773
|}
774
| style="width: 5px;text-align: right;" | (53)
775
|}
776
777
where <math display="inline">A_{M}^{0}</math> is the area (in the original configuration) of the central element in the patch.
778
779
Substituting Eq.(53) into (52) and integrating by parts the area integral gives the curvature vector within the element in terms of the following line integral
780
781
<span id="eq-54"></span>
782
{| class="formulaSCP" style="width: 100%; text-align: left;" 
783
|-
784
| 
785
{| style="text-align: left; margin:auto;" 
786
|-
787
| style="text-align: center;" | <math>{\boldsymbol \kappa }=\left\{ \begin{array}{c} \kappa _{11}\\ \kappa _{22}\\ 2\kappa _{12} \end{array} \right\} =\frac{1}{A_{M}^{0}}   {\displaystyle \oint _{\Gamma _{M}^{0}}}  \left[ \begin{array}{cc} -n_{1} & 0\\ 0 & -n_{2}\\ -n_{2} & -n_{1} \end{array} \right] \left[ \begin{array}{c} \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }1}\\ \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }2} \end{array} \right] d\Gamma </math>
788
|}
789
| style="width: 5px;text-align: right;" | (54)
790
|}
791
792
where <math display="inline">n_{i}</math> are the components (in the local system) of the normals to the element sides in the initial configuration <math display="inline">\Gamma _{M}^{0}</math>.
793
794
For the definition of the normal vector <math display="inline">\mathbf{t}_{3}</math>, the linear interpolation over the central element is used. In this case the tangent plane components are
795
796
{| class="formulaSCP" style="width: 100%; text-align: left;" 
797
|-
798
| 
799
{| style="text-align: left; margin:auto;" 
800
|-
801
| style="text-align: center;" | <math>{\boldsymbol \varphi }_{^{\prime }\alpha } = \sum _{i=1}^{3} L_{i,\alpha }^M {\boldsymbol \varphi }_{i} \quad ,\quad \alpha=1,2 </math>
802
|}
803
| style="width: 5px;text-align: right;" | (55)
804
|}
805
806
<span id="eq-56"></span>
807
{| class="formulaSCP" style="width: 100%; text-align: left;" 
808
|-
809
| 
810
{| style="text-align: left; margin:auto;" 
811
|-
812
| style="text-align: center;" | <math>\mathbf{t}_{3}=\frac{{\boldsymbol \varphi }_{\prime{1}}\times{\boldsymbol \varphi }_{\prime{2}}}{\left\vert {\boldsymbol \varphi }_{\prime{1}}\times{\boldsymbol \varphi }_{\prime{2}}\right\vert }=\lambda \;_{1}\times{\boldsymbol \varphi }_{2} </math>
813
|}
814
| style="width: 5px;text-align: right;" | (56)
815
|}
816
817
From these expressions it is also possible to compute in the original configuration the element area <math display="inline">A^{0}_{M}</math>, the outer normals <math display="inline">\left( n_{1} ,n_{2}\right) ^{i}</math> at each side and the side lengths <math display="inline">l_{i}^{M}</math>. Eq.([[#eq-56|56]]) also allows to evaluate the thickness ratio <math display="inline">\lambda </math> in the deformed configuration and the actual normal <math display="inline">\mathbf{t}_{3}</math>.
818
819
In order to compute the line integral of equation ([[#eq-54|54]]) the averaging procedure described in Section 2 is used. Hence along each side of the triangle the average value of <math display="inline">{\boldsymbol \varphi }_{^{\prime }\alpha }</math> between the main triangle and the adjacent one is taken leading to
820
821
{| class="formulaSCP" style="width: 100%; text-align: left;" 
822
|-
823
| 
824
{| style="text-align: left; margin:auto;" 
825
|-
826
| style="text-align: center;" | <math>{\boldsymbol \kappa }=\frac{1}{A^{0}_{M}}\sum _{I=1}^{3}\left[ \begin{array}{cc} -n_{1}^{i} & 0\\ 0 & n_{2}^{i}\\ n_{2}^{i} & -n_{1}^{i} \end{array} \right] \left[ \begin{array}{c} \mathbf{t}_{3}\cdot \frac{1}{2}\left( \mathbf{\boldsymbol \varphi }_{^{\prime }1} ^{M}+\mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}\right)\\ \mathbf{t}_{3}\cdot \frac{1}{2}\left( \mathbf{\boldsymbol \varphi }_{^{\prime }2} ^{M}+\mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}\right) \end{array} \right] </math>
827
|}
828
| style="width: 5px;text-align: right;" | (57)
829
|}
830
831
where the sum extends over the three elements adjacent to the central triangle <math display="inline">M</math>.
832
833
Noting that <math display="inline">\mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }\alpha }^{M}=0</math> in the main triangle it can be found <span id='citeF-24'></span>[[#cite-24|[24]]]
834
835
<span id="eq-58"></span>
836
{| class="formulaSCP" style="width: 100%; text-align: left;" 
837
|-
838
| 
839
{| style="text-align: left; margin:auto;" 
840
|-
841
| style="text-align: center;" | <math>{\boldsymbol \kappa }=\sum _{i=1}^{3}\left[ \begin{array}{cc} L_{i,1}^M & 0 \\         0 & L_{i,2}^M \\ L_{i,2}^M & L_{i,1}^M  \end{array} \right] \left[ \begin{array}{c} \mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}\\ \mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i} \end{array} \right] </math>
842
|}
843
| style="width: 5px;text-align: right;" | (58)
844
|}
845
846
This can be seen as the projection of the local derivatives in the adjacent triangles <math display="inline">\mathbf{\boldsymbol \varphi }_{^{\prime }\alpha }^{i}</math> (where index <math display="inline">i</math> denotes values associated to the adjacent elements) over the normal to the main triangle <math display="inline">\mathbf{t}_{3}</math>. As the triangles have a common side, <math display="inline">\mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }s}^{i}=0</math>, where <math display="inline">\mathbf{\boldsymbol \varphi }_{^{\prime } s}^{i}</math> is the derivative along the side. Hence only the derivative along the side normal (<math display="inline">\mathbf{\boldsymbol \varphi }_{^{\prime }n}^{i}</math>) has non-zero component over <math display="inline">\mathbf{t}_{3}</math>. This gives
847
848
<span id="eq-59"></span>
849
{| class="formulaSCP" style="width: 100%; text-align: left;" 
850
|-
851
| 
852
{| style="text-align: left; margin:auto;" 
853
|-
854
| style="text-align: center;" | <math>\left[ \begin{array}{c} \mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}\\ \mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i} \end{array} \right] =\left( \mathbf{t}_{3}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime  }n}^{i}\right)\mathbf{n}^{i} </math>
855
|}
856
| style="width: 5px;text-align: right;" | (59)
857
|}
858
859
An alternative form to express the curvatures, which is useful when their variations are needed, is to define the vectors
860
861
<span id="eq-60"></span>
862
{| class="formulaSCP" style="width: 100%; text-align: left;" 
863
|-
864
| 
865
{| style="text-align: left; margin:auto;" 
866
|-
867
| style="text-align: center;" | <math>\mathbf{h}_{ij}=\sum _{k=1}^{3}\frac{1}{2}\left( L_{k,i}^{M}{\boldsymbol \varphi  }_{^{\prime }j}^{k}+L_{k,j}^{M}{\boldsymbol \varphi }_{\prime i}^{k}\right) </math>
868
|}
869
| style="width: 5px;text-align: right;" | (60)
870
|}
871
872
This gives
873
874
<span id="eq-61"></span>
875
{| class="formulaSCP" style="width: 100%; text-align: left;" 
876
|-
877
| 
878
{| style="text-align: left; margin:auto;" 
879
|-
880
| style="text-align: center;" | <math>\kappa _{ij}=\mathbf{h}_{ij}\cdot \mathbf{t}_{3}</math>
881
|}
882
| style="width: 5px;text-align: right;" | (61)
883
|}
884
885
The last expression allows to interpret the curvatures as the projections of the vectors <math display="inline">\mathbf{h}_{ij}</math> over the normal of the central element. The variation of the curvatures can be obtained as
886
887
<span id="eq-62"></span>
888
{| class="formulaSCP" style="width: 100%; text-align: left;" 
889
|-
890
| 
891
{| style="text-align: left; margin:auto;" 
892
|-
893
| style="text-align: center;" | <math>\delta{\boldsymbol \kappa }=\sum _{i=1}^{3}\left\{ \left[ \begin{array}{cc} L_{i,1}^{M} & 0\\ 0 & L_{i,2}^{M}\\ L_{i,2}^{M} & L_{i,1}^{M} \end{array} \right] \sum _{J=1}^{3}\left[ \begin{array}{c} L_{j,1}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}_{j}^{i})\\ N_{j,2}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}_{j}^{i}) \end{array} \right] -2\left[ \begin{array}{c} (L_{i,1}^{M}\rho _{11}^{1}+L_{i,2}^{M}\rho _{11}^{2})\\ (L_{i,1}^{M}\rho _{22}^{1}+L_{i,2}^{M}\rho _{22}^{2})\\ (L_{i,1}^{M}\rho _{12}^{1}+L_{i,2}^{M}\rho _{12}^{2}) \end{array} \right] (\mathbf{t}_{3}\cdot \delta \mathbf{u}_{i}^{M})\right\} </math>
894
|}
895
| style="width: 5px;text-align: right;" | (62)
896
|}
897
898
where the projections of the vectors <math display="inline">\mathbf{h}_{ij}</math> over the contravariant base vectors <math display="inline">\tilde{\boldsymbol \varphi }_{^{\prime }\alpha }</math> have been included
899
900
<span id="eq-63"></span>
901
{| class="formulaSCP" style="width: 100%; text-align: left;" 
902
|-
903
| 
904
{| style="text-align: left; margin:auto;" 
905
|-
906
| style="text-align: center;" | <math>\rho _{ij}^{\alpha }=\mathbf{h}_{ij}\cdot \tilde{\boldsymbol \varphi }_{^{\prime }\alpha } \quad ,\quad \alpha ,i,j=1,2</math>
907
|}
908
| style="width: 5px;text-align: right;" | (63)
909
|}
910
911
with
912
913
{| class="formulaSCP" style="width: 100%; text-align: left;" 
914
|-
915
| 
916
{| style="text-align: left; margin:auto;" 
917
|-
918
| style="text-align: right;" | <math>\mathbf{\tilde{\boldsymbol \varphi }}_{^{\prime }1}  </math>
919
| <math>  =\lambda \;\mathbf{\boldsymbol \varphi }_{^{\prime }2}\times \mathbf{t}_{3}</math>
920
|}
921
| style="width: 5px;text-align: right;" | (64)
922
|}
923
924
{| class="formulaSCP" style="width: 100%; text-align: left;" 
925
|-
926
| 
927
{| style="text-align: left; margin:auto;" 
928
|-
929
| style="text-align: right;" | <math> \mathbf{\tilde{\boldsymbol \varphi }}_{^{\prime }2}  </math>
930
| <math>  =-\lambda \;\mathbf{\boldsymbol \varphi  }_{^{\prime }1}\times \mathbf{t}_{3} </math>
931
|}
932
| style="width: 5px;text-align: right;" | (65)
933
|}
934
935
936
In above expressions superindexes in <math display="inline">L_{j}^k</math> and <math display="inline">\delta \mathbf{u}_{j}^k</math> refer to element numbers whereas subscripts denote node numbers. As usual the superindex <math display="inline">M</math> denotes values in the central triangle (Figure [[#img-1|1]]). Note that as expected the curvatures (and their variations) in the central element are a function of the nodal displacements of the six nodes in the four elements patch. Note also that
937
938
{| class="formulaSCP" style="width: 100%; text-align: left;" 
939
|-
940
| 
941
{| style="text-align: left; margin:auto;" 
942
|-
943
| style="text-align: center;" | <math>\lambda ={\frac{h}{h^{0}}}={\frac{A_{M}^{0}}{A_{M}}} </math>
944
|}
945
| style="width: 5px;text-align: right;" | (66)
946
|}
947
948
Details of the derivation of Eq.([[#eq-62|62]]) can be found in [27].
949
950
Eq.([[#eq-62|62]]) can be rewritten in the form
951
952
{| class="formulaSCP" style="width: 100%; text-align: left;" 
953
|-
954
| 
955
{| style="text-align: left; margin:auto;" 
956
|-
957
| style="text-align: center;" | <math>\delta{\boldsymbol \kappa }=\mathbf{B}_{b}\delta \mathbf{a}^{p} </math>
958
|}
959
| style="width: 5px;text-align: right;" | (67)
960
|}
961
962
where
963
964
<span id="eq-68"></span>
965
{| class="formulaSCP" style="width: 100%; text-align: left;" 
966
|-
967
| 
968
{| style="text-align: left; margin:auto;" 
969
|-
970
| style="text-align: center;" | <math>\begin{array}{c} \\ \delta \mathbf{a}^{p}\\ 18\times{1} \end{array} =[\delta \mathbf{u}_{1}^{T},\delta \mathbf{u}_{2}^{T},\delta \mathbf{u}_{3} ^{T},\delta \mathbf{u}_{4}^{T},\delta \mathbf{u}_{5}^{T},\delta \mathbf{u}_{6}^{T}]^{T}</math>
971
|}
972
| style="width: 5px;text-align: right;" | (68)
973
|}
974
975
is the virtual displacement vector of the patch
976
977
<span id="eq-69"></span>
978
{| class="formulaSCP" style="width: 100%; text-align: left;" 
979
|-
980
| 
981
{| style="text-align: left; margin:auto;" 
982
|-
983
| style="text-align: center;" | <math>\mathbf{B}_{b}=[\mathbf{B}_{b1},\mathbf{B}_{b2}\cdots ,\mathbf{B}_{b6}]</math>
984
|}
985
| style="width: 5px;text-align: right;" | (69)
986
|}
987
988
is the curvature matrix relating the virtual curvatures within the central element and the 18 virtual displacements of the six nodes in the patch.
989
990
The form of matrix <math display="inline">\mathbf{B}_{b}</math> is given in the Appendix.
991
992
==5 ENHANCED BASIC SHELL TRIANGLE==
993
994
An enhanced version of the BST element (termed EBST) has been recently proposed by Flores and Oñate [27]. The main features of the element formulation are the following:
995
996
<ol>
997
998
<li>The geometry of the patch formed by the central element and the three adjacent elements is ''quadratically interpolated'' from the position of the six nodes in the patch. </li>
999
1000
<li>The membrane strains are assumed to vary ''linearly'' within the central triangle and are expressed in terms of the (continuous) values of the deformation gradient at the mid side points of the triangle. </li>
1001
1002
<li>The assumed ''constant curvature'' field within the central triangle is obtained by expression ([[#eq-54|54]]) using now twice the values of the (continuous) deformation gradient at the mid side points. </li>
1003
1004
</ol>
1005
1006
Details of the derivation of the EBST element are given below.
1007
1008
===5.1 Definition of the element geometry and computation of membrane strains===
1009
1010
As mentioned above a quadratic approximation of the geometry of the four elements patch is chosen using the position of the six nodes in the patch. It is useful to define the patch in the isoparametric space using the nodal positions given in the Table 2 (see also Figure 2).
1011
1012
1013
<div class="center" style="font-size: 75%;">
1014
'''Table 2'''. Isoparametric coordinates of the six nodes in the patch of Figure 2.</div>
1015
1016
{| class="wikitable" style="text-align: left; margin: 1em auto;"
1017
|- style="border-top: 2px solid;"
1018
|
1019
| 1 
1020
| 2 
1021
| 3 
1022
| 4 
1023
| 5 
1024
| 6
1025
|- style="border-top: 2px solid;"
1026
| <math display="inline">\xi </math> 
1027
| 0 
1028
| 1 
1029
| 0 
1030
| 1 
1031
| -1 
1032
| 1
1033
|- style="border-top: 2px solid;border-bottom: 2px solid;"
1034
| <math display="inline">\eta </math> 
1035
| 0 
1036
| 0 
1037
| 1 
1038
| 1 
1039
| 1 
1040
| -1
1041
1042
|}
1043
1044
The quadratic interpolation is defined by
1045
1046
<span id="eq-70"></span>
1047
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1048
|-
1049
| 
1050
{| style="text-align: left; margin:auto;" 
1051
|-
1052
| style="text-align: center;" | <math>{\boldsymbol \varphi }=\sum _{i=1}^{6}N_{i}{\boldsymbol \varphi }_{i}</math>
1053
|}
1054
| style="width: 5px;text-align: right;" | (70)
1055
|}
1056
1057
with (<math display="inline">\zeta=1-\xi-\eta</math>)
1058
1059
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1060
|-
1061
| 
1062
{| style="text-align: left; margin:auto;" 
1063
|-
1064
| style="text-align: center;" | <math>\begin{array}{ccc} N_{1}=\zeta{+\xi}\eta &  & N_{4}=\frac{\zeta }{2}\left( \zeta{-1}\right) \\ N_{2}=\xi{+\eta}\zeta &  & N_{5}=\frac{\xi }{2}\left( \xi{-1}\right) \\ N_{3}=\eta{+\zeta}\xi &  & N_{6}=\frac{\eta }{2}\left( \eta{-1}\right) \end{array} </math>
1065
|}
1066
| style="width: 5px;text-align: right;" | (71)
1067
|}
1068
1069
This interpolation allows to compute the displacement gradients at selected points in order to use an assumed strain approach. The computation of the gradients is performed at the mid side points of the central element of the patch denoted by <math display="inline">G_{1}</math>, <math display="inline">G_{2}</math> and <math display="inline">G_{3}</math> in Figure [[#img-2|2]]. This choice has the following advantages.
1070
1071
<div id='img-2'></div>
1072
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1073
|-
1074
|[[Image:draft_Samper_226033773-fig2.png|300px|Patch of elements in the isoparametric space.]]
1075
|- style="text-align: center; font-size: 75%;"
1076
| colspan="1" | '''Figure 2:''' Patch of elements in the isoparametric space.
1077
|}
1078
1079
* Gradients at the three mid side points depend only on the nodes belonging to the two elements adjacent to each side. This can be easily verified by sampling the derivatives of the shape functions at each mid-side point.
1080
1081
* When gradients are computed at the common mid-side point of two adjacent elements, the same values are obtained, as the coordinates of the same four points are used. This in practice means that the gradients at the mid-side points are independent of the element where they are computed. A side-oriented implementation of the finite element will therefore lead to a unique evaluation of the gradients per side.
1082
1083
The cartesian derivatives of the shape functions are computed at the original configuration by the standard expression
1084
1085
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1086
|-
1087
| 
1088
{| style="text-align: left; margin:auto;" 
1089
|-
1090
| style="text-align: center;" | <math>\left[ \begin{array}{c} N_{i,1}\\ N_{i,2} \end{array} \right] =\mathbf{J}^{-1}\left[ \begin{array}{c} N_{i,\xi } \\ N_{i,\eta } \end{array} \right] </math>
1091
|}
1092
| style="width: 5px;text-align: right;" | (72)
1093
|}
1094
1095
where the Jacobian matrix at the original configuration is
1096
1097
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1098
|-
1099
| 
1100
{| style="text-align: left; margin:auto;" 
1101
|-
1102
| style="text-align: center;" | <math>\mathbf{J=}\left[ \begin{array}{cc} \mathbf{\boldsymbol \varphi }_{^{\prime }\xi }^{0}\cdot \mathbf{t}_{1} & \mathbf{\boldsymbol \varphi  }_{^{\prime }\eta }^{0}\cdot \mathbf{t}_{1}\\ \mathbf{\boldsymbol \varphi }_{^{\prime }\xi }^{0}\cdot \mathbf{t}_{2} & \mathbf{\boldsymbol \varphi  }_{^{\prime }\eta }^{0}\cdot \mathbf{t}_{2} \end{array} \right] </math>
1103
|}
1104
| style="width: 5px;text-align: right;" | (73)
1105
|}
1106
1107
The deformation gradients on the middle surface, associated to an arbitrary spatial cartesian system and to the material cartesian system defined on the middle surface are related by
1108
1109
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1110
|-
1111
| 
1112
{| style="text-align: left; margin:auto;" 
1113
|-
1114
| style="text-align: center;" | <math>\left[ {\boldsymbol \varphi }_{^{\prime }1},\mathbf{\boldsymbol \varphi }_{^{\prime }2}\right] =\left[ \mathbf{\boldsymbol \varphi }_{^{\prime }\xi },\mathbf{\boldsymbol \varphi }_{^{\prime }\eta }\right]  \mathbf{J}^{-1} </math>
1115
|}
1116
| style="width: 5px;text-align: right;" | (74)
1117
|}
1118
1119
The Green-Lagrange membrane strains within the central triangle are now obtained using a linear assumed membrane strain field <math display="inline">\hat{\boldsymbol \varepsilon }_{m}</math>, i.e.
1120
1121
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1122
|-
1123
| 
1124
{| style="text-align: left; margin:auto;" 
1125
|-
1126
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_{m}=\hat{\boldsymbol \varepsilon }_{m} </math>
1127
|}
1128
| style="width: 5px;text-align: right;" | (75)
1129
|}
1130
1131
with
1132
1133
<span id="eq-76"></span>
1134
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1135
|-
1136
| 
1137
{| style="text-align: left; margin:auto;" 
1138
|-
1139
| style="text-align: center;" | <math>\hat{\boldsymbol \varepsilon }_{m}=(1-2\zeta ){\boldsymbol \varepsilon }_{m}^{1}+(1-2\xi ){\boldsymbol \varepsilon  }_{m}^{2}+(1-2\eta ){\boldsymbol \varepsilon }_{m}^{3}=\sum _{i=1}^{3}\bar{N}_{i} {\boldsymbol \varepsilon }_{m}^{i}</math>
1140
|}
1141
| style="width: 5px;text-align: right;" | (76)
1142
|}
1143
1144
where <math display="inline">{\boldsymbol \varepsilon }_{m}^{i}</math> are the membrane strains computed at the three mid side points <math display="inline">G_{i}</math> (<math display="inline">i=1,2,3</math>  see Figure [[#img-2|2]]). In Eq.([[#eq-76|76]]) <math display="inline">\bar{N}_{1}=(1-2\zeta )</math>, etc.
1145
1146
The gradient at each mid side point is computed from the quadratic interpolation ([[#eq-70|70]]):
1147
1148
<span id="eq-77"></span>
1149
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1150
|-
1151
| 
1152
{| style="text-align: left; margin:auto;" 
1153
|-
1154
| style="text-align: center;" | <math>\left( {\boldsymbol \varphi }_{^{\prime }\alpha }\right) _{G_{i}}={\boldsymbol \varphi }_{^{\prime  }\alpha }^{i}=\left[ \sum _{j=1}^{3}N_{j,\alpha }^{i}{\boldsymbol \varphi }_{j}\right] +N_{i+3,\alpha }^{i}{\boldsymbol \varphi }_{i+3}\quad ,\quad \alpha=1,2\quad ,\quad  i=1,2,3</math>
1155
|}
1156
| style="width: 5px;text-align: right;" | (77)
1157
|}
1158
1159
Substituting Eq.([[#eq-22|22]]) into ([[#eq-76|76]]) and using Eq.([[#eq-22|22]]) gives the membrane strain vector as
1160
1161
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1162
|-
1163
| 
1164
{| style="text-align: left; margin:auto;" 
1165
|-
1166
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_{m}=\sum _{i=1}^{3}\frac{1}{2}\bar{N}_{i}\left\{ \begin{array}{c} {\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}-1\\ {\boldsymbol \varphi }_{^{\prime }2}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}-1\\ 2{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i} \end{array} \right\} </math>
1167
|}
1168
| style="width: 5px;text-align: right;" | (78)
1169
|}
1170
1171
and the virtual membrane strains as
1172
1173
<span id="eq-79"></span>
1174
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1175
|-
1176
| 
1177
{| style="text-align: left; margin:auto;" 
1178
|-
1179
| style="text-align: center;" | <math>\delta{\boldsymbol \varepsilon }_{m}=\sum _{i=1}^{3}\bar{N}_{i}\left\{ \begin{array}{c} {\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \delta \mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}\\ {\boldsymbol \varphi }_{2}^{i}\cdot \delta \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}\\ \delta{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2} ^{i}+{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \delta \mathbf{\boldsymbol \varphi }_{2}^{i} \end{array} \right\} </math>
1180
|}
1181
| style="width: 5px;text-align: right;" | (79)
1182
|}
1183
1184
We note that the gradient at each mid side point <math display="inline">G_{i}</math> depends only on the coordinates of the three nodes of the central triangle and on those of an additional node in the patch, associated to the side <math display="inline">i</math> where the gradient is computed.
1185
1186
Combining Eqs.([[#eq-79|79]]) and ([[#eq-77|77]]) gives
1187
1188
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1189
|-
1190
| 
1191
{| style="text-align: left; margin:auto;" 
1192
|-
1193
| style="text-align: center;" | <math>\delta{\boldsymbol \varepsilon }_{m}=\mathbf{B}_{m}\delta \mathbf{a}^{p} </math>
1194
|}
1195
| style="width: 5px;text-align: right;" | (80)
1196
|}
1197
1198
where <math display="inline">\delta \mathbf{a}^{p}</math> is the patch displacement vector (see Eq.([[#eq-68|68]])) and <math display="inline">\mathbf{B}_{m}</math> is the membrane strain matrix. An explicit form of this matrix is given in the Appendix.
1199
1200
Differently from the original BST element the membrane strains within the EBST element are now a function of the displacements of the six patch nodes.
1201
1202
===5.2 Computation of curvatures===
1203
1204
The constant curvature field assumed for the BST element is chosen again here. The numerical evaluation of the line  integral in Eq.([[#eq-54|54]]) results in a sum over the integration points at the element boundary which are, in fact, the same points used for evaluating the gradients when computing the membrane strains. As one integration point is used over each side, it is not necessary to distinguish between sides (<math display="inline">i</math>) and integration points (<math display="inline">G_{i}</math>). In this way the curvatures can be computed by
1205
1206
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1207
|-
1208
| 
1209
{| style="text-align: left; margin:auto;" 
1210
|-
1211
| style="text-align: center;" | <math>{\boldsymbol \kappa }=2\sum _{i=1}^{3}\left[ \begin{array}{cc} L_{i,1}^M & 0\\ 0         & L_{i,2}^M \\ L_{i,2}^M & L_{i,1}^M  \end{array} \right] \left[ \begin{array}{c} \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }1}^{i}\\ \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }2}^{i} \end{array} \right] </math>
1212
|}
1213
| style="width: 5px;text-align: right;" | (81)
1214
|}
1215
1216
In the standard BST element <span id='citeF-21'></span><span id='citeF-24'></span>[[#cite-21|[21,24]]] the gradient <math display="inline">\mathbf{\boldsymbol \varphi  }_{\prime \alpha }^{i}</math> is computed as the average of the linear approximations over the two adjacent elements (see Section 4.3). In the enhanced version, the gradient is evaluated at each side <math display="inline">G_{i}</math> from the quadratic interpolation
1217
1218
<span id="eq-82"></span>
1219
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1220
|-
1221
| 
1222
{| style="text-align: left; margin:auto;" 
1223
|-
1224
| style="text-align: center;" | <math>\left[ \begin{array}{c} {\boldsymbol \varphi }_{\prime{1}}^{i}\\ {\boldsymbol \varphi }_{\prime{2}}^{i} \end{array} \right] =\left[ \begin{array}{cccc} N_{1,1}^{i} & N_{2,1}^{i} & N_{3,1}^{i} & N_{i+3,1}^{i}\\ N_{1,2}^{i} & N_{2,2}^{i} & N_{3,2}^{i} & N_{i+3,2}^{i} \end{array} \right] \left[ \begin{array}{c} {\boldsymbol \varphi }_{1}\\ {\boldsymbol \varphi }_{2}\\ {\boldsymbol \varphi }_{3}\\ {\boldsymbol \varphi }_{i+3} \end{array} \right]  </math>
1225
|}
1226
| style="width: 5px;text-align: right;" | (82)
1227
|}
1228
1229
Note again than at each side the gradients depend only on the positions of the three nodes of the central triangle and of an extra node (<math display="inline">i+3</math>), associated precisely to the side (<math display="inline">G_{i}</math>) where the gradient is computed.
1230
1231
Direction '''t'''<math display="inline">_{3}</math> in Eq.([[#eq-82|82]]) can be seen as a reference direction. If a different direction than that given by Eq.([[#eq-56|56]]) is chosen, at an angle <math display="inline">\theta </math> with the former, this has an influence of order <math display="inline">\theta ^{2}</math> in the projection. This justifies Eq.([[#eq-56|56]]) for the definition of '''t'''<math display="inline">_{3}</math> as a function exclusively of the three nodes of the central triangle, instead of using the 6-node isoparametric interpolation.
1232
1233
The variation of the curvatures can be obtained as
1234
1235
<span id="eq-89"></span>
1236
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1237
|-
1238
| 
1239
{| style="text-align: left; margin:auto;" 
1240
|-
1241
| style="text-align: right;" | <math>\delta{\boldsymbol \kappa } </math>
1242
| <math>  =2\sum _{i=1}^{3}\left[ \begin{array}{cc} L_{i,1}^M & 0\\ 0         & L_{i,2}^M\\ L_{i,2}^M & L_{i,1}^M \end{array} \right] \left\{ \sum _{i=1}^{3}\left[ \begin{array}{c} N_{j,1}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}_{j})\\ N_{j,2}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}_{j}) \end{array} \right] +\left[ \begin{array}{c} N_{i+3,1}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}^{i+3})\\ N_{i+3,2}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}^{i+3}) \end{array} \right] \right\} -</math>
1243
|-
1244
| style="text-align: right;" | 
1245
| <math>  -\sum _{i=1}^{3}\left[ \begin{array}{c} (L_{i,1}^M\rho _{11}^{1}+L_{i,2}^M\rho _{11}^{2})\\ (L_{i,1}^M\rho _{22}^{1}+L_{i,2}^M\rho _{22}^{2})\\ (L_{i,1}^M\rho _{12}^{1}+L_{i,2}^M\rho _{12}^{2}) \end{array} \right] (\mathbf{t}_{3}\cdot \delta \mathbf{u}_{i})=\mathbf{B}_{b} \delta \mathbf{a}^{p}</math>
1246
|}
1247
| style="width: 5px;text-align: right;" | (83)
1248
|}
1249
1250
where the definitions ([[#eq-61|61]]) and ([[#eq-63|63]]) still hold but with the new definition of <math display="inline">\mathbf{h}_{ij}</math> given by <span id='citeF-27'></span>[[#cite-27|[27]]]
1251
1252
<span id="eq-90"></span>
1253
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1254
|-
1255
| 
1256
{| style="text-align: left; margin:auto;" 
1257
|-
1258
| style="text-align: center;" | <math>\mathbf{h}_{ij}=\sum _{k=1}^{3}\left( L_{k,i}^M{\boldsymbol \varphi }_{\prime j}^{k} +L_{k,j}^M{\boldsymbol \varphi }_{^{\prime }i}^{k}\right) </math>
1259
|}
1260
| style="width: 5px;text-align: right;" | (84)
1261
|}
1262
1263
In Eq.([[#eq-89|83]])
1264
1265
<span id="eq-91"></span>
1266
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1267
|-
1268
| 
1269
{| style="text-align: left; margin:auto;" 
1270
|-
1271
| style="text-align: center;" | <math>\mathbf{B}_{b}=[\mathbf{B}_{b_{1}},\mathbf{B}_{b_{2}},\cdots ,\mathbf{B}_{b_{6}}]</math>
1272
|}
1273
| style="width: 5px;text-align: right;" | (85)
1274
|}
1275
1276
The expression of the curvature matrix <math display="inline">\mathbf{B}_b</math> is given in the Appendix. Details of the derivation of Eq.([[#eq-89|89]]) can be found in [27].
1277
1278
===5.3 The EBST1 element===
1279
1280
A simplified and yet very effective version of the EBST element can be obtained by using ''one point quadrature'' for the computation of all the element integrals. This element is termed EBST1. Note that this only affects the membrane stiffness matrices and it is equivalent to using a assumed constant membrane strain field defined by an average of the metric tensors computed at each side.
1281
1282
Numerical experiments have shown that both the EBST and the EBST1 elements are free of spurious energy modes.
1283
1284
==6 BOUNDARY CONDITIONS==
1285
1286
Elements at the domain boundary, where an adjacent element does not exist, deserve a special attention. The treatment of essential boundary conditions associated to translational constraints is straightforward, as they are the natural degrees of freedom of the element. The conditions associated to the normal vector are crucial in this formulation for bending. For clamped sides or symmetry planes, the normal vector <math display="inline">\mathbf{t}_{3}</math> must be kept fixed (clamped case), or constrained to move in the plane of symmetry (symmetry case). The former case can be seen as a special case of the latter, so we will consider symmetry planes only. This restriction can be imposed through the definition of the tangent plane at the boundary, including the normal to the plane of symmetry <math display="inline">\boldsymbol \varphi _{^{\prime }n}^{0}</math> that does not change during the process.
1287
1288
<div id='img-3'></div>
1289
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1290
|-
1291
|
1292
[[File:Draft_Samper_226033773_5749_Fig3.jpeg|600px|Local cartesian system for the treatment of symmetry boundary conditions]]
1293
|- style="text-align: center; font-size: 75%;"
1294
| colspan="1" | '''Figure 3:''' Local cartesian system for the treatment of symmetry boundary conditions
1295
|}
1296
1297
The tangent plane at the boundary (mid-side point) is expressed in terms of two orthogonal unit vectors referred to a local-to-the-boundary Cartesian system (see Figure [[#img-3|3]]) defined as
1298
1299
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1300
|-
1301
| 
1302
{| style="text-align: left; margin:auto;" 
1303
|-
1304
| style="text-align: center;" | <math>\left[\boldsymbol \varphi _{^{\prime }n}^{0},\;\bar{\boldsymbol \varphi }_{^{\prime }s}\right] </math>
1305
|}
1306
| style="width: 5px;text-align: right;" | (86)
1307
|}
1308
1309
where vector <math display="inline">\boldsymbol \varphi _{^{\prime }n}^{0}</math> is fixed during the process while direction <math display="inline">\bar{\boldsymbol \varphi }_{^{\prime }s}</math> emerges from the intersection of the symmetry plane with the plane defined by the central element (<math display="inline">M</math>). The plane (gradient) defined by the central element in the selected original convective Cartesian system (<math display="inline">\mathbf{t}_{1},\mathbf{t}_{2} </math>) is
1310
1311
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1312
|-
1313
| 
1314
{| style="text-align: left; margin:auto;" 
1315
|-
1316
| style="text-align: center;" | <math>\left[\boldsymbol \varphi _{^{\prime }1}^{M},\;\boldsymbol \varphi _{^{\prime  }2}^{M}\right] </math>
1317
|}
1318
| style="width: 5px;text-align: right;" | (87)
1319
|}
1320
1321
the intersection line (side <math display="inline">i</math>) of this plane with the plane of symmetry can be written in terms of the position of the nodes that define the side (<math display="inline">j </math> and <math display="inline">k</math>) and the original length of the side <math display="inline">l_{i}^{M}</math>, i.e.
1322
1323
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1324
|-
1325
| 
1326
{| style="text-align: left; margin:auto;" 
1327
|-
1328
| style="text-align: center;" | <math>\boldsymbol \varphi _{^{\prime }s}^{i}=\frac{1}{l_{i}^{M}}\left(\boldsymbol  \varphi _{k}-\boldsymbol \varphi _{j}\right) </math>
1329
|}
1330
| style="width: 5px;text-align: right;" | (88)
1331
|}
1332
1333
That together with the outer normal to the side <math display="inline">\mathbf{n}^{i} = \left[n_{1},n_{2}\right]^{T}=\left[\mathbf{n\cdot t}_{1},\mathbf{n\cdot t}_{2}\right]^{T}</math> (resolved in the selected original convective Cartesian system) leads to
1334
1335
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1336
|-
1337
| 
1338
{| style="text-align: left; margin:auto;" 
1339
|-
1340
| style="text-align: center;" | <math>\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }1}^{iT} \\ \boldsymbol \varphi _{^{\prime }2}^{iT} \end{array}\right]=\left[ \begin{array}{cc}n_{1} & -n_{2} \\ n_{2} & n_{1} \end{array}\right]\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }n}^{iT} \\ \boldsymbol \varphi _{^{\prime }s}^{iT} \end{array}\right] </math>
1341
|}
1342
| style="width: 5px;text-align: right;" | (89)
1343
|}
1344
1345
where, noting  that <math display="inline">\lambda </math> is the determinant of the gradient, the normal component of the gradient <math display="inline">\boldsymbol \varphi _{^{\prime }n}^{i}</math> can be approximated by
1346
1347
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1348
|-
1349
| 
1350
{| style="text-align: left; margin:auto;" 
1351
|-
1352
| style="text-align: center;" | <math>\boldsymbol \varphi _{^{\prime }n}^{i}=\frac{\boldsymbol \varphi _{^{\prime }n}^{0} }{\lambda |\boldsymbol \varphi _{^{\prime }s}^{i}|} </math>
1353
|}
1354
| style="width: 5px;text-align: right;" | (90)
1355
|}
1356
1357
In this way the contribution of the gradient at side <math display="inline">i</math> to vectors <math display="inline">\mathbf{ h}_{\alpha \beta }</math> (equations [[#eq-60|60]] and [[#eq-90|84]]) results in
1358
1359
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1360
|-
1361
| 
1362
{| style="text-align: left; margin:auto;" 
1363
|-
1364
| style="text-align: center;" | <math>\left[ \begin{array}{c}\mathbf{h}_{11}^{T} \\ \mathbf{h}_{22}^{T} \\ 2\mathbf{h}_{12}^{T} \end{array}\right]^{i}=2\left[ \begin{array}{cc}L_{i,1}^{M} & 0 \\ 0 & L_{i,2}^{M} \\ L_{i,2}^{M} & L_{i,1}^{M} \end{array}\right]\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }1}^{iT} \\ \boldsymbol \varphi _{^{\prime }2}^{iT} \end{array}\right]=2\left[ \begin{array}{cc}L_{i,1}^{M} & 0 \\ 0 & L_{i,2}^{M} \\ L_{i,2}^{M} & L_{i,1}^{M} \end{array}\right]\left[ \begin{array}{cc}n_{1} & -n_{2} \\ n_{2} & n_{1} \end{array}\right]\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }n}^{iT} \\ \boldsymbol \varphi _{^{\prime }s}^{iT} \end{array}\right] </math>
1365
|}
1366
| style="width: 5px;text-align: right;" | (91)
1367
|}
1368
1369
For the computation of the curvature variations, the contribution from the gradient at side <math display="inline">i</math> is now (see Ref. <span id='citeF-27'></span>[[#cite-27|[27]]])
1370
1371
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1372
|-
1373
| 
1374
{| style="text-align: left; margin:auto;" 
1375
|-
1376
| style="text-align: center;" | <math> \delta \left[ \begin{array}{c} \mathbf{h}_{11}^{T} \\ \mathbf{h}_{22}^{T} \\ 2\mathbf{h}_{12}^{T} \end{array} \right]^{i} =2\left[ \begin{array}{cc} L_{i,1}^{M} & 0 \\ 0 & L_{i,2}^{M} \\ L_{i,2}^{M} & L_{i,1}^{M} \end{array} \right]\left[ \begin{array}{cc} n_{1} & -n_{2} \\ n_{2} & n_{1} \end{array} \right]\left[ \begin{array}{c} \mathbf{0} \\ \frac{1}{L_{o}}\left[\delta \mathbf{u}_{k}-\delta \mathbf{u}_{j}\right]^{T} \end{array} \right]</math>
1377
|}
1378
| style="width: 5px;text-align: right;" |  (92a)
1379
|}
1380
1381
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1382
|-
1383
| 
1384
{| style="text-align: left; margin:auto;" 
1385
|-
1386
| style="text-align: center;" | <math> =\frac{2}{l_{i}^{M}}\left[ \begin{array}{c} -L_{i,1}^{M}n_{2} \\ L_{i,2}^{M}n_{1} \\ L_{i,1}^{M}n_{1}-L_{i,2}^{M}n_{2} \end{array} \right]\left[\delta \mathbf{u}_{k}-\delta \mathbf{u}_{j}\right]^{T}</math>
1387
|}
1388
| style="width: 5px;text-align: right;" |  (92b)
1389
|}
1390
1391
where the influence of variations in the length of vector <math display="inline">\boldsymbol \varphi  _{^{\prime }n}</math> has been neglected.
1392
1393
For a simple supported (hinged) side, the problem is not completely defined. The simplest choice is to neglect the contribution to the side rotations from the adjacent element missing in the patch in the evaluation of the curvatures via eq.([[#eq-54|54]]) <span id='citeF-29'></span>[[#cite-29|29]]. This is equivalent to assume that the gradient at the side is equal to the gradient in the central element, i.e.
1394
1395
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1396
|-
1397
| 
1398
{| style="text-align: left; margin:auto;" 
1399
|-
1400
| style="text-align: center;" | <math>\left[\boldsymbol \varphi _{^{\prime }1}^{i},\;\boldsymbol \varphi _{^{\prime }2}^{i} \right]=\left[\boldsymbol \varphi _{^{\prime }1}^{M},\;\boldsymbol \varphi  _{^{\prime }2}^{M}\right] </math>
1401
|}
1402
| style="width: 5px;text-align: right;" |  (93)
1403
|}
1404
1405
More precise changes can be however introduced to account for the different natural boundary conditions. One may assume that the curvature normal to the side is zero, and consider a contribution of the missing side to introduce this constraint. As the change of curvature parallel to the side is zero along the hinged side, both things lead to zero curvatures in both directions. Denoting the contribution to curvatures of the existing sides (<math display="inline">j </math> and <math display="inline">k</math>) by
1406
1407
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1408
|-
1409
| 
1410
{| style="text-align: left; margin:auto;" 
1411
|-
1412
| style="text-align: center;" | <math>\left[ \begin{array}{c}\kappa _{11} \\ \kappa _{22} \\ \kappa _{12} \end{array} \right]^{j-k} </math>
1413
|}
1414
|}
1415
1416
It can be easily shown that to set the normal curvature to zero the contribution of the simple supported side (<math display="inline">i</math>) should be
1417
1418
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1419
|-
1420
| 
1421
{| style="text-align: left; margin:auto;" 
1422
|-
1423
| style="text-align: center;" | <math>\left[ \begin{array}{c}\kappa _{11} \\ \kappa _{22} \\ \kappa _{12} \end{array} \right]^{i}=-\left[ \begin{array}{ccc}\left(n_{1}\right)^{4} & \left(n_{1}\right)^{2}\left(n_{2}\right)^{2} & \left(n_{1}\right)^{3}n_{2} \\ \left(n_{1}\right)^{2}\left(n_{2}\right)^{2} & \left(n_{2}\right)^{4} & n_{1}\left(n_{2}\right)^{3} \\ 2\left(n_{1}\right)^{3}n_{2} & 2n_{1}\left(n_{2}\right)^{3} & 2\left( n_{1}\right)^{2}\left(n_{2}\right)^{2} \end{array} \right]\left[ \begin{array}{c}\kappa _{11} \\ \kappa _{22} \\ \kappa _{12} \end{array} \right]^{j-k} </math>
1424
|}
1425
| style="width: 5px;text-align: right;" | (94)
1426
|}
1427
1428
For the case of a triangle with two sides associated to hinged sides, the normal curvatures to both sides must be zero. Denoting by <math display="inline">\mathbf{n}^{i}</math> and <math display="inline">\mathbf{n}^{j}</math> the normal to the sides, and by <math display="inline">\mathbf{m}^{i}</math> and <math display="inline"> \mathbf{m}^{j}</math> the dual base (associated to base <math display="inline">\mathbf{n}^{i}-</math> <math display="inline">\mathbf{ n}^{j}</math>), the contribution from the hinged sides (<math display="inline">i</math> and <math display="inline">j</math>) can be written as a function of the contribution of the only existing side (<math display="inline">k</math>):
1429
1430
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1431
|-
1432
| 
1433
{| style="text-align: left; margin:auto;" 
1434
|-
1435
| style="text-align: center;" | <math>\left[ \begin{array}{c}\kappa _{11} \\ \kappa _{22} \\ \kappa _{12} \end{array} \right]^{i-j}=-\left[ \begin{array}{c}m_{1}^{i}m_{1}^{j} \\ m_{2}^{i}m_{2}^{j} \\ m_{1}^{i}m_{2}^{j}+m_{2}^{i}m_{1}^{j} \end{array} \right]\left[ \begin{array}{ccc}2n_{1}^{i}n_{1}^{j} & 2n_{2}^{i}n_{2}^{j} & n_{1}^{i}n_{2}^{j}+n_{2}^{i}n_{1}^{j} \end{array} \right]\left[ \begin{array}{c}\kappa _{11} \\ \kappa _{22} \\ \kappa _{12} \end{array} \right]^{k} </math>
1436
|}
1437
| style="width: 5px;text-align: right;" | (95)
1438
|}
1439
1440
For a free edge the same approximation can be used but due to Poisson's effect this will lead to some error. The curvature variations of these contributions can be easily computed.
1441
1442
For the membrane formulation of element EBST, the gradient at the mid-side point of the boundary is assumed equal to the gradient of the main triangle.
1443
1444
==7 IMPLICIT SOLUTION SCHEME==
1445
1446
For a step <math display="inline">n</math> the configuration <math display="inline">\mathbf{\boldsymbol \varphi }^{n}</math> and the plastic strains <math display="inline">{\boldsymbol \varepsilon }_{p}^{n}</math> are known. The configuration <math display="inline">\mathbf{\boldsymbol \varphi }^{n}</math> is obtained by adding the total displacements to the original configuration <math display="inline"> \mathbf{\boldsymbol \varphi }^{n}=\mathbf{\boldsymbol \varphi }^{0} +\mathbf{u}^{n}</math>. The stresses are computed at each triangle using a single sampling (integration) point at the center and <math display="inline">N_{L}</math> integration points (layers) through the thickness. The plane stress state condition of the classical thin shell theory is assumed, so that for every layer three stress components are computed, (<math display="inline">\sigma _{11}</math>,<math display="inline">\sigma _{22}</math>, and <math display="inline">\sigma _{12}</math>) referred to the local cartesian system.
1447
1448
The computation of the incremental stresses is as follows:
1449
1450
<ol>
1451
1452
<li>Evaluate the incremental displacements: <math display="inline">\Delta \mathbf{u}^{n}=\mathbf{K}_{T}^{n}\mathbf{r}^{n}</math> where <math display="inline">\mathbf{K}_{T}</math> is the tangent stiffness matrix and '''r''' is the residual force vector  defined by for each element
1453
1454
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1455
|-
1456
| 
1457
{| style="text-align: left; margin:auto;" 
1458
|-
1459
| style="text-align: center;" | <math>\mathbf{r}^e_i =\int \int _A L_i {t}\, dA - \int \int _{A^\circ } ({B}_{m_i}^T {\boldsymbol \sigma }_m + {B}_{b_i}^T {\boldsymbol \sigma }_b)dA </math>
1460
|}
1461
| style="width: 5px;text-align: right;" | (96)
1462
|}</li>
1463
1464
The expression of the tangent stiffness matrix for the element is given below. Details of the derivation can be found in <span id='citeF-24'></span>[[#cite-24|[24]]],<span id='citeF-27'></span>[[#cite-27|[27]]].
1465
1466
<li>Generate the actual configuration <math display="inline">\mathbf{\boldsymbol \varphi }^{n+1} =\mathbf{\boldsymbol \varphi }^{n}+\Delta \mathbf{u}^{n}</math> </li>
1467
1468
<li>Compute the metric tensor <math display="inline">a_{\alpha \beta }^{n+1}\mathbf{ }</math>and the curvatures <math display="inline">\kappa _{\alpha \beta }^{n+1}</math> </li>
1469
1470
<li>Compute the total and elastic deformations at each layer <math display="inline">k</math> </li>
1471
1472
<span id="eq-96"></span>
1473
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1474
|-
1475
| 
1476
{| style="text-align: left; margin:auto;" 
1477
|-
1478
| style="text-align: right;" | <math>{\boldsymbol \varepsilon }_{k}^{n+1} </math>
1479
| <math>  ={\boldsymbol \varepsilon }_{m}^{n+1}+z_{k}{\boldsymbol \chi } ^{n+1}</math>
1480
| style="width: 5px;text-align: right;" | (97)
1481
|-
1482
| style="text-align: right;" | <math> \left[ {\boldsymbol \varepsilon }_{e}\right] _{k}^{n+1} </math>
1483
| <math>  ={\boldsymbol \varepsilon  }_{k}^{n+1}-\left[ {\boldsymbol \varepsilon }_{p}\right] _{k}^{n} </math>
1484
|}
1485
|}
1486
1487
<li>Compute the trial elastic stresses at each layer <math display="inline">k</math>
1488
1489
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1490
|-
1491
| 
1492
{| style="text-align: left; margin:auto;" 
1493
|-
1494
| style="text-align: center;" | <math>{\boldsymbol \sigma } _{k}^{n+1}=\mathbf{C}\left[ {\boldsymbol \varepsilon }_{e}\right] _{k}^{n+1} </math>
1495
|}
1496
| style="width: 5px;text-align: right;" | (98)
1497
|}</li>
1498
1499
<li>Check the plasticity condition and return to the plasticity surface. If necessary correct the plastic strains <math display="inline">\left[{\boldsymbol \varepsilon } _{p}\right] _{k}^{n+1}</math> at each layer (small strain plasticity) </li>
1500
1501
<li>Compute the generalized stresses
1502
1503
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1504
|-
1505
| 
1506
{| style="text-align: left; margin:auto;" 
1507
|-
1508
| style="text-align: right;" | <math>{\boldsymbol \sigma }^{n+1}_{m}  </math>
1509
| <math>  =\frac{h^{0}}{N_{L}}\sum _{k=1}^{N_{L}}\boldsymbol \sigma _{k} ^{n+1} w_{k}</math>
1510
|-
1511
| style="text-align: right;" | <math> {\boldsymbol \sigma }^{n+1}_{b}  </math>
1512
| <math>  =\frac{h^{0}}{N_{L}}\sum _{k=1}^{N_{L}}\boldsymbol \sigma _{k} ^{n+1}z_{k} w_{k} </math>
1513
|}
1514
| style="width: 5px;text-align: right;" | (99)
1515
|}</li>
1516
1517
Where <math display="inline"> w_{k}</math> is the weight of the through-the-thickness integration point. Recall that <math display="inline">z_{k}</math> is the current distance of the layer to the mid-surface and not the original distance. However, for small strain plasticity this distinction is not important.
1518
1519
This computation of stresses is adequate for an implicit scheme independent of the step size and it is exact for an elastic problem.
1520
1521
<li>Compute the residual force vector. The contribution for the <math display="inline">M</math>th element is given by
1522
1523
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1524
|-
1525
| 
1526
{| style="text-align: left; margin:auto;" 
1527
|-
1528
| style="text-align: center;" | <math>(\mathbf{r}^{M})^{n+1}=-A_{M}^{0}\left[ \begin{array}{cc} \mathbf{B}_{m}^{T} & \mathbf{B}_{b}^{T} \end{array} \right] ^{n+1}\left[ \begin{array}{c} \boldsymbol \sigma _{m}\\ \boldsymbol \sigma _{b} \end{array} \right] ^{n+1} </math>
1529
|}
1530
| style="width: 5px;text-align: right;" | (100)
1531
|}</li>
1532
1533
</ol>
1534
1535
===7.1 Tangent stiffness matrix===
1536
1537
As usual the tangent stiffness matrix is split into material and geometric components. The material tangent stiffness matrix is simply computed from the integral
1538
1539
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1540
|-
1541
| 
1542
{| style="text-align: left; margin:auto;" 
1543
|-
1544
| style="text-align: center;" | <math>\mathbf{K}^{M}=\int \int _{A^{M}}\mathbf{B}^{T}\mathbf{C}\mathbf{B}dA </math>
1545
|}
1546
| style="width: 5px;text-align: right;" | (101)
1547
|}
1548
1549
where <math display="inline">\mathbf{B}=\mathbf{B}_{m}+\mathbf{B}_{b}</math> includes:<br/>
1550
1551
* '''-'''  a membrane contribution <math display="inline">\mathbf{B}_{m}</math> given by Eq.([[#eq-51|51]]) or Eq.(80).
1552
1553
* '''-'''  a bending contribution <math display="inline">\mathbf{B}_{b}</math> given by Eq.([[#eq-69|69]]) or Eq.([[#eq-91|85]])  which is constant over the element.
1554
1555
<br/>
1556
1557
A three point quadrature is used for integrating the stiffness terms <math display="inline">\mathbf{B}_{m}^{T}\mathbf{C}\mathbf{B}_{m}</math> (recall that for the EBST element the membrane strains vary linearly within the element) whereas one point quadrature is chosen for the rest of the terms in <math display="inline">\mathbf{K}^{M}</math>.
1558
1559
===7.2 Geometric tangent stiffness matrix===
1560
1561
The geometric stiffness is written as
1562
1563
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1564
|-
1565
| 
1566
{| style="text-align: left; margin:auto;" 
1567
|-
1568
| style="text-align: center;" | <math>\mathbf{K}^{G}=\mathbf{K}_{m}^{G}+\mathbf{K}_{b}^{G} </math>
1569
|}
1570
| style="width: 5px;text-align: right;" | (102)
1571
|}
1572
1573
where subscripts <math display="inline">m</math> and <math display="inline">b</math> denote as usual membrane and bending contributions. For the BST element the membrane part is the same than for the standard constant strain triangle, leading to
1574
1575
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1576
|-
1577
| 
1578
{| style="text-align: left; margin:auto;" 
1579
|-
1580
| style="text-align: right;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{m}^{G}\mathbf{\;}\Delta \mathbf{u} </math>
1581
| <math> =A^{M}\sum _{i=1}^{3}\sum _{j=1}^{3}\left\{ \delta \mathbf{u}_{i}\;\left[ \begin{array}{cc} L_{i,1}^{M} & L_{i,2}^{M} \end{array} \right] \left[ \begin{array}{cc} N_{11} & N_{12}\\ N_{21} & N_{22} \end{array} \right] \left[ \begin{array}{c} L_{j,1}^{M}\\ L_{j,2}^{M} \end{array} \right] \Delta \mathbf{u}_{j}\right\} </math>
1582
|-
1583
| style="text-align: right;" | 
1584
| 
1585
|}
1586
| style="width: 5px;text-align: right;" | (103)
1587
|}
1588
1589
While for the EBST element the membrane part is computed as the sum of the contributions of the three sides, i.e.
1590
1591
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1592
|-
1593
| 
1594
{| style="text-align: left; margin:auto;" 
1595
|-
1596
| style="text-align: right;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{m}^{G}\mathbf{\;}\Delta \mathbf{u} </math>
1597
| <math> =\frac{A^{M}}{3}\sum _{k=1}^{3}\sum _{i=1}^{6}\sum _{j=1}^{6}\left\{ \delta \mathbf{u}_{i}\;\left[ \begin{array}{cc} N_{i,1}^{k} & N_{i,2}^{k} \end{array} \right] \left[ \begin{array}{cc} N_{11}^{k} & N_{12}^{k}\\ N_{21}^{k} & N_{22}^{k} \end{array} \right] \left[ \begin{array}{c} N_{j,1}^{k}\\ N_{j,2}^{k} \end{array} \right] \Delta \mathbf{u}_{j}\right\} </math>
1598
|-
1599
| style="text-align: right;" | 
1600
| 
1601
|}
1602
| style="width: 5px;text-align: right;" | (104)
1603
|}
1604
1605
where <math display="inline">N_{ij}={\sigma _{m}}_{ij}</math> are the axial forces defined in Eq.(29).
1606
1607
The geometric stiffness associated to bending moments is much more involved and can be found in  [27]. Numerical experiments have shown that the bending part of the geometric stiffness is not so important and can be disregarded in the iterative process.
1608
1609
Again three and one point quadratures are used for computing the membrane and bending contributions to the geometric stiffness matrix. We note that for elastic-plastic problems a uniform one point quadrature has been chosen for integration of both the membrane and bending stiffness matrices.
1610
1611
==8 EXPLICIT SOLUTION SCHEME==
1612
1613
For simulations including large non-linearities, such as frictional contact conditions on complex geometries or large instabilities in membranes, convergence is difficult to achieve with implicit schemes. In those cases an explicit solution algorithm is typically most advantageous. This scheme provides the solution for dynamic problems and also for static problems if an adequate damping is chosen.
1614
1615
The dynamic equations of motion to solve are of the form
1616
1617
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1618
|-
1619
| 
1620
{| style="text-align: left; margin:auto;" 
1621
|-
1622
| style="text-align: center;" | <math>\mathbf{r}(\mathbf{u}) + \mathbf{C} \dot{\mathbf{u}} + \mathbf{M} \ddot{\mathbf{u}} = 0 </math>
1623
|}
1624
| style="width: 5px;text-align: right;" | (108)
1625
|}
1626
1627
where <math display="inline">\mathbf{M}</math> is the mass matrix, <math display="inline">\mathbf{C}</math> is the damping matrix and the dot means the time derivative. The solution is performed using the ''central difference method''. To make the method competitive a diagonal (lumped) <math display="inline">\mathbf{M}</math> matrix is typically used and <math display="inline">\mathbf{C}</math> is taken proportional to <math display="inline">\mathbf{M}</math>. As usual, mass lumping is performed by assigning, one third of the triangular element mass to each node in the central element.
1628
1629
The explicit solution scheme can be summarized as follows. At each time step <math display="inline">n</math> where displacements have been computed:
1630
1631
<ol>
1632
1633
<li>Compute the internal forces <math display="inline">\mathbf{r}^{n}</math>. This simply follows the same steps (2-8) described for the implicit scheme in the previous section. </li>
1634
1635
<li>Compute the accelerations at time <math display="inline">t_{n}</math>
1636
1637
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1638
|-
1639
| 
1640
{| style="text-align: left; margin:auto;" 
1641
|-
1642
| style="text-align: center;" | <math>\ddot{\mathbf{u}}^{n} = {M}_d^{-1} [ \mathbf{r}^{n} - \mathbf{C} \dot{\mathbf{u}}^{n-1/2} ]  </math>
1643
|}
1644
|}</li>
1645
1646
where <math display="inline">{M}_d</math> is the diagonal (lumped) mass matrix.
1647
1648
<li>Compute the velocities at time <math display="inline">t_{n+1/2}</math>
1649
1650
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1651
|-
1652
| 
1653
{| style="text-align: left; margin:auto;" 
1654
|-
1655
| style="text-align: center;" | <math>\dot{\mathbf{u}}^{n+1/2} = \dot{\mathbf{u}}^{n-1/2} \ddot{\mathbf{u}}^{n} \delta t  </math>
1656
|}
1657
|}</li>
1658
1659
<li>Compute the displacements at  time <math display="inline">t_{n+1}</math>
1660
1661
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1662
|-
1663
| 
1664
{| style="text-align: left; margin:auto;" 
1665
|-
1666
| style="text-align: center;" | <math>\mathbf{u}^{n+1} = \mathbf{u}^{n} +\dot{\mathbf{u}}^{n+1/2} \delta t  </math>
1667
|}
1668
|}</li>
1669
<li>Update the shell geometry </li>
1670
<li>Check frictional contact conditions </li>
1671
1672
</ol>
1673
1674
Further details of the implementation of the standard BST element within an explicit solution scheme can be found in [26].
1675
1676
==9 EXAMPLES==
1677
1678
In this section several examples are presented to show the good performance of the rotation-free shell elements (BST, EBST and EBST1). The first five static examples are solved using an implicit code. The rest of the examples are computed using the explicit dynamic scheme. For the explicit scheme the  EBST element is always integrated using one integration point per element (EBST1 version) although not indicated.
1679
1680
===9.1 Patch tests===
1681
1682
The three elements considered (BST, EBST and EBST1) satisfy the membrane patch test defined in Figure [[#img-4|4]]. A uniform axial tensile stress is obtained in all cases.
1683
1684
<div id='img-4'></div>
1685
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1686
|-
1687
|
1688
[[File:Draft_Samper_226033773_1361_Fig4.jpeg|300px|Patch test for uniform tensile stress]]
1689
|- style="text-align: center; font-size: 75%;"
1690
| colspan="1" | '''Figure 4:''' Patch test for uniform tensile stress
1691
|}
1692
1693
<div id='img-5'></div>
1694
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1695
|-
1696
|
1697
[[File:Draft_Samper_226033773_8502_Fig5.jpeg|400px|Patch test for uniform torsion]]
1698
|- style="text-align: center; font-size: 75%;"
1699
| colspan="1" | '''Figure 5:''' Patch test for uniform torsion
1700
|}
1701
1702
The element bending formulation does not allow to apply external bending moments (there are not rotational DOFs). Hence it is not possible to analyse a patch of elements under loads leading to a uniform bending moment. A uniform torsion can be considered if a point load is applied at the corner of a rectangular plate with two consecutive free sides and two simple supported sides. Figure [[#img-5|5]] shows three patches leading to correct results both in displacements and stresses. All three patches are structured meshes. When the central node in the third patch is shifted from the center, the results obtained with the EBST and EBST1 elements are not correct. This however does not seems to preclude the excellent performance of these elements, as proved in the rest of the examples analyzed. On the other hand, the BST element  gives correct results in all torsion patch tests if natural boundary conditions are imposed in the formulation. If this is not the case, incorrect results are obtained even with structured meshes.
1703
1704
===9.2 Cook's membrane problem===
1705
1706
This example is used to assess the membrane performance of the EBST and EBST1 elements and to compare it with the standard linear triangle (constant strain) and the quadratic triangle (linear strain). This example involves important shear energy and was proposed to assess the distortion capability of elements. Figure [[#img-6|6]].a shows the geometry and the applied load. Figure [[#img-6|6]].b plots the vertical displacement of the upper vertex as a function of the number of nodes in the mesh. Results obtained with other isoparametric elements have also been  plotted for comparison. They include the constant strain triangle (CST), the bilinear quadrilateral (QUAD4) and the linear strain triangle (LST). Note that as this is a pure  membrane problem  the BST and the CST elements give identical results.
1707
1708
<div id='img-6'></div>
1709
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1710
|-
1711
|[[Image:draft_Samper_226033773-fig6a.png|300px|]]
1712
|
1713
[[File:Draft_Samper_226033773_1597_Fig6b.jpg|300px|Cook membrane problem (a) Geometry (b) Results]]
1714
|-
1715
|style="text-align: center; font-size: 75%;padding:10px;"|(a)
1716
|style="text-align: center; font-size: 75%;padding:10px;"|(b)
1717
|- style="text-align: center; font-size: 75%;"
1718
| colspan="2" | '''Figure 6:''' Cook membrane problem (a) Geometry (b) Results
1719
|}
1720
1721
From the plot shown it can be seen that the enhanced element with three integration points (EBST) gives values slightly better that the constant strain triangle for the most coarse mesh (only two elements). However, when the mesh is refined, a performance similar to the linear strain triangle is obtained that is dramatically superior that the former. On the other hand, if a one point quadrature is used (EBST1) the convergence in the reported displacement is notably better that for the rest of the elements.
1722
1723
===9.3 Cylindrical roof===
1724
1725
In this example an effective membrane interpolation is of primary importance. The geometry is a cylindrical roof supported by a rigid diaphragm at both ends and it is loaded by a uniform dead weight (see Figure [[#img-7|7]].a.). Only one quarter of the structure is modelled due to symmetry conditions. Unstructured and structured meshes are considered. In the latter case two orientations are possible (Figure [[#img-7|7]] shows orientation B).
1726
1727
Tables [[#table-3|3]], [[#table-4|4]] and [[#table-5|5]] present the normalized vertical displacements at the crown (point A) and at the midpoint of the free side (point B) for the two orientations of the structured meshes and for the non-structured mesh. Values used for normalization are <math display="inline">u_{A}=0.5407</math> y <math display="inline">u_{B}=-3.610</math> that are quoted in reference [32].
1728
1729
<div id='img-7'></div>
1730
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1731
|-
1732
|
1733
[[File:Draft_Samper_226033773_7308_Fig7a.jpeg|400px|]]
1734
|[[Image:draft_Samper_226033773-fig7b.png|400px|]]
1735
|- style="text-align: center; font-size: 75%;"
1736
| colspan="2" | '''Figure 7:''' Cylindrical roof under dead weight. <math>E=3 \times 10^{6}</math>, <math>\nu=0.0</math>, Thickness =3.0, shell weight =0.625 per unit area.
1737
|}
1738
1739
1740
{| class="wikitable" style="text-align: left; margin: 1em auto;"
1741
|+ <span id='table-3'></span>Table. 3 Cylindrical roof under dead weight. Normalized vertical displacements for mesh orientation A
1742
|- style="border-top: 2px solid;"
1743
| [c]|r|r|r|r|r|r|r| 
1744
| colspan='3' style="text-align: (null);border-left: 2px solid;border-right: 2px solid;" | Point-A
1745
| colspan='3' style="text-align: (null);border-left: 2px solid;border-right: 2px solid;" | Point-B
1746
|- style="border-top: 2px solid;"
1747
|  NDOFs 
1748
| EBST 
1749
| EBST1 
1750
| BST 
1751
| CBST 
1752
| EBST1 
1753
| BST
1754
|-
1755
| 16 
1756
| 0.65724 
1757
| 0.91855 
1758
| 0.74161 
1759
| 0.40950 
1760
| 0.70100 
1761
| 1.35230
1762
|-
1763
| 56 
1764
| 0.53790 
1765
| 1.03331 
1766
| 0.74006 
1767
| 0.54859 
1768
| 1.00759 
1769
| 0.75590
1770
|-
1771
| 208 
1772
| 0.89588 
1773
| 1.04374 
1774
| 0.88491 
1775
| 0.91612 
1776
| 1.02155 
1777
| 0.88269
1778
|-
1779
| 800 
1780
| 0.99658 
1781
| 1.01391 
1782
| 0.96521 
1783
| 0.99263 
1784
| 1.00607 
1785
| 0.96393
1786
|- style="border-bottom: 2px solid;"
1787
| 3136 
1788
| 1.00142 
1789
| 1.00385 
1790
| 0.99105 
1791
| 0.99881 
1792
| 1.00102 
1793
| 0.98992
1794
1795
|}
1796
1797
1798
{| class="wikitable" style="text-align: left; margin: 1em auto;"
1799
|+ <span id='table-4'></span>Table. 4 Cylindrical roof under dead weight. Normalized vertical displacements for mesh orientation B
1800
|- style="border-top: 2px solid;"
1801
| [c]|r|r|r|r|r|r|r| 
1802
| colspan='3' style="text-align: (null);border-left: 2px solid;border-right: 2px solid;" | Point-A
1803
| colspan='3' style="text-align: (null);border-left: 2px solid;border-right: 2px solid;" | Point-B
1804
|- style="border-top: 2px solid;"
1805
|  NDOFs 
1806
| EBST 
1807
| EBST1 
1808
| BST 
1809
| CBST 
1810
| EBST1 
1811
| BST
1812
|-
1813
| 16 
1814
| 0.26029 
1815
| 0.83917 
1816
| 0.40416 
1817
| 0.52601 
1818
| 0.86133 
1819
| 0.89778
1820
|-
1821
| 56 
1822
| 0.81274 
1823
| 1.10368 
1824
| 0.61642 
1825
| 0.67898 
1826
| 0.93931 
1827
| 0.68238
1828
|-
1829
| 208 
1830
| 0.97651 
1831
| 1.04256 
1832
| 0.85010 
1833
| 0.93704 
1834
| 1.00255 
1835
| 0.86366
1836
|-
1837
| 800 
1838
| 1.00085 
1839
| 1.01195 
1840
| 0.95626 
1841
| 0.99194 
1842
| 1.00211 
1843
| 0.95864
1844
|- style="border-bottom: 2px solid;"
1845
| 3136 
1846
| 1.00129 
1847
| 1.00337 
1848
| 0.98879 
1849
| 0.99828 
1850
| 1.00017 
1851
| 0.98848
1852
1853
|}
1854
1855
1856
{| class="wikitable" style="text-align: left; margin: 1em auto;"
1857
|+ <span id='table-5'></span>Table. 5 Cylindrical roof under dead weight. Normalized vertical displacements for non-structured mesh
1858
|- style="border-top: 2px solid;"
1859
| [c]|r|r|r|r|r|r|r| 
1860
| colspan='3' style="text-align: (null);border-left: 2px solid;border-right: 2px solid;" | Point-A
1861
| colspan='3' style="text-align: (null);border-left: 2px solid;border-right: 2px solid;" | Point-B
1862
|- style="border-top: 2px solid;"
1863
|  NDOFs 
1864
| EBST 
1865
| EBST1 
1866
| BST 
1867
| EBST 
1868
| EBST1 
1869
| BST
1870
|-
1871
| 851 
1872
| 0.97546 
1873
| 0.8581 
1874
| 0.97598 
1875
| 0.97662 
1876
| 1.0027 
1877
| 0.97194
1878
|-
1879
| 3311 
1880
| 0.98729 
1881
| 0.9682 
1882
| 0.98968 
1883
| 0.98476 
1884
| 1.0083 
1885
| 0.98598
1886
|- style="border-bottom: 2px solid;"
1887
| 13536 
1888
| 0.99582 
1889
| 0.9992 
1890
| 1.00057 
1891
| 0.99316 
1892
| 0.9973 
1893
| 0.99596
1894
1895
|}
1896
1897
Plots in Figure [[#img-7|7]].b show the normalized displacement of point-B for structured meshes as a function of the number of degrees of freedom for each case studied. An excellent convergence for the EBST element can be seen. The version with only one integration point (EBST1) presents a behavior a little more flexible and converges from above for structured meshes. Table [[#table-5|5]] shows that both the EBST and the EBST1 elements have an excellent behavior for non structured meshes.
1898
1899
===9.4 Open semi-spherical dome with point loads===
1900
1901
The main problem of finite elements with initially curved geometry is the so called membrane locking. The EBST element  has a quadratic interpolation of the geometry, then it may suffer from this problem. To assess this we resort to an example of inextensional bending. This is an hemispherical shell of radius <math display="inline">r=10</math> and thickness <math display="inline">h=0.04</math> with an 18<math display="inline">^{o}</math> hole in the pole and free at all boundaries, subjected to two inward and two outward forces 90<math display="inline">^{o}</math> apart. Material properties are <math display="inline">E=6.825\times{10}^{7}</math> and <math display="inline">\nu=0.3</math>. Figure [[#img-8|8]].a shows the discretized geometry (only one quarter of the geometry is considered due to symmetry).
1902
1903
<div id='img-8'></div>
1904
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1905
|-
1906
|
1907
[[File:Draft_Samper_226033773_7798_Fig8a.jpeg|400px|]]
1908
|[[Image:draft_Samper_226033773-fig8b.png|400px|Pinched hemispherical shell with a hole, (a)geometry, (b)normalized displacement]]
1909
|- style="text-align: center; font-size: 75%;"
1910
| colspan="2" | '''Figure 8:''' Pinched hemispherical shell with a hole, (a)geometry, (b)normalized displacement
1911
|}
1912
1913
In Figure [[#img-8|8]].b the displacements of the points under the loads have been plotted versus the number of nodes used in the discretization. Due to the orientation of the meshes chosen, the displacement of the point under the inward load is not the same as the displacement under the outward load, so in the figure an average (the absolute values) has been used. Results obtained with other elements have been included for comparison: two membrane locking free elements, namely the original linear BST element and a transverse shear-deformable quadrilateral (QUAD) [33]; a transverse shear deformable quadratic triangle (TRIA) [2] that is vulnerable to locking and an assumed strain quadratic triangle (TRIC) [3] that does not exhibit membrane locking.
1914
1915
From the plotted results it can be seen that the EBST element presents slight membrane locking in bending dominated problems with initially curved geometries. This locking is much less severe than in a standard quadratic triangle. Membrane locking disappears when only one integration point is used (EBST1 element).
1916
1917
===9.5 Inflation of a sphere===
1918
1919
The example is the inflation of a spherical shell under internal pressure. An incompressible Mooney-Rivlin constitutive material has been considered. The Ogden parameters are <math display="inline">N=2</math>, <math display="inline">\alpha _{1}=2</math>, <math display="inline">\mu _{1}=40</math>, <math display="inline">\alpha _{2}=-2</math>, <math display="inline">\mu _{2}=-20</math>. Due to the simple geometry an analytical solution exists [34] (with <math display="inline">\gamma =R/R^{0}</math>):
1920
1921
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1922
|-
1923
| 
1924
{| style="text-align: left; margin:auto;" 
1925
|-
1926
| style="text-align: center;" | <math> p=\frac{h^{0}}{R^{0}\gamma ^{2}}\frac{dW}{d\gamma }=\frac{8h^{0} }{R^{0}\gamma ^{2}} \left( \gamma ^{6}-1\right) \left( \mu _{1}-\mu _{2}\gamma ^{2}\right) </math>
1927
|}
1928
|}
1929
1930
In this numerical simulation the same geometric and material parameters used in Ref. <span id='citeF-23'></span>[[#cite-23|[23]]] have been adopted: <math display="inline">R^{0}=1</math> and <math display="inline">h^{0}=0.02</math>. The three meshes of EBST1 element considered to evaluate convergence are shown in Figure [[#img-9|9]].a. The value of the actual radius as a function of the internal pressure is plotted in Figure [[#img-9|9]].b for the different meshes and is also compared with the analytical solution. It can be seen that with a few degrees of freedom it is possible to obtain an excellent agreement for the range of strains considered. The final value corresponds to a  ratio of <math display="inline">h/R=0.00024</math>.
1931
1932
<div id='img-9'></div>
1933
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1934
|-
1935
|[[File:Draft_Samper_226033773_5233_Fig9.jpeg|500px|]]
1936
|-
1937
|[[File:Draft_Samper_226033773_9697_Fig9b.jpeg|300px|Inflation of sphere of Mooney-Rivlin material. (a) Meshes of EBST1 elements used in the analysis (b) Change of radius as a function of the internal pressure.]]
1938
|- style="text-align: center; font-size: 75%;"
1939
| '''Figure 9:''' Inflation of sphere of Mooney-Rivlin material. (a) Meshes of EBST1 elements used in the analysis (b) Change of radius as a function of the internal pressure.
1940
|}
1941
1942
===9.6 Clamped spherical dome under impulse pressure loading===
1943
1944
The geometry of the dome and the material properties chosen are shown in Figure [[#img-10|10]]. A uniform pressure load of 600 psi is applied to the upper surface of the dome. The different meshes used in the analysis are shown in Figure [[#img-11|11]]. One fourth of the dome is considered only due to symmetry. Two different analyses under elastic and elastic-plastic conditions were carried out. The number of thickness layers in eq.([[#eq-96|96]]) is four. Numerical experiments show that this suffice to provide an accurate solution for large elastic-plastic problems [26]. Results are obtained using the explicit scheme.
1945
1946
<div id='img-10'></div>
1947
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1948
|-
1949
|[[Image:draft_Samper_226033773-fig10.png|400px|Spherical dome under impulse pressure. Geometry and material]]
1950
|- style="text-align: center; font-size: 75%;"
1951
| colspan="1" | '''Figure 10:''' Spherical dome under impulse pressure. Geometry and material
1952
|}
1953
1954
<div id='img-11'></div>
1955
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1956
|-
1957
|
1958
[[File:Draft_Samper_226033773_1877_Fig11.jpeg|600px|Spherical dome under impulse pressure. Meshes used in the analysis. Mesh-1 with 338 elements, Mesh-2 with 1250 elements, and Mesh-3 with 2888 elements]]
1959
|- style="text-align: center; font-size: 75%;"
1960
| colspan="1" | '''Figure 11:''' Spherical dome under impulse pressure. Meshes used in the analysis. Mesh-1 with 338 elements, Mesh-2 with 1250 elements, and Mesh-3 with 2888 elements
1961
|}
1962
1963
<div id='img-12'></div>
1964
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1965
|-
1966
|[[Image:draft_Samper_226033773-fig12.png|600px|Spherical dome under impulse pressure. History of central deflection for elastic material]]
1967
|- style="text-align: center; font-size: 75%;"
1968
| colspan="1" | '''Figure 12:''' Spherical dome under impulse pressure. History of central deflection for elastic material
1969
|}
1970
1971
<div id='img-13'></div>
1972
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1973
|-
1974
|[[Image:draft_Samper_226033773-fig13.png|600px|Spherical dome under impulse pressure. History of central deflection for elastic-plastic material]]
1975
|- style="text-align: center; font-size: 75%;"
1976
| colspan="1" | '''Figure 13:''' Spherical dome under impulse pressure. History of central deflection for elastic-plastic material
1977
|}
1978
1979
Figure [[#img-12|12]] shows results for the time history of the central deflection using different meshes and ''elastic material properties'' for both  BST and EBST1 elements. Results are almost identical for mesh-2 and mesh-3, showing the excellent convergence properties. The coarsest mesh shows some differences between both elements, but for the finer meshes the results are almost identical. Figure [[#img-13|13]] shows similar results but now for an ''elastic-plastic material''. The excellent convergence  of the BST and EBST elements is again noticeable.
1980
1981
Results obtained with the present elements compare very well with published results using fine meshes. See for example ABAQUS Explicit example problems manual <span id='citeF-35'></span>[[#cite-35|[35]]] and WHAMS-3D manual [36], showing plots comparing results using different shell elements.
1982
1983
A summary of results for the central deflection at significant times is given in Tables [[#table-6|6]] and [[#table-7|7]]. Further details on the solution of this problem with the standard  BST element can be found in [26].
1984
1985
1986
{| class="wikitable" style="text-align: right; margin: 1em auto;"
1987
|+ <span id='table-6'></span>Table. 6 Spherical dome. Elastic material. Comparison of the central deflection values at the mid point obtained with the BST and EBST1  elements for different meshes
1988
|- style="border-top: 2px solid;"
1989
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Element/mesh 
1990
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t = 0.2 ms</math>
1991
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t = 0.4 ms</math>
1992
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t = 0.6 ms</math>
1993
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t = 0.8 ms</math>
1994
|- style="border-top: 2px solid;"
1995
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |   BST Coarse  
1996
| style="border-left: 2px solid;border-right: 2px solid;" | -0.05155 
1997
| style="border-left: 2px solid;border-right: 2px solid;" | -0.09130 
1998
| style="border-left: 2px solid;border-right: 2px solid;" | 0.04414 
1999
| style="border-left: 2px solid;border-right: 2px solid;" | -0.08945 
2000
|-
2001
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | BST Medium  
2002
| style="border-left: 2px solid;border-right: 2px solid;" | -0.04542 
2003
| style="border-left: 2px solid;border-right: 2px solid;" | -0.09177 
2004
| style="border-left: 2px solid;border-right: 2px solid;" | 0.03863 
2005
| style="border-left: 2px solid;border-right: 2px solid;" | -0.08052 
2006
|-
2007
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | BST Fine    
2008
| style="border-left: 2px solid;border-right: 2px solid;" | -0.04460 
2009
| style="border-left: 2px solid;border-right: 2px solid;" | -0.09022 
2010
| style="border-left: 2px solid;border-right: 2px solid;" | 0.03514 
2011
| style="border-left: 2px solid;border-right: 2px solid;" | -0.08132 
2012
|- style="border-top: 2px solid;"
2013
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  EBST1 Coarse  
2014
| style="border-left: 2px solid;border-right: 2px solid;" | -0.05088 
2015
| style="border-left: 2px solid;border-right: 2px solid;" | -0.08929 
2016
| style="border-left: 2px solid;border-right: 2px solid;" | 0.04348 
2017
| style="border-left: 2px solid;border-right: 2px solid;" | -0.08708 
2018
|-
2019
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | EBST1 Medium  
2020
| style="border-left: 2px solid;border-right: 2px solid;" | -0.04527 
2021
| style="border-left: 2px solid;border-right: 2px solid;" | -0.09134 
2022
| style="border-left: 2px solid;border-right: 2px solid;" | 0.03865 
2023
| style="border-left: 2px solid;border-right: 2px solid;" | -0.07979 
2024
|- style="border-bottom: 2px solid;"
2025
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | EBST1 Fine    
2026
| style="border-left: 2px solid;border-right: 2px solid;" | -0.04453 
2027
| style="border-left: 2px solid;border-right: 2px solid;" | -0.09004 
2028
| style="border-left: 2px solid;border-right: 2px solid;" | 0.03510 
2029
| style="border-left: 2px solid;border-right: 2px solid;" | -0.08099 
2030
2031
|}
2032
2033
2034
{| class="wikitable" style="text-align: left; margin: 1em auto;"
2035
|+ <span id='table-7'></span>Table. 7 Spherical dome. Elastic-plastic material. Comparison of the central deflection values at the mid point obtained with the BST and EBST1  elements for different meshes
2036
|- style="border-top: 2px solid;"
2037
| [c]|l|r|r|r|r| Element/mesh 
2038
| <math>t = 0.2 ms</math>
2039
| <math>t = 0.4 ms</math>
2040
| <math>t = 0.6 ms</math>
2041
| <math>t = 0.8 ms</math>
2042
|- style="border-top: 2px solid;"
2043
|   BST Coarse  
2044
| -0.05888 
2045
| -0.05869 
2046
| -0.02938 
2047
| -0.06521 
2048
|-
2049
| BST Medium  
2050
| -0.05376 
2051
| -0.06000 
2052
| -0.02564 
2053
| -0.06098 
2054
|-
2055
| BST Fine    
2056
| -0.05312 
2057
| -0.05993 
2058
| -0.02464 
2059
| -0.06105 
2060
|- style="border-top: 2px solid;"
2061
|  EBST1 Coarse  
2062
| -0.05827 
2063
| -0.05478 
2064
| -0.02792 
2065
| -0.06187 
2066
|-
2067
| EBST1 Medium  
2068
| -0.05374 
2069
| -0.05884 
2070
| -0.02543 
2071
| -0.06080 
2072
|- style="border-bottom: 2px solid;"
2073
| EBST1 Fine    
2074
| -0.05317 
2075
| -0.05935 
2076
| -0.02458 
2077
| -0.06085 
2078
2079
|}
2080
2081
===9.7 Cylindrical panel under impulse loading===
2082
2083
The geometry of the cylinder and the material properties are shown in Figure [[#img-14|14]]. A prescribed initial normal velocity of <math display="inline">v_{o}=-5650</math> in/sec is applied to the points in the region shown modelling the effect of the detonation of an explosive layer. The panel is assumed clamped along all the boundary. One half of the cylinder is discretized only due to symmetry conditions. Three different meshes of <math display="inline">6\times{12}</math>, <math display="inline">12\times{32}</math> and <math display="inline">18\times{48}</math>  BST elements were used for the analysis. The deformed configurations for <math display="inline">time =1 msec</math> are shown for the three meshes in Figure [[#img-15|15]].
2084
2085
<div id='img-14'></div>
2086
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2087
|-
2088
|[[Image:draft_Samper_226033773-fig14.png|600px|Cylindrical panel under impulse loading. Geometry and material properties]]
2089
|- style="text-align: center; font-size: 75%;"
2090
| colspan="1" | '''Figure 14:''' Cylindrical panel under impulse loading. Geometry and material properties
2091
|}
2092
2093
<div id='img-15'></div>
2094
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2095
|-
2096
|[[Image:draft_Samper_226033773-fig15.png|600px|Impulsively loaded cylindrical panel. Deformed meshes for time =1 msec]]
2097
|- style="text-align: center; font-size: 75%;"
2098
| colspan="1" | '''Figure 15:''' Impulsively loaded cylindrical panel. Deformed meshes for <math>time =1 msec</math>
2099
|}
2100
2101
<div id='img-16'></div>
2102
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2103
|-
2104
|[[Image:draft_Samper_226033773-fig16.png|600px|Cylindrical panel under impulse loading. Time evolution of the displacement of two points along the crown line. Comparison of results obtained with BST and EBST1 elements (mesh 1: 6×12 elements and mesh 3: 18×48 elements) and experimental values ]]
2105
|- style="text-align: center; font-size: 75%;"
2106
| colspan="1" | '''Figure 16:''' Cylindrical panel under impulse loading. Time evolution of the displacement of two points along the crown line. Comparison of results obtained with BST and EBST1 elements (mesh 1: <math>6\times{12}</math> elements and mesh 3: <math>18\times{48}</math> elements) and experimental values 
2107
|}
2108
2109
The analysis was performed assuming an elastic-perfect plastic material behaviour (<math display="inline">\sigma _y = k</math> <math display="inline">k'=0</math>). A study of the convergence of the solution with the number of thickness layers showed again that four layers suffice to capture accurately the non linear material response [26].
2110
2111
A comparison of the results obtained with both elements using the coarse mesh and the finer mesh is shown in Figure [[#img-16|16]] where experimental results reported in <span id='citeF-37'></span>[[#cite-37|[37]]] have also been plotted for comparison purposes. Good agreement between the numerical and experimental results is obtained. Figures [[#img-16|16]] show the time evolution of the vertical displacement of two reference points along the center line located at <math display="inline">y=6.28</math>in and <math display="inline">y=9.42</math>in, respectively. For the finer mesh results between both elements are almost identical. For the coarse mesh it can been seen again that the element BST is more flexible than element EBST1.
2112
2113
The numerical values of the vertical displacement at the two reference points obtained with the BST and EBST1  elements after a time of 0.4 ms using the <math display="inline">16\times{32}</math> mesh are compared in Table [[#table-8|8]]  with a numerical solution obtained by Stolarski ''et al.'' [38] using a curved triangular shell element and the <math display="inline">16\times{32}</math> mesh. Experimental results reported in [37] are also given for comparison. It is interesting to note the reasonable agreement of the results for <math display="inline">y=6.28</math>in. and the discrepancy of present and other published numerical solutions with the experimental value for <math display="inline">y=9.42</math>in.
2114
2115
2116
{| class="wikitable" style="text-align: left; margin: 1em auto;"
2117
|+ <span id='table-8'></span>Table. 8 Cylindrical panel under impulse load. Comparison of vertical displacement values of two central points for <math>t=0.4</math> ms
2118
|- style="border-top: 2px solid;"
2119
| [c]|l|r|r|r|r| 
2120
| colspan='2' style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Vertical displacement (in.)
2121
|- style="border-top: 2px solid;"
2122
|  element/mesh                
2123
| <math>y=6.28</math>in 
2124
| <math>y=9.42</math>in 
2125
|- style="border-top: 2px solid;"
2126
|  BST  (<math display="inline"> 6\times 12</math> el.)    
2127
| -1.310     
2128
| -0.679      
2129
|-
2130
| BST  (<math display="inline">18\times 48</math> el.)    
2131
| -1.181     
2132
| -0.587      
2133
|-
2134
| EBST1 (<math display="inline"> 6\times 12</math> el.)    
2135
| -1.147     
2136
| -0.575      
2137
|-
2138
| EBST1 (<math display="inline">18\times 48</math> el.)    
2139
| -1.171     
2140
| -0.584      
2141
|-
2142
| Stolarski ''et al.'' [38] 
2143
| -1.183     
2144
| -0.530      
2145
|- style="border-bottom: 2px solid;"
2146
| Experimental [37] 
2147
| -1.280     
2148
| -0.700      
2149
2150
|}
2151
2152
<div id='img-17'></div>
2153
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2154
|-
2155
|[[Image:draft_Samper_226033773-fig17.png|600px|Cylindrical panel under impulse loading. Final deformation (t=1 msec) of the panel at the cross section y=6.28 in Comparison with experimental values. ]]
2156
|- style="text-align: center; font-size: 75%;"
2157
| colspan="1" | '''Figure 17:''' Cylindrical panel under impulse loading. Final deformation (<math>t=1 msec</math>) of the panel at the cross section <math>y=6.28 in</math> Comparison with experimental values. 
2158
|}
2159
2160
<div id='img-18'></div>
2161
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2162
|-
2163
|[[Image:draft_Samper_226033773-fig18.png|600px|Cylindrical panel under impulse loading. Final deformation (t=1 msec) of the panel at the crown line (x=0.00 in). Comparison with experimental values. ]]
2164
|- style="text-align: center; font-size: 75%;"
2165
| colspan="1" | '''Figure 18:''' Cylindrical panel under impulse loading. Final deformation (<math>t=1 msec</math>) of the panel at the crown line (<math>x=0.00 in</math>). Comparison with experimental values. 
2166
|}
2167
2168
The deformed shapes of the transverse section for <math display="inline">y=6.28</math>in. and the longitudinal section for <math display="inline">x=0</math> obtained with the both elements for the coarse and the fine meshes after 1ms. are compared with the experimental results in Figures [[#img-17|17]] and [[#img-18|18]].  Excellent agreement is observed for the fine mesh for both elements.
2169
2170
===9.8 Airbag Membranes===
2171
2172
===Inflation/deflation of a circular airbag===
2173
2174
This example has been taken from Ref.[23] where it is shown that the final configuration is mesh dependent due to the strong instabilities leading to a non-uniqueness of the solution. In [23]  it is also discussed the important regularizing properties of the bending energy, that when disregarded leads to massive wrinkling in the compressed zones.
2175
2176
The airbag geometry is initially circular with an undeformed radius of <math display="inline">0.35</math>.  The constitutive material is a linear isotropic elastic one with modulus of elasticity <math display="inline">E=6\times 10^{7}</math>Pa, Poisson's ratio <math display="inline">\nu =0.3</math> and density <math display="inline">\rho = 2000</math>kg/m<math display="inline">^3</math>.  Arbitrarily only one quarter of the geometry has been modelled.  Only the normal displacement to the original plane is constrained along the boundaries.  The thickness considered is <math display="inline">h=0.0004</math>m and the inflation pressure is <math display="inline">5000</math>Pa. Pressure is linearly increased from 0 to the final value in <math display="inline">t=0.15</math> sec.
2177
2178
Figure 19 shows the final deformed configurations for a mesh with 10201 nodes and 20000 EBST1 elements.  The figure on the left (a) corresponds to an analysis including full bending effects and the right figure (b) is a pure membrane analysis.
2179
2180
We note that when the bending energy is included a more regular final pattern is obtained.  Also the final pattern is rather independent of the discretization (note that the solution is non unique due to the strong instabilities), and a massive wrinkling appears in the center of the modelled region.  On the other hand, the pure membrane solution shows a wrinkling pattern where the width of the wrinkle is the length of the element.
2181
2182
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2183
|-
2184
|
2185
[[File:Draft_Samper_226033773_2615_Fig19.jpeg|600px|Inflation of a circular airbag. Deformed configurations for final pressure. (a) bending formulation; (b) membrane formulation.]]
2186
|- style="text-align: center; font-size: 75%;"
2187
| colspan="1" | '''Figure 19:''' Inflation of a circular airbag. Deformed configurations for final pressure. (a) bending formulation; (b) membrane formulation.
2188
|}
2189
2190
<div id='img-20'></div>
2191
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2192
|-
2193
|
2194
[[File:Draft_Samper_226033773_9407_Fig20.jpeg|600px|Inflation and deflation of a circular air-bag.]]
2195
|- style="text-align: center; font-size: 75%;"
2196
| colspan="1" | '''Figure 20:''' Inflation and deflation of a circular air-bag.
2197
|}
2198
2199
<div id='img-21'></div>
2200
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2201
|-
2202
|
2203
[[File:Draft_Samper_226033773_9850_Fig21.jpeg|600px|Inflation and deflation of a closed  tube. L=5, D=1, h=5×10⁻³.]]
2204
|- style="text-align: center; font-size: 75%;"
2205
| colspan="1" | '''Figure 21:''' Inflation and deflation of a closed  tube. <math>L=5</math>, <math>D=1</math>, <math>h=5\times 10^{-3}</math>.
2206
|}
2207
2208
<div id='img-22'></div>
2209
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2210
|-
2211
|
2212
[[File:Draft_Samper_226033773_8861_Fig22.jpeg|600px|Inflation and deflation of a closed  tube. L=6, D=1, h=3×10⁻³.]]
2213
|- style="text-align: center; font-size: 75%;"
2214
| colspan="1" | '''Figure 22:''' Inflation and deflation of a closed  tube. <math>L=6</math>, <math>D=1</math>, <math>h=3\times 10^{-3}</math>.
2215
|}
2216
2217
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2218
|-
2219
|
2220
[[File:Draft_Samper_226033773_6495_fig23.jpg|600px|Inflation of a square airbag against an spherical object. Deformed configurations for different times. Left figure: results obtained with the full bending formulation. Right figure: results obtained with a pure membrane solution.]]
2221
|- style="text-align: center; font-size: 75%;"
2222
| colspan="1" | '''Figure 23:''' Inflation of a square airbag against an spherical object. Deformed configurations for different times. Left figure: results obtained with the full bending formulation. Right figure: results obtained with a pure membrane solution.
2223
|}
2224
2225
Figure 20 shows the results obtained for the de-inflation process.  Note that the spherical membrane falls down due to the self weight.  The final configuration is of course non-unique.
2226
2227
The next problem is the study of the inflating and de-inflating of a tube with a semi-spherical end cap.  The tube diameter is <math display="inline">D=1</math>, its total length is <math display="inline">L=5</math>m and the thickness <math display="inline">h=5\times 10^{-3}</math>m.  The material has the following properties <math display="inline">E=4\times 10^{8}</math>Pa, <math display="inline">\nu =0.35 </math>, <math display="inline">\rho =5\times 10^{4}</math>kg/m<math display="inline">^3</math>.  The tube is inflated fast until a pressure of <math display="inline">10^4</math> and then is de-inflated under self weight.  The analysis is performed with a mesh of 4176 EBST1 elements and 2163 nodes modelling a quarter of the geometry.  The evolution of the tube walls during the de-inflating process can be seen in Figure 21.  Note that the central part collapses as expected, while the semi-spherical cap remains unaltered.
2228
2229
The same analysis is repeated for a longer and thinner tube (<math display="inline">L=6</math>m and <math display="inline">h=3\times 10^{-3}</math>m).  The same material than in the previous case was chosen with a higher density (<math display="inline">\rho =7.5\times 10^{4}</math>kg/m<math display="inline">^3</math>).  The evolution of the tube walls is shown in Figure 22.  Note that the central part collapses again but in a less smoother manner due to the smaller thickness.
2230
2231
The last example of this kind is the inflation of a square airbag supporting a spherical object.  This example resembles a problem studied (numerically and experimentally) in Ref.[39], where fluid-structure interaction is the main subject.  Here the fluid is not modelled, and a uniform pressure is applied over all the internal surfaces.  The lower surface part of the airbag is limited by a rigid plane and on the upper part a spherical dummy object is set to study the interaction between the airbag and the object.
2232
2233
The airbag geometry is initially square with an undeformed side length of 0.643m.  The constitutive material chosen is a linear isotropic elastic one with <math display="inline">E=5.88\times 10^8</math>Pa, <math display="inline">\nu =0.4</math> and a density of <math display="inline">\rho = 1000</math> kg/m<math display="inline">^3</math>.  Only one quarter of the geometry has been modelled due to symmetry.  The thickness <math display="inline">h=0.00075</math>m and the inflation pressure is 250000Pa.  Pressure is linearly incremented from 0 to the final value in <math display="inline">t=0.15</math>sec.  The spherical object has a radius <math display="inline">r=0.08</math>m and a mass of 4.8kg (one quarter), and is subjected to gravity load during all the process.
2234
2235
The mesh includes 8192 EBST1 elements and 4225 nodes on each surface of the airbag.  Figure 23 shows the deformed configurations for three different times.  The sequence on the left of the figure corresponds to an analysis including full bending effects and the sequence on the right is the result of a pure membrane analysis.  A standard penalty formulation is used for frictionless contact.
2236
2237
===9.9 S-rail sheet stamping===
2238
2239
The final problem corresponds to one of the sheet stamping benchmark tests proposed in NUMISHEET'96 <span id='citeF-40'></span>[[#cite-40|[40]]].  The analysis comprises two parts, namely, stamping of a S-rail sheet component and springback computations once the stamping tools are removed.  Figure [[#img-24|24]] shows the deformed sheet after springback.
2240
2241
<div id='img-24'></div>
2242
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2243
|-
2244
|[[Image:draft_Samper_226033773-fig_n1.png|600px|Stamping of a S-rail. Final deformation of the sheet after springback obtained in the simulation. The triangular mesh of the deformed sheet is also shown]]
2245
|- style="text-align: center; font-size: 75%;"
2246
| colspan="1" | '''Figure 24:''' Stamping of a S-rail. Final deformation of the sheet after springback obtained in the simulation. The triangular mesh of the deformed sheet is also shown
2247
|}
2248
2249
The detailed geometry and material data can be found in the proceedings of the conference <span id='citeF-40'></span>[[#cite-40|[40]]] or in the web <span id='citeF-41'></span>[[#cite-41|[41]]]. The mesh used for the sheet has 6000 three  node triangular elements and 3111 points (Figure 24). The tools are treated as rigid bodies. The meshes used for the sheet and the tools are those provided by the  benchmark organizers. The material considered here is a mild steel (IF) with Young Modulus <math display="inline">E=2.06 GPa</math> and Poisson ratio <math display="inline">\nu=0.3</math>. Mises yield criterion was used for plasticity behaviour with non-linear isotropic hardening defined by <math display="inline">k(e^p) = 545(0.13+e^p)^{0.267} [MPa]</math>. A uniform friction of 0.15 was used for all the tools. A low (10kN) blank holder force was considered in this simulation.
2250
2251
Figure [[#img-25|25]] compares the punch force during the stamping stage obtained with both BST and EBST1 elements for the simulation and experimental values. Also for reference the average values of the simulations presented in the conference are included. Explicit and implicit simulations are considered as different curves. There is a remarkable coincidence between the experimental values and the results obtained with BST and EBST1 elements.
2252
2253
<div id='img-25'></div>
2254
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2255
|-
2256
|[[Image:draft_Samper_226033773-fig_n2.png|600px|Stamping of a S-rail. Punch force versus punch travel. Average of explicit and implicit results reported at the benchmark are also shown. ]]
2257
|- style="text-align: center; font-size: 75%;"
2258
| colspan="1" | '''Figure 25:''' Stamping of a S-rail. Punch force versus punch travel. Average of explicit and implicit results reported at the benchmark are also shown. 
2259
|}
2260
2261
Figure [[#img-26|26]] plots the <math display="inline">Z</math> coordinate along line B"&#8211;G" after springback stage. The top surface of the sheet does not remain plane due to some instabilities for the low blank holder force used. Results obtained with the simulations compare very well with the experimental results.
2262
2263
<div id='img-26'></div>
2264
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2265
|-
2266
|[[Image:draft_Samper_226033773-fig_n3.png|600px|Stamping of a S-rail. Z-coordinate along line B''&#8211;-G'' after springback. Average of explicit and implicit results reported at the benchmark are also shown. ]]
2267
|- style="text-align: center; font-size: 75%;"
2268
| colspan="1" | '''Figure 26:''' Stamping of a S-rail. Z-coordinate along line B''&#8211;-G'' after springback. Average of explicit and implicit results reported at the benchmark are also shown. 
2269
|}
2270
2271
==10 CONCLUDING REMARKS==
2272
2273
We have presented in the paper two alternative formulations for the rotation-free basic shell triangle (BST) using an assumed strain approach.  The simplest element of the family is based on an assumed constant curvature field expressed in terms of the nodal deflections of a patch of four elements and a constant membrane field computed from the standard linear interpolation of the displacements within each triangle. An enhanced version of the element is obtained by using a quadratic interpolation of the geometry in terms of the six patch nodes.  This allows to compute an assumed linear membrane strain field which improves the in-plane behaviour of the original element.  A simple and economic version of the element using a single integration point has been presented.  The efficiency of the different rotation-free shell triangles has been demonstrated in many examples of application including linear and non linear analysis of shells under static and dynamic loads, the inflation and de-inflation of membranes and a sheet stamping problem.
2274
2275
The enhanced rotation-free basic shell triangle element with a single integration point (the EBST1 element) has proven to be an excellent candidate for solving practical engineering shell and membrane problems involving complex geometry, dynamics, material non linearity and frictional contact conditions.
2276
2277
==ACKNOWLEDGEMENTS==
2278
2279
The problems analyzed with the explicit formulation were solved with the computer code STAMPACK <span id='citeF-42'></span>[[#cite-42|[42]]] where the rotation-free elements here presented have been implemented.  The support of the company QUANTECH (www.quantech.es) providing the code STAMPACK is gratefully acknowledged.
2280
2281
===BIBLIOGRAPHY===
2282
2283
<div id="cite-1"></div>
2284
'''[[#citeF-1|[1]]]''' E. Oñate. A review of some finite element families for thick and thin plate and shell analysis. Publication CIMNE N.53, May 1994.
2285
2286
<div id="cite-2"></div>
2287
'''[2]''' F.G. Flores, E. Oñate and F. Zárate. New assumed strain triangles for non-linear shell analysis. ''Computational Mechanics'', '''17''', 107&#8211;114, 1995.
2288
2289
<div id="cite-3"></div>
2290
'''[3]''' J.H. Argyris, M. Papadrakakis, C. Aportolopoulun and S. Koutsourelakis. The TRIC element. Theoretical and numerical investigation. ''Comput. Meth. Appl. Mech. Engrg.'', '''182''', 217&#8211;245, 2000.
2291
2292
<div id="cite-4"></div>
2293
'''[4]''' O.C. Zienkiewicz and R.L. Taylor. ''The finite element method. Solid Mechanics''. Vol II, Butterworth-Heinemann, 2000.
2294
2295
<div id="cite-5"></div>
2296
'''[5]'''  H. Stolarski, T. Belytschko and S.-H. Lee. A review of shell finite elements and corotational theories. ''Computational Mechanics Advances'', Vol. 2 N.2, North-Holland, 1995.
2297
2298
<div id="cite-6"></div>
2299
'''[[#citeF-6|[6]]]'''  E. Ramm and W.A. Wall. Shells in advanced computational environment. In ''V World Congress on Computational Mechanics'', J. Eberhardsteiner, H. Mang and F. Rammerstorfer (Eds.), Vienna, Austria, July 7&#8211;12, 2002. http://wccm.tuwien.ac.at.
2300
2301
<div id="cite-7"></div>
2302
'''[[#citeF-7|[7]]]'''  D. Bushnell and B.O. Almroth, “Finite difference energy method for non linear shell analysis”, ''J. Computers and Structures'', Vol. '''1''', 361, 1971.
2303
2304
<div id="cite-8"></div>
2305
'''[8]''' S.P. Timoshenko. ''Theory of Plates and Shells'', McGraw Hill, New York, 1971.
2306
2307
<div id="cite-9"></div>
2308
'''[[#citeF-9|[9]]]''' A.C. Ugural. ''Stresses in  Plates and Shells'', McGraw Hill, New York, 1981.
2309
2310
<div id="cite-10"></div>
2311
'''[10]''' R.A. Nay and S. Utku. An alternative to the finite element method. ''Variational Methods Eng.'', Vol. 1, 1972.
2312
2313
<div id="cite-11"></div>
2314
'''[11]''' J.K. Hampshire, B.H.V. Topping and H.C.  Chan. Three node triangular elements with one degree of freedom per node. ''Engng. Comput''. Vol. '''9''', pp. 49&#8211;62, 1992.
2315
2316
<div id="cite-12"></div>
2317
'''[12]''' R. Phaal and C.R. Calladine. A simple class of finite elements for plate and shell problems. I: Elements for beams and thin plates. ''Int. J. Num. Meth. Engng.'', Vol. '''35''', pp. 955&#8211;977, 1992.
2318
2319
<div id="cite-13"></div>
2320
'''[13]''' R. Phaal and C.R.  Calladine. A simple class of finite elements for plate and shell problems. II: An element for thin shells with only translational degrees of freedom. ''Int. J. Num. Meth. Engng.'', Vol. '''35''',  pp. 979&#8211;996, 1992.
2321
2322
<div id="cite-14"></div>
2323
'''[14]''' E. Oñate and Cervera M. Derivation of thin plate bending elements with one degree of freedom per node. ''Engineering Computations'', Vol. '''10''', pp 553&#8211;561, 1993.
2324
2325
<div id="cite-15"></div>
2326
'''[15]''' M. Brunet and F. Sabourin. Prediction of necking and wrinkles with a simplified shell element in sheet forming. ''Int. Conf. of Metal Forming Simulation in Industry'', Vol. II, pp. 27&#8211;48, B. Kröplin (Ed.), 1994.
2327
2328
<div id="cite-16"></div>
2329
'''[16]''' G. Rio, B. Tathi and H. Laurent. A new efficient finite element model of shell with only three degrees of freedom per node. Applications to industrial deep drawing test. in ''Recent Developments in Sheet Metal Forming Technology'', Ed. M.J.M. Barata Marques, 18th IDDRG Biennial Congress, Lisbon, 1994.
2330
2331
<div id="cite-17"></div>
2332
'''[17]''' G. Rio, B. Tathi and  H. Laurent. A new efficient finite element model of shell with only three degrees of freedom per node. Applications to industrial deep drawing test. in ''Recent Developments in Sheet Metal Forming Technology'', Ed. M.J.M. Barata Marques, 18th IDDRG Biennial Congress, Lisbon, 1994.
2333
2334
<div id="cite-18"></div>
2335
'''[18]'''  J. Rojek and E. Oñate. Sheet springback analysis using a simple shell triangle with translational degrees of freedom only. ''Int. J. of Forming Processes'', Vol. '''1''', No. 3, 275&#8211;296, 1998.
2336
2337
<div id="cite-19"></div>
2338
'''[19]'''  J. Rojek, E. Oñate and E. Postek. Application of explicit FE codes to simulation of sheet and bulk forming processes. ''J. of Materials Processing Technology'', Vols. '''80-81''', 620&#8211;627, 1998.
2339
2340
<div id="cite-20"></div>
2341
'''[20]'''  J. Jovicevic and E. Oñate. ''Analysis of beams and shells using a rotation-free finite element-finite volume formulation'', Monograph 43, CIMNE, Barcelona, 1999.
2342
2343
<div id="cite-21"></div>
2344
'''[[#citeF-21|[21]]]''' E. Oñate and F. Zárate. Rotation-free plate and shell triangles. ''Int. J. Num. Meth. Engng.'', '''47''', pp. 557&#8211;603, 2000.
2345
2346
<div id="cite-22"></div>
2347
'''[22]''' F. Cirak and M. Ortiz. Subdivision surfaces: A new paradigm for thin-shell finite element analysis. ''Int. J. Num. Meths in Engng'',  Vol. 47, 2000, 2039-2072.
2348
2349
<div id="cite-23"></div>
2350
'''[[#citeF-23|[23]]]''' F. Cirak and M. Ortiz. Fully <math display="inline">C^{1}</math>-conforming subdivision elements for finite deformations thin-shell analysis. ''Int. J. Num. Meths in Engng'', vol. 51, 2001, 813-833.
2351
2352
<div id="cite-24"></div>
2353
'''[[#citeF-24|[24]]]''' F.G. Flores and E. Oñate. A basic thin shell triangle with only translational DOFs for large strain plasticity. ''Int. J. Num. Meths in Engng'', Vol. '''51''', pp 57-83, 2001.
2354
2355
<div id="cite-25"></div>
2356
'''[25]''' G. Engel, K. Garikipati, T.J.R. Hughes, M.G. Larson, L. Mazzei and R.L. Taylor. Continuous/discontinuous finite element approximation of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity. ''Comput. Methods Appl. Mech. Engrg.'', Vol. 191, 3669&#8211;3750, 2002.
2357
2358
<div id="cite-26"></div>
2359
'''[26]'''  E. Oñate, P. Cendoya and J. Miquel. Non linear explicit dynamic analysis of shells using the BST rotation-free triangle. ''Engineering Computations'', '''19''' (6), 662&#8211;706, 2002.
2360
2361
<div id="cite-27"></div>
2362
'''[[#citeF-27|[27]]]''' F.G. Flores and E. Oñate. Improvements in the membrane behaviour of the three node rotation-free BST shell triangle  using an assumed strain approach. ''Computer Methods in Applied Mechanics and Engineering'', in press, 2003.
2363
2364
<div id="cite-28"></div>
2365
'''[28]'''  O.C. Zienkiewicz and E. Oñate. Finite Elements vs. Finite Volumes. Is there really a choice?. ''Nonlinear Computational Mechanics''. State of the Art. (Eds. P. Wriggers and R. Wagner). Springer Verlag, Heidelberg, 1991.
2366
2367
<div id="cite-29"></div>
2368
'''[[#citeF-29|[29]]]''' E. Oñate, M.  Cervera  and O.C. Zienkiewicz. A finite volume format for structural mechanics. ''Int. J. Num. Meth. Engng.'', '''37''', pp. 181&#8211;201, 1994.
2369
2370
<div id="cite-30"></div>
2371
'''[30]''' R. Hill. A Theory of the Yielding and Plastic Flow of Anisotropic Metals. ''Proc. Royal Society London'', Vol. '''A193''', pp. 281, 1948.
2372
2373
<div id="cite-31"></div>
2374
'''[31]''' R.W. Ogden. Large deformation isotropic elasticity: on the correlation of theory and experiments for incompressible rubberlike solids. ''Proceedings of the Royal Society of London'', Vol. '''A326''', pp. 565&#8211;584, 1972.
2375
2376
<div id="cite-32"></div>
2377
'''[32]''' H.C. Huang, ''Static and Dynamic Analysis of Plates and Shells'', page 40, Springer-Verlag, Berlin, 1989.
2378
2379
<div id="cite-33"></div>
2380
'''[33]''' E.N. Dvorkin and K.J. Bathe. A continuum mechanics based four node shell element for general non-linear analysis. ''Eng. Comp.'', '''1''', 77&#8211;88, 1984.
2381
2382
<div id="cite-34"></div>
2383
'''[34]''' A. Needleman. Inflation of spherical rubber ballons. ''Int. J. of Solids and Structures'', '''13''', 409&#8211;421, 1977.
2384
2385
<div id="cite-35"></div>
2386
'''[[#citeF-35|[35]]]'''  Hibbit, Karlson and Sorensen Inc. ABAQUS, version 5.8, Pawtucket, USA, 1998.
2387
2388
<div id="cite-36"></div>
2389
'''[36]''' WHAMS-3D. An explicit 3D finite element program. KBS2  Inc., Willow Springs, Illinois 60480, USA.
2390
2391
<div id="cite-37"></div>
2392
'''[[#citeF-37|[37]]]'''  H.A. Balmer and E.A. Witmer. Theoretical experimental correlation of large dynamic and permanent deformation of impulsively loaded simple structures. ''Air force flight Dynamic Lab. Rep. FDQ-TDR-64-108'', Wright-Patterson AFB, Ohio, USA, 1964.
2393
2394
<div id="cite-38"></div>
2395
'''[38]'''  H. Stolarski, T. Belytschko and N. Carpenter. A simple triangular curved shell element. ''Eng. Comput.'', Vol. 1, 210&#8211;218, 1984.
2396
2397
<div id="cite-39"></div>
2398
'''[39]'''  P.O. Marklund and L. Nilsson. Simulation of airbag inflation processes using a coupled fluid structure approach. ''Computational Mechanics'', '''29''', 289&#8211;297, 2002.
2399
2400
<div id="cite-40"></div>
2401
'''[[#citeF-40|[40]]]''' NUMISHEET'96, ''Third International Conference and Workshop on Numerical Simulation of 3D Sheet Forming Processes, NUMISHEET'96'', E.H. Lee, G.L. Kinzel and R.H. Wagoner (Eds.), Dearbon-Michigan, USA, 1996.
2402
2403
<div id="cite-41"></div>
2404
'''[[#citeF-41|[41]]]'''  <code>http://rclsgi.eng.ohio-state.edu/%Elee-j-k/numisheet96/</code>
2405
2406
<div id="cite-42"></div>
2407
'''[[#citeF-42|[42]]]''' STAMPACK. ''A General Finite Element System for Sheet Stamping and Forming Problems'', Quantech ATZ, Barcelona, Spain, 2003 (www.quantech.es).
2408
2409
==APPENDIX==
2410
2411
==11 Curvature matrix for the BST element==
2412
2413
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2414
|-
2415
| 
2416
{| style="text-align: left; margin:auto;" 
2417
|-
2418
| style="text-align: center;" | <math>\delta{\boldsymbol \kappa }=\mathbf{B}_{b} \times \mathbf{t}_3 \delta \mathbf{a}^{p}  </math>
2419
|}
2420
|}
2421
2422
with
2423
2424
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2425
|-
2426
| 
2427
{| style="text-align: left; margin:auto;" 
2428
|-
2429
| style="text-align: center;" | <math>\begin{array}{c} \\ \delta \mathbf{a}^{p}\\ 18\times{1} \end{array} =[\delta \mathbf{u}_{1}^{T},\delta \mathbf{u}_{2}^{T},\delta \mathbf{u}_{3} ^{T},\delta \mathbf{u}_{4}^{T},\delta \mathbf{u}_{5}^{T},\delta \mathbf{u}_{6}^{T}]^{T}  </math>
2430
|}
2431
|}
2432
2433
and 
2434
2435
<math>\mathbf{B}_{b}^{T}=</math>
2436
2437
2438
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2439
|- style="border-top: 2px solid;"
2440
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">+L^{M}_{2,1} L^{2}_{2,1}    +L^{M}_{3,1} L^{3}_{3,1} </math> 
2441
| style="border-left: 2px solid;border-right: 2px solid;" | <math>+L^{M}_{2,2} L^{2}_{2,2}    +L^{M}_{3,2} L^{3}_{3,2} </math>
2442
| style="border-left: 2px solid;border-right: 2px solid;" | <math>+L^{M}_{2,2} L^{2}_{2,1} +L^{M}_{2,1} L^{2}_{2,2}    +L^{M}_{3,2} L^{3}_{3,1} +L^{M}_{3,1} L^{3}_{3,2} </math>
2443
|- style="border-top: 2px solid;"
2444
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> L^{M}_{1,1} L^{1}_{3,1}    +L^{M}_{3,1} L^{3}_{2,1} </math> 
2445
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{1,2} L^{1}_{3,2}    +L^{M}_{3,2} L^{3}_{2,2} </math>
2446
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{1,2} L^{1}_{3,1} +L^{M}_{1,1} L^{1}_{3,2}    +L^{M}_{3,2} L^{3}_{2,1} +L^{M}_{3,1} L^{3}_{2,2} </math>
2447
|- style="border-top: 2px solid;"
2448
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> L^{M}_{1,1} L^{1}_{2,1}    +L^{M}_{2,1} L^{2}_{3,1} </math> 
2449
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{1,2} L^{1}_{2,2}    +L^{M}_{2,2} L^{2}_{3,2} </math>
2450
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{1,2} L^{1}_{2,1} +L^{M}_{1,1} L^{1}_{j,3}    +L^{M}_{2,2} L^{2}_{3,1} +L^{M}_{2,1} L^{2}_{3,2} </math>
2451
|- style="border-top: 2px solid;"
2452
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L^{M}_{1,1} L^{1}_{1,1} </math> 
2453
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L^{M}_{1,2} L^{1}_{1,2} </math>
2454
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L^{M}_{1,2} L^{1}_{1,1} +L^{M}_{1,1} L^{1}_{1,3} </math>
2455
|- style="border-top: 2px solid;"
2456
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L^{M}_{2,1} L^{2}_{1,1} </math> 
2457
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L^{M}_{2,2} L^{2}_{1,2} </math>
2458
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L^{M}_{2,2} L^{2}_{1,1} +L^{M}_{2,1} L^{2}_{1,3} </math>
2459
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2460
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L^{M}_{3,1} L^{3}_{1,1} </math> 
2461
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L^{M}_{3,2} L^{3}_{1,2} </math>
2462
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L^{M}_{3,2} L^{3}_{1,1} +L^{M}_{3,1} L^{3}_{1,3} </math>
2463
2464
|}
2465
2466
<math display="inline">-2</math> 
2467
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2468
|- style="border-top: 2px solid;"
2469
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> L^{M}_{1,1}\rho _{11}^{1}+L^{M}_{1,2}\rho _{11}^{2} </math> 
2470
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{1,1}\rho _{22}^{1}+L^{M}_{i,2}\rho _{22}^{2} </math>
2471
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{1,1}\rho _{12}^{1}+L^{M}_{1,2}\rho _{12}^{2} </math>
2472
|- style="border-top: 2px solid;"
2473
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> L^{M}_{2,1}\rho _{11}^{1}+L^{M}_{2,2}\rho _{11}^{2} </math> 
2474
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{2,1}\rho _{22}^{1}+L^{M}_{2,2}\rho _{22}^{2} </math>
2475
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{2,1}\rho _{12}^{1}+L^{M}_{2,2}\rho _{12}^{2} </math>
2476
|- style="border-top: 2px solid;"
2477
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> L^{M}_{3,1}\rho _{11}^{1}+L^{M}_{3,2}\rho _{11}^{2} </math> 
2478
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{3,1}\rho _{22}^{1}+L^{M}_{3,2}\rho _{22}^{2} </math>
2479
| style="border-left: 2px solid;border-right: 2px solid;" | <math> L^{M}_{3,1}\rho _{12}^{1}+L^{M}_{3,2}\rho _{12}^{2} </math>
2480
|- style="border-top: 2px solid;"
2481
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">0</math> 
2482
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2483
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2484
|- style="border-top: 2px solid;"
2485
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">0</math> 
2486
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2487
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2488
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2489
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">0</math> 
2490
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2491
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2492
2493
|}
2494
2495
<br/><br/>
2496
2497
==12 Membrane strain matrix and curvature matrix for the EBST element==
2498
2499
===12.1 Membrane strain matrix===
2500
2501
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2502
|-
2503
| 
2504
{| style="text-align: left; margin:auto;" 
2505
|-
2506
| style="text-align: center;" | <math>\delta {\boldsymbol \varepsilon }_m ={B}_m \delta {a}^p  </math>
2507
|}
2508
|}
2509
2510
<math>\mathbf{B}_{m}^{T}=\frac{1}{3}</math>
2511
2512
2513
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2514
|- style="border-top: 2px solid;"
2515
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{1,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}   + N^{2}_{1,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}   + N^{3}_{1,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1} </math> 
2516
| style="border-left: 2px solid;border-right: 2px solid;" | <math> N^{1}_{1,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2}   + N^{2}_{1,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2}   + N^{3}_{1,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math>
2517
|- style="border-top: 2px solid;"
2518
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{2,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}   + N^{2}_{2,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}   + N^{3}_{2,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1} </math> 
2519
| style="border-left: 2px solid;border-right: 2px solid;" | <math> N^{1}_{2,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2}   + N^{2}_{2,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2}   + N^{3}_{2,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math>
2520
|- style="border-top: 2px solid;"
2521
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{3,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}   + N^{2}_{3,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}   + N^{3}_{3,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1} </math> 
2522
| style="border-left: 2px solid;border-right: 2px solid;" | <math> N^{1}_{3,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2}   + N^{2}_{3,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2}   + N^{3}_{3,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math>
2523
|- style="border-top: 2px solid;"
2524
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{4,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1} </math> 
2525
| style="border-left: 2px solid;border-right: 2px solid;" | <math> N^{1}_{4,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2} </math>
2526
|- style="border-top: 2px solid;"
2527
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{2}_{5,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1} </math> 
2528
| style="border-left: 2px solid;border-right: 2px solid;" | <math> N^{2}_{5,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2} </math>
2529
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2530
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{3}_{6,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1} </math> 
2531
| style="border-left: 2px solid;border-right: 2px solid;" | <math> N^{3}_{6,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math>
2532
2533
|}
2534
2535
2536
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2537
|- style="border-top: 2px solid;"
2538
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{1,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}    +N^{1}_{1,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2}    +N^{2}_{1,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}    +N^{2}_{1,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2}    +N^{3}_{1,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1}    +N^{3}_{1,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math> 
2539
|- style="border-top: 2px solid;"
2540
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{2,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}    +N^{1}_{2,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2}   + N^{2}_{2,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}    +N^{2}_{2,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2}   + N^{3}_{2,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1}    +N^{3}_{2,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math> 
2541
|- style="border-top: 2px solid;"
2542
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{3,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}    +N^{1}_{3,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2}    +N^{2}_{3,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}    +N^{2}_{3,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2}    +N^{3}_{3,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1}    +N^{3}_{3,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math> 
2543
|- style="border-top: 2px solid;"
2544
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{1}_{4,2}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }1}    +N^{1}_{4,1}\mathbf{\boldsymbol \varphi }^{1}_{^{\prime }2} </math> 
2545
|- style="border-top: 2px solid;"
2546
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{2}_{5,2}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }1}    +N^{2}_{5,1}\mathbf{\boldsymbol \varphi }^{2}_{^{\prime }2} </math> 
2547
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2548
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline"> N^{3}_{6,2}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }1}    +N^{3}_{6,1}\mathbf{\boldsymbol \varphi }^{3}_{^{\prime }2} </math> 
2549
2550
|}
2551
2552
===12.2 Curvature matrix===
2553
2554
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2555
|-
2556
| 
2557
{| style="text-align: left; margin:auto;" 
2558
|-
2559
| style="text-align: center;" | <math>\delta {\boldsymbol \kappa } ={B}_b \times \mathbf{t}_3 \delta {a}^p  </math>
2560
|}
2561
|}
2562
2563
<math>\mathbf{B}_{b}^{T}=2</math>
2564
2565
2566
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2567
|- style="border-top: 2px solid;"
2568
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,1}\left(N_{1,1}\right)_{G_{1}}   +L_{2,1}\left(N_{1,1}\right)_{G_{2}}   +L_{3,1}\left(N_{1,1}\right)_{G_{3}}</math> 
2569
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L_{1,2}\left(N_{1,2}\right)_{G_{1}}   +L_{2,2}\left(N_{1,2}\right)_{G_{2}}   +L_{3,2}\left(N_{1,2}\right)_{G_{3}}</math>
2570
|- style="border-top: 2px solid;"
2571
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,1}\left(N_{2,1}\right)_{G_{1}}   +L_{2,1}\left(N_{2,1}\right)_{G_{2}}   +L_{3,1}\left(N_{2,1}\right)_{G_{3}}</math> 
2572
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L_{1,2}\left(N_{2,2}\right)_{G_{1}}   +L_{2,2}\left(N_{2,2}\right)_{G_{2}}   +L_{3,2}\left(N_{2,2}\right)_{G_{3}}</math>
2573
|- style="border-top: 2px solid;"
2574
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,1}\left(N_{3,1}\right)_{G_{1}}   +L_{2,1}\left(N_{3,1}\right)_{G_{2}}   +L_{3,1}\left(N_{3,1}\right)_{G_{3}}</math> 
2575
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L_{1,2}\left(N_{3,2}\right)_{G_{1}}   +L_{2,2}\left(N_{3,2}\right)_{G_{2}}   +L_{3,2}\left(N_{3,2}\right)_{G_{3}}</math>
2576
|- style="border-top: 2px solid;"
2577
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,1}\left(N_{4,1}\right)_{G_{1}}</math> 
2578
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L_{1,2}\left(N_{4,2}\right)_{G_{1}}</math>
2579
|- style="border-top: 2px solid;"
2580
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{2,1}\left(N_{5,1}\right)_{G_{2}}</math> 
2581
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L_{2,2}\left(N_{5,2}\right)_{G_{2}}</math>
2582
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2583
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{3,1}\left(N_{6,1}\right)_{G_{3}}</math> 
2584
| style="border-left: 2px solid;border-right: 2px solid;" | <math>L_{3,2}\left(N_{6,2}\right)_{G_{3}}</math>
2585
2586
|}
2587
2588
2589
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2590
|- style="border-top: 2px solid;"
2591
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,2}\left(N_{1,1}\right)_{G_{1}}+L_{1,1}\left(N_{1,2}\right)_{G_{1}}   +L_{2,2}\left(N_{1,1}\right)_{G_{2}}+L_{2,1}\left(N_{1,2}\right)_{G_{2}}   +L_{3,2}\left(N_{1,1}\right)_{G_{3}}+L_{3,1}\left(N_{1,2}\right)_{G_{3}}</math> 
2592
|- style="border-top: 2px solid;"
2593
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,2}\left(N_{2,1}\right)_{G_{1}}+L_{1,1}\left(N_{2,2}\right)_{G_{1}}   +L_{2,2}\left(N_{2,1}\right)_{G_{2}}+L_{2,1}\left(N_{2,2}\right)_{G_{2}}   +L_{3,2}\left(N_{2,1}\right)_{G_{3}}+L_{3,1}\left(N_{2,2}\right)_{G_{3}}</math> 
2594
|- style="border-top: 2px solid;"
2595
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,2}\left(N_{3,1}\right)_{G_{1}}+L_{1,1}\left(N_{j,3}\right)_{G_{1}}   +L_{2,2}\left(N_{3,1}\right)_{G_{2}}+L_{2,1}\left(N_{j,3}\right)_{G_{2}}   +L_{3,2}\left(N_{3,1}\right)_{G_{3}}+L_{3,1}\left(N_{j,3}\right)_{G_{3}}</math> 
2596
|- style="border-top: 2px solid;"
2597
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{1,2}\left(N_{4,1}\right)_{G_{1}}+L_{1,1}\left(N_{4,3}\right)_{G_{1}}</math> 
2598
|- style="border-top: 2px solid;"
2599
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{2,2}\left(N_{5,1}\right)_{G_{2}}+L_{2,1}\left(N_{5,3}\right)_{G_{2}}</math> 
2600
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2601
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">L_{3,2}\left(N_{6,1}\right)_{G_{3}}+L_{3,1}\left(N_{6,3}\right)_{G_{6}}</math> 
2602
2603
|}
2604
2605
<math>-2</math>
2606
2607
2608
{| class="wikitable" style="text-align: center; margin: 1em auto;"
2609
|- style="border-top: 2px solid;"
2610
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\left(L_{1,1}\rho _{11}^{1}+L_{1,2}\rho _{11}^{2}\right)</math> 
2611
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\left(L_{1,1}\rho _{22}^{1}+L_{i,2}\rho _{22}^{2}\right)</math>
2612
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\left(L_{1,1}\rho _{12}^{1}+L_{1,2}\rho _{12}^{2}\right)</math>
2613
|- style="border-top: 2px solid;"
2614
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\left(L_{2,1}\rho _{11}^{1}+L_{2,2}\rho _{11}^{2}\right)</math> 
2615
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\left(L_{2,1}\rho _{22}^{1}+L_{2,2}\rho _{22}^{2}\right)</math>
2616
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\left(L_{2,1}\rho _{12}^{1}+L_{2,2}\rho _{12}^{2}\right)</math>
2617
|- style="border-top: 2px solid;"
2618
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\left(L_{3,1}\rho _{11}^{1}+L_{3,2}\rho _{11}^{2}\right)</math> 
2619
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\left(L_{3,1}\rho _{22}^{1}+L_{3,2}\rho _{22}^{2}\right)</math>
2620
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\left(L_{3,1}\rho _{12}^{1}+L_{3,2}\rho _{12}^{2}\right)</math>
2621
|- style="border-top: 2px solid;"
2622
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">0</math> 
2623
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2624
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2625
|- style="border-top: 2px solid;"
2626
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">0</math> 
2627
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2628
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2629
|- style="border-top: 2px solid;border-bottom: 2px solid;"
2630
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">0</math> 
2631
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2632
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
2633
2634
|}
2635
2636
In this last expression <math display="inline">L_{i,j} =L_{i,j}^{M}</math>
2637

Return to Onate Flores 2003a.

Back to Top

Document information

Published on 01/01/2005

DOI: 10.1016/j.cma.2004.07.039
Licence: CC BY-NC-SA license

Document Score

0

Times cited: 43
Views 80
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?