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
2
3
==Abstract==
4
5
In this paper an assumed strain approach is presented in order to improve the  membrane behaviour of a thin shell triangular element. The so called Basic Shell Triangle (BST) has three  nodes with only translational degrees of freedom and is based on a Total  Lagrangian Formulation. As in the original BST element the curvatures are computed  resorting to the surrounding elements (patch of four elements). Membrane  strains are now also computed from the same patch of elements which leads to a  non-conforming  membrane behaviour. Despite this non-conformity the element  passes the patch test. Large strain plasticity is considered using a logarithmic  strain-stress pair. A plane stress behaviour with an additive decomposition of  elastic and plastic strains is assumed. A hyperplastic law is considered for  the elastic part while for the plastic part an anisotropic quadratic (Hill)  yield function with non-linear isotropic hardening is adopted. The element,  termed CBST, has been implemented  in an explicit (hydro-)code adequate to simulate sheet-stamping processes and  in an implicit static/dynamic code. Several examples are given showing the good performance of the enhnaced rotation-free shell triangle.
6
7
8
9
==1 Introduction==
10
11
Development of efficient and robust shell finite elements for modelling large-scale  industrial applications has been a field of constant research in the last  thirty years. The continuining efforts of many scientists has focused in the direction of deriving simple elements (preferably triangles) applicable to large scale computation, typical of practical shell problems. Some recent successful shell triangles are reported in <span id='citeF-1'></span>[[#cite-1|[1,2]]]. For a comprehensive review of shell elements see <span id='citeF-3'></span>[[#cite-3|[3]]]. From the theoretical point of view Kirchhoff hypothesis are  enough for the simulation of the essential aspects of most shells problems,  including structures in civil, mechanical and naval engineering, sheet metal  forming, crash-worthiness analysis and others. Due the well-known difficulties  associated to <math display="inline">C^{1}</math> continuity, “thick shell” approaches that need <math display="inline">C^{0}</math>  continuity only have been extensively used, although they are computationally  more expensive due to the number of degrees of freedom necessary to obtain thin   solutions with similar precision. The development of shell elements without  transversal shear strains then collides with the need for C<math display="inline">^{1}</math> continuity, which is not  simple to satisfy in general conditions and has led to the development of  non-conforming approaches. Among the many options of this kind, the use of shell elements with only  displacements as degrees of freedom (rotation-free) is very attractive from  the computational point of view.
12
13
Few shell elements exists in the literature with only translations as degrees of  freedom. They are generally based on non-conforming approaches. A state of the  art can be found in Ref. <span id='citeF-4'></span>[[#cite-4|[4]]]. Recently Cirak y  Ortiz <span id='citeF-5'></span>[[#cite-5|[5]]] have developed a conforming thin shell element,  but it departs in some aspects from the standard finite element formulation. The  use of rotation-free shell elements in large scale simulations of sheet metal forming processes  is growing <span id='citeF-6'></span>[[#cite-6|[6]]]-<span id='citeF-9'></span>[[#cite-9|[9]]], specially in  explicit integration codes.
14
15
In Ref. <span id='citeF-10'></span>[[#cite-10|[10]]] we presented a thin shell triangular element with only  displacements as degrees of freedom, based on a total lagrangian formulation.  The element is an extension of the rotation-free Basic Shell Triangle (BST) originally developed by Oñate y  Zárate <span id='citeF-4'></span>[[#cite-4|[4]]] using an updated lagrangian formulation  The element has three nodes and a linear approximation of the displacement field. The membrane part is identical to the standard  constant strain triangle. The bending part is also constant and is computed  from an integration over the element boundary using information from the  deflection gradients of the adjacent elements. This leads to a non-conforming element for  the bending part, but it converges to the right solution and is robust. The BST element  has  a very good performance in bending for structured meshes, but it is less accurate  for non-structured grids. On the other hand, for membrane dominated problems,  as is the case for most sheet metal forming simulations, it requires  fine  discretizations associated to the constant strain triangle behaviour.
16
17
In this paper an extension of the rotation-free BST element is presented.  The displacement gradient terms needed for the membrane and bending strains are  computed from the nodal displacements of a patch of element using a new procedure. This allows to improve the membrane  behaviour and obtain a smoother curvature field. The outline of this  paper is as follows. In sections 2 and 3 we introduce the kinematics of the  shell and the constitutive models used in the numerical simulations. Sections  4 to 6 describe the finite element approximation, including the new  proposal for the gradient evaluation from a patch of elements, the computation of  the metric tensor and the curvatures and the derivation of the necessary  element expressions for the computer implementation. Sections 7 and 8 sumarize  the numerical experiments performed in the linear and non-linear range. Numerical results   demonstrate the excellent performance of the enhanced three node rotation-free shell triangle.
18
19
==2 Shell kinematics==
20
21
A  summary of the most relevant hypothesis related to the kinematic  behaviour of the shell are presented. Further details may be found in the wide  literature dedicated to this field <span id='citeF-3'></span>[[#cite-3|[3]]].
22
23
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 any point in the undeformed and the deformed configurations  can be respectively written as a function of the position of the middle  surface <math display="inline">\mathbf{\varphi }</math> and the normal <math display="inline">\mathbf{t}_{3}</math> at the point as
24
25
{| class="formulaSCP" style="width: 100%; text-align: left;" 
26
|-
27
| 
28
{| style="text-align: left; margin:auto;width: 100%;" 
29
|-
30
| style="text-align: center;" | <math>\mathbf{x}^{0}\left( \xi _{1},\xi _{2},\zeta \right)    =\mathbf{\varphi }^{0}\left( \xi _{1},\xi _{2}\right) +\lambda \mathbf{t}_{3}^{0}</math>
31
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
32
|-
33
| style="text-align: center;" | <math> \mathbf{x}\left( \xi _{1},\xi _{2},\zeta \right)    =\mathbf{\varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3} </math>
34
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
35
|}
36
|}
37
38
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.
39
40
{| class="formulaSCP" style="width: 100%; text-align: left;" 
41
|-
42
| 
43
{| style="text-align: left; margin:auto;width: 100%;" 
44
|-
45
| style="text-align: center;" | <math>  \lambda =\frac{h}{h^{0}} </math>
46
|}
47
|}
48
49
A convective system is defined at each point as
50
51
{| class="formulaSCP" style="width: 100%; text-align: left;" 
52
|-
53
| 
54
{| style="text-align: left; margin:auto;width: 100%;" 
55
|-
56
| style="text-align: center;" | <math>\mathbf{g}_{i}\left( \mathbf{\xi }\right) =\frac{\partial \mathbf{x}}{\partial \xi _{i}} </math>
57
|}
58
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
59
|}
60
61
{| class="formulaSCP" style="width: 100%; text-align: left;" 
62
|-
63
| 
64
{| style="text-align: left; margin:auto;width: 100%;" 
65
|-
66
| style="text-align: center;" | <math>\mathbf{g}_{\alpha }\left( \mathbf{\xi }\right)    =\frac{\partial \left( \mathbf{\varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \xi _{\alpha }}=\mathbf{\varphi }_{^{\prime }\alpha } +\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }\alpha }</math>
67
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
68
|-
69
| style="text-align: center;" | <math> \mathbf{g}_{3}\left( \mathbf{\xi }\right)    =\frac{\partial \left( \mathbf{\varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \zeta }=\lambda \mathbf{t}_{3} </math>
70
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
71
|}
72
|}
73
74
This can be particularized for the points on the middle surface as
75
76
{| class="formulaSCP" style="width: 100%; text-align: left;" 
77
|-
78
| 
79
{| style="text-align: left; margin:auto;width: 100%;" 
80
|-
81
| style="text-align: center;" | <math>\mathbf{a}_{\alpha }    =\mathbf{g}_{\alpha }\left( \zeta=0\right) =\mathbf{\varphi }_{^{\prime }\alpha }</math>
82
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
83
|-
84
| style="text-align: center;" | <math> \mathbf{a}_{3}    =\mathbf{g}_{3}\left( \zeta=0\right) =\lambda  \mathbf{t}_{3} </math>
85
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
86
|}
87
|}
88
89
This allows to introduce the covariant (first fundamental form) and  contravariant metric tensors of the middle surface as
90
91
{| class="formulaSCP" style="width: 100%; text-align: left;" 
92
|-
93
| 
94
{| style="text-align: left; margin:auto;width: 100%;" 
95
|-
96
| style="text-align: center;" | <math>a_{\alpha \beta }=\mathbf{a}_{\alpha }\cdot \mathbf{a}_{\beta } </math>
97
|}
98
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
99
|}
100
101
{| class="formulaSCP" style="width: 100%; text-align: left;" 
102
|-
103
| 
104
{| style="text-align: left; margin:auto;width: 100%;" 
105
|-
106
| style="text-align: center;" | <math>a^{\alpha \beta }=\mathbf{a}^{\alpha }\cdot \mathbf{a}^{\beta }=\mathbf{\tilde  {\varphi }}_{^{\prime }\alpha }\cdot \mathbf{\tilde{\varphi }}_{^{\prime }\beta } </math>
107
|}
108
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
109
|}
110
111
and the curvatures (second fundamental form) of the middle surface as
112
113
{| class="formulaSCP" style="width: 100%; text-align: left;" 
114
|-
115
| 
116
{| style="text-align: left; margin:auto;width: 100%;" 
117
|-
118
| style="text-align: center;" | <math>\kappa _{\alpha \beta }=\frac{1}{2}\left( \mathbf{\varphi }_{^{\prime }\alpha  }\cdot \mathbf{t}_{3^{\prime }\beta }+\mathbf{\varphi }_{^{\prime }\beta } \cdot \mathbf{t}_{3^{\prime }\alpha }\right) </math>
119
|}
120
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
121
|}
122
123
The deformation gradient can be written as
124
125
{| class="formulaSCP" style="width: 100%; text-align: left;" 
126
|-
127
| 
128
{| style="text-align: left; margin:auto;width: 100%;" 
129
|-
130
| style="text-align: center;" | <math>\mathbf{F=}\left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }1}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }1} & \mathbf{\varphi }_{^{\prime }2}+\zeta \left( \lambda  \mathbf{t}_{3}\right) _{^{\prime }2} & \lambda \mathbf{t}_{3} \end{array}  \right] </math>
131
|}
132
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
133
|}
134
135
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) is
136
137
{| class="formulaSCP" style="width: 100%; text-align: left;" 
138
|-
139
| 
140
{| style="text-align: left; margin:auto;width: 100%;" 
141
|-
142
| style="text-align: center;" | <math>\mathbf{U}^{2}    =\left[ \begin{array}{c} \mathbf{\varphi }_{^{\prime }1}^{T}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }1}^{T}\\
143
 \mathbf{\varphi }_{^{\prime }2}^{T}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }2}^{T}\\
144
 \lambda \mathbf{t}_{3}^{T} \end{array}  \right] \left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }1}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }1} & \mathbf{\varphi }_{^{\prime }2}+\zeta \left( \lambda  \mathbf{t}_{3}\right) _{^{\prime }2} & \lambda \mathbf{t}_{3} \end{array}  \right]</math>
145
|-
146
| style="text-align: center;" | <math>   =\left[ \begin{array}{ccc} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} & 0\\
147
 \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} &  \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2} & 0\\
148
 0 & 0 & \lambda ^{2} \end{array}  \right]</math>
149
|-
150
| style="text-align: center;" | <math>   +\zeta \lambda \left[ \begin{array}{ccc} 2\mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{t}_{3^{\prime }1} & \mathbf{\varphi  }_{^{\prime }1}\cdot \mathbf{t}_{3^{\prime }2}+\mathbf{\varphi }_{^{\prime }2} \cdot \mathbf{t}_{3^{\prime }1} & 0 \\
151
 \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{t}_{3^{\prime }2}+\mathbf{\varphi  }_{^{\prime }2}\cdot \mathbf{t}_{3^{\prime }1} & \mathbf{\varphi }_{^{\prime } 2}\cdot \mathbf{t}_{3^{\prime }2} & 0\\
152
 0 & 0 & 0  \end{array}  \right]</math>
153
|-
154
| style="text-align: center;" | <math>   +\zeta ^{2}\lambda ^{2}\left[ \begin{array}{ccc} \mathbf{t}_{3^{\prime }1}\cdot \mathbf{t}_{3^{\prime }1} & \mathbf{t}_{3^{\prime  }1}\cdot \mathbf{t}_{3^{\prime }2} & 0\\
155
 \mathbf{t}_{3^{\prime }1}\cdot \mathbf{t}_{3^{\prime }2} & \mathbf{t}_{3^{\prime  }2}\cdot \mathbf{t}_{3^{\prime }2} & 0\\
156
 0 & 0 & 0  \end{array}  \right] </math>
157
|}
158
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
159
|}
160
161
where the derivatives of the thickness ratio <math display="inline">\lambda _{^{\prime }a}</math> have been  neglected. Neglecting also the term associated to <math display="inline">\zeta ^{2}</math> and introducing  the definition of the covariant metric tensor <math display="inline">a_{\alpha \beta }</math> at the middle  surface and the curvatures <math display="inline">k_{\alpha \beta }</math> gives
162
163
{| class="formulaSCP" style="width: 100%; text-align: left;" 
164
|-
165
| 
166
{| style="text-align: left; margin:auto;width: 100%;" 
167
|-
168
| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{ccc} a_{11}+2k_{11}\zeta \lambda & a_{12}+2k_{12}\zeta \lambda & 0\\
169
 a_{12}+2k_{12}\zeta \lambda & a_{22}+2k_{22}\zeta \lambda & 0\\
170
 0 & 0 & \lambda ^{2} \end{array}  \right] </math>
171
|}
172
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
173
|}
174
175
Above expression shows that <math display="inline">\mathbf{U}^{2}</math>   is non-zero 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
176
177
{| class="formulaSCP" style="width: 100%; text-align: left;" 
178
|-
179
| 
180
{| style="text-align: left; margin:auto;width: 100%;" 
181
|-
182
| style="text-align: center;" | <math>\chi _{ij}=\kappa _{ij}-\kappa _{ij}^{0} </math>
183
|}
184
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
185
|}
186
187
For computational convenience the following approximate expression (which is  exact for initially flat surfaces) will be adopted
188
189
{| class="formulaSCP" style="width: 100%; text-align: left;" 
190
|-
191
| 
192
{| style="text-align: left; margin:auto;width: 100%;" 
193
|-
194
| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{ccc} a_{11}+2\chi _{11}\zeta \lambda & a_{12}+2\chi _{12}\zeta \lambda & 0\\
195
 a_{12}+2\chi _{12}\zeta \lambda & a_{22}+2\chi _{22}\zeta \lambda & 0\\
196
 0 & 0 & \lambda ^{2} \end{array}  \right] </math>
197
|}
198
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
199
|}
200
201
Above expression of <math display="inline">\mathbf{U}^{2}</math> 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
202
203
{| class="formulaSCP" style="width: 100%; text-align: left;" 
204
|-
205
| 
206
{| style="text-align: left; margin:auto;width: 100%;" 
207
|-
208
| style="text-align: center;" | <math>\mathbf{U=}\sum _{\alpha=1}^{3}\lambda _{\alpha } \mathbf{r}_{\alpha } \otimes \mathbf{r}_{\alpha } </math>
209
|}
210
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
211
|}
212
213
where <math display="inline">\lambda _{\alpha }</math> and <math display="inline">\mathbf{r}_{\alpha }</math> are the eigenvalues and  eigenvectors of the right stretch tensor <math display="inline">\mathbf{U}</math>.
214
215
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
216
217
{| class="formulaSCP" style="width: 100%; text-align: left;" 
218
|-
219
| 
220
{| style="text-align: left; margin:auto;width: 100%;" 
221
|-
222
| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=}\left[ \begin{array}{ccc} \varepsilon _{11} & \varepsilon _{21} & 0\\
223
 \varepsilon _{12} & \varepsilon _{22} & 0\\
224
 0 & 0 & \varepsilon _{33} \end{array}  \right] =\sum _{\alpha=1}^{3}\ln \left( \lambda _{\alpha }\right) \mathbf{r}_{\alpha }\otimes \mathbf{r}_{\alpha } </math>
225
|}
226
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
227
|}
228
229
Consistently, the Hencky stress tensor <math display="inline">\mathbf{T}</math>  will be  used as the stress measure. Using a total lagrangian formulation it is convenient to work with  the the second Piola-Kirchhoff stresses (<math display="inline">\mathbf{S}</math>) for the evaluation of  the residual forces. We define the rotated tensors as
230
231
{| class="formulaSCP" style="width: 100%; text-align: left;" 
232
|-
233
| 
234
{| style="text-align: left; margin:auto;width: 100%;" 
235
|-
236
| style="text-align: center;" | <math>\mathbf{T}_{L}    =\mathbf{R}_{L}^{T}\;\mathbf{T\;R}_{L}</math>
237
|-
238
| style="text-align: center;" | <math> \mathbf{S}_{L}    =\mathbf{R}_{L}^{T}\;\mathbf{S\;R}_{L} </math>
239
|}
240
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
241
|}
242
243
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
244
245
{| class="formulaSCP" style="width: 100%; text-align: left;" 
246
|-
247
| 
248
{| style="text-align: left; margin:auto;width: 100%;" 
249
|-
250
| 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>
251
|}
252
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
253
|}
254
255
The relationship between the Hencky and Piola-Kirchhoff stresses is
256
257
{| class="formulaSCP" style="width: 100%; text-align: left;" 
258
|-
259
| 
260
{| style="text-align: left; margin:auto;width: 100%;" 
261
|-
262
| style="text-align: center;" | <math>\left[ S_{L}\right] _{\alpha \alpha }    =\frac{1}{\lambda _{\alpha }^{2} }\left[ T_{L}\right] _{\alpha \alpha }</math>
263
|-
264
| style="text-align: center;" | <math> \left[ S_{L}\right] _{\alpha \beta }    =\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>
265
|}
266
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
267
|}
268
269
The second Piola-Kirchhoff stress tensor can be finally computed by
270
271
{| class="formulaSCP" style="width: 100%; text-align: left;" 
272
|-
273
| 
274
{| style="text-align: left; margin:auto;width: 100%;" 
275
|-
276
| style="text-align: center;" | <math>\mathbf{S}=\mathbf{R}_{L}\;\mathbf{S}_{L}\mathbf{\;R}_{L}^{T} </math>
277
|}
278
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
279
|}
280
281
The generalized stresses (forces and moments) are obtained by integrating across the  original thickness the second Piola-Kirchhoff stress tensor using the actual  distance to the middle surface for the evaluation of the bending moments, i.e.
282
283
{| class="formulaSCP" style="width: 100%; text-align: left;" 
284
|-
285
| 
286
{| style="text-align: left; margin:auto;width: 100%;" 
287
|-
288
| style="text-align: center;" | <math>\mathbf{N}=\int _{h^{0}}\mathbf{S} d\zeta </math>
289
| style="width: 5px;text-align: right;white-space: nowrap;" | (22.a)
290
|-
291
| style="text-align: center;" | <math> \mathbf{M}=\int _{h^{0}}\mathbf{S }\lambda \zeta  d\zeta  </math>
292
| style="width: 5px;text-align: right;white-space: nowrap;" | (22.b)
293
|}
294
|}
295
296
With these values the weak form of the equilibrium equations can be written as
297
298
<span id="eq-23"></span>
299
{| class="formulaSCP" style="width: 100%; text-align: left;" 
300
|-
301
| 
302
{| style="text-align: left; margin:auto;width: 100%;" 
303
|-
304
| style="text-align: center;" | <math>\delta \Pi =\int _{\Omega ^{0}}\left[ \delta \mathbf{E}:\mathbf{N} +\delta {\boldsymbol K}:\mathbf{M}\right] d\Omega ^{0}+\delta \Pi _{ext}=0   </math>
305
|}
306
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
307
|}
308
309
where <math display="inline">{\boldsymbol K}</math> is the virtual curvature vector and <math display="inline">\mathbf{E}</math> is the virtual  Green-Lagrange strain vector of the middle surface, with
310
311
{| class="formulaSCP" style="width: 100%; text-align: left;" 
312
|-
313
| 
314
{| style="text-align: left; margin:auto;width: 100%;" 
315
|-
316
| style="text-align: center;" | <math>E_{ij}={1\over 2} (a_{ij}-\delta _{ij})  </math>
317
|}
318
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
319
|}
320
321
==3 Constitutive models==
322
323
In this section, a brief description  of the constitutive models used  in the numerical examples presented in Sections 7 and 8 is given. Two types of  materials are described: an elastic-plastic material associated to thin rolled  metal sheets and a hyper-elastic material for rubbers.
324
325
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
326
327
{| class="formulaSCP" style="width: 100%; text-align: left;" 
328
|-
329
| 
330
{| style="text-align: left; margin:auto;width: 100%;" 
331
|-
332
| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=E}_{\ln }^{e}+\mathbf{E}_{\ln }^{p} </math>
333
|}
334
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
335
|}
336
337
A constant linear relationship between (plane) stresses and elastic strains is  also adopted giving
338
339
{| class="formulaSCP" style="width: 100%; text-align: left;" 
340
|-
341
| 
342
{| style="text-align: left; margin:auto;width: 100%;" 
343
|-
344
| style="text-align: center;" | <math>\mathbf{T}=\mathbf{CE}_{\ln }^{e} </math>
345
|}
346
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
347
|}
348
349
These constitutive equations are integrated using a standard return algorithm.  The following Mises-Hill <span id='citeF-11'></span>[[#cite-11|[11]]] yield function with non-linear isotropic  hardening is adopted
350
351
{| class="formulaSCP" style="width: 100%; text-align: left;" 
352
|-
353
| 
354
{| style="text-align: left; margin:auto;width: 100%;" 
355
|-
356
| 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}=\kappa \left( \varepsilon _{0}+e^{p}\right) ^{n} </math>
357
|}
358
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
359
|}
360
361
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">\kappa ,</math> <math display="inline">\varepsilon _{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>.
362
363
The Mises-Hill yield function has the advantage of simplicity and it allows, as a  first approximation, to treat rolled thin metal sheets with  planar and transversal anisotropy.
364
365
For the case of rubbers, the Ogden <span id='citeF-12'></span>[[#cite-12|[12]]] model extended to  the compressible range is considered. The material behaviour is characterized  by the strain energy density per unit undeformed volume of the form
366
367
{| class="formulaSCP" style="width: 100%; text-align: left;" 
368
|-
369
| 
370
{| style="text-align: left; margin:auto;width: 100%;" 
371
|-
372
| style="text-align: center;" | <math>\psi =\frac{K}{2}\left( \ln J\right) ^{2}+\sum _{p=1}^{N}\frac{\mu _{p}}{\alpha _{p}}\left[ J^{-\dfrac{\alpha _{p}}{3}}\left( \sum _{i=1}^{3} \lambda _{i}^{\alpha _{p}-1}\right) -3\right] </math>
373
|}
374
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
375
|}
376
377
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.
378
379
The stress measures associated to the principal logarithmic strains are  denoted by <math display="inline">\beta _{i}</math>. They can be computed noting that
380
381
{| class="formulaSCP" style="width: 100%; text-align: left;" 
382
|-
383
| 
384
{| style="text-align: left; margin:auto;width: 100%;" 
385
|-
386
| style="text-align: center;" | <math>\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^{-\dfrac{\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>
387
|}
388
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
389
|}
390
391
we define now
392
393
{| class="formulaSCP" style="width: 100%; text-align: left;" 
394
|-
395
| 
396
{| style="text-align: left; margin:auto;width: 100%;" 
397
|-
398
| style="text-align: center;" | <math>a_{p}=\sum _{j=1}^{3}\lambda _{j}^{\alpha _{p}} </math>
399
|}
400
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
401
|}
402
403
which gives
404
405
{| class="formulaSCP" style="width: 100%; text-align: left;" 
406
|-
407
| 
408
{| style="text-align: left; margin:auto;width: 100%;" 
409
|-
410
| style="text-align: center;" | <math>\beta _{i}=K\left( \ln J\right) +\sum _{p=1}^{N}\mu _{p}J^{-\dfrac{\alpha _{p} }{3}}\left( \lambda _{i}^{\alpha _{p}}-\frac{1}{3}a_{p}\right) </math>
411
|}
412
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
413
|}
414
415
The values of <math display="inline">\beta _i</math>, expressed in the principal strains directions, allow to evaluate  the stresses  in the convective coordinate system as
416
417
{| class="formulaSCP" style="width: 100%; text-align: left;" 
418
|-
419
| 
420
{| style="text-align: left; margin:auto;width: 100%;" 
421
|-
422
| style="text-align: center;" | <math>\mathbf{T}=\sum _{i=1}^{3}\beta _{i}\;\mathbf{r}_{i}\otimes \mathbf{r}_{i} </math>
423
|}
424
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
425
|}
426
427
We note that the Hencky stress tensor <math display="inline">\mathbf{T}</math> can be easily  particularized for the plane stress case.
428
429
==4 Interpolation functions and gradient evaluation==
430
431
The starting point of the enhanced rotation-free of BST element is  to discretize the shell surface with a standard 3-node  triangular mesh. The difference with a standard finite element method is that for the  computation of strains within an element, the configuration of the three adjacent triangular  elements is  used. Then at each triangle a patch of four elements, formed by the central triangle and the three adjacent ones is considered (see Figure [[#img-1|1]].a). This allows to define a  quadratic interpolation of the geometry from the position  of the 6 nodes. In the isoparametric space, we keep the vertices of the  3-node master element (standard linear triangle) with the positions
432
433
* node 1: <math display="inline">\left( \xi ,\eta \right) =\left( 0,0\right) </math>
434
435
* node 2: <math display="inline">\left( \xi ,\eta \right) =\left( 1,0\right) </math>
436
437
* node 3: <math display="inline">\left( \xi ,\eta \right) =\left( 0,1\right) </math>
438
439
and the three additional nodes that complete the patch with the positions
440
441
* node 4: <math display="inline">\left( \xi ,\eta \right) =\left( 1,1\right) </math>
442
443
* node 5: <math display="inline">\left( \xi ,\eta \right) =\left( -1,1\right) </math>
444
445
* node 6: <math display="inline">\left( \xi ,\eta \right) =\left( 1,-1\right) </math>
446
447
<div id='img-1'></div>
448
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
449
|-
450
|[[Image:Draft_Samper_340530075-test-fig1.png|600px|Patch of four elements  (a)in spacial coordinates (b)in natural coordinates]]
451
|- style="text-align: center; font-size: 75%;"
452
| colspan="1" | '''Figure 1:''' Patch of four elements  (a)in spacial coordinates (b)in natural coordinates
453
|}
454
455
It is possible to define now a set of (non-standard) quadratic shape functions over  over the six node element as:
456
457
{| class="formulaSCP" style="width: 100%; text-align: left;" 
458
|-
459
| 
460
{| style="text-align: left; margin:auto;width: 100%;" 
461
|-
462
| style="text-align: center;" | <math>\begin{array}{ccc} N^{1}=\zeta{+\xi}\eta &  & N^{4}=\frac{\zeta }{2}\left( \zeta{-1}\right)\\
463
 N^{2}=\xi{+\eta}\zeta &  & N^{5}=\frac{\xi }{2}\left( \xi{-1}\right)\\
464
 N^{3}=\eta{+\zeta}\xi &  & N^{6}=\frac{\eta }{2}\left( \eta{-1}\right) \end{array}  </math>
465
|}
466
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
467
|}
468
469
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 three points over the boundary of the central  element of the patch. These points are located at the mid-point of each side,  denoted by G1, G2 and G3 in Figure [[#img-1|1]].b. This election has the  following characteristics.
470
471
* Gradients at these points depend only on the nodes pertaining to  the two elements adjacent to the side. This can be easily verified by sampling the derivatives of  the shape functions  at each mid-side point.
472
473
* When gradients are computed at the common mid-side point from two  adjacent elements, the same values are  obtained, as the coordinates of  the same four points are used. That is, the gradients at the mid-side points are  independent of the element where they are computed. A side-oriented  implementation of the finite element could lead to a unique evaluation of the  gradients per side.
474
475
Next, some details of the implementation chosen are presented. We denote by <math display="inline">\mathbf{t}_{1}</math> and <math display="inline">\mathbf{t}_{2}</math> the orthogonal unit vectors over the  tangent plane (undeformed configuration) associated to a conveniently selected  local cartesian system. The cartesian derivatives of the shape functions are computed at the original configuration by the standard expression
476
477
{| class="formulaSCP" style="width: 100%; text-align: left;" 
478
|-
479
| 
480
{| style="text-align: left; margin:auto;width: 100%;" 
481
|-
482
| style="text-align: center;" | <math>\left[ \begin{array}{c} N_{^{\prime }1}^{I}\\
483
 N_{^{\prime }2}^{I} \end{array}  \right] =\mathbf{J}^{-1}\left[ \begin{array}{c} N_{^{\prime }\xi }^{I}\\
484
 N_{^{\prime }\eta }^{I} \end{array}  \right] </math>
485
|}
486
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
487
|}
488
489
where the Jacobian matrix at the original configuration is
490
491
{| class="formulaSCP" style="width: 100%; text-align: left;" 
492
|-
493
| 
494
{| style="text-align: left; margin:auto;width: 100%;" 
495
|-
496
| style="text-align: center;" | <math>\mathbf{J=}\left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }\xi }^{\left( 0\right) }\cdot \mathbf{t}_{1} &  \mathbf{\varphi }_{^{\prime }\eta }^{\left( 0\right) }\cdot \mathbf{t}_{1}\\
497
 \mathbf{\varphi }_{^{\prime }\xi }^{\left( 0\right) }\cdot \mathbf{t}_{2} &  \mathbf{\varphi }_{^{\prime }\eta }^{\left( 0\right) }\cdot \mathbf{t}_{2} \end{array}  \right] </math>
498
|}
499
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
500
|}
501
502
With the previous definitions the deformation gradient on the middle surface,  associated to an arbitrary spatial cartesian system and to the material  cartesian system defined on the middle surface, is
503
504
{| class="formulaSCP" style="width: 100%; text-align: left;" 
505
|-
506
| 
507
{| style="text-align: left; margin:auto;width: 100%;" 
508
|-
509
| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1},\mathbf{\varphi }_{^{\prime }2}\right] =\left[ \mathbf{\varphi }_{^{\prime }\xi },\mathbf{\varphi }_{^{\prime }\eta  }\right]  \mathbf{J}^{-1} </math>
510
|}
511
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
512
|}
513
514
The covariant metric tensor can be computed from above expression as
515
516
{| class="formulaSCP" style="width: 100%; text-align: left;" 
517
|-
518
| 
519
{| style="text-align: left; margin:auto;width: 100%;" 
520
|-
521
| style="text-align: center;" | <math>\mathbf{g=}\left[ \begin{array}{cc} a_{11} & a_{12}\\
522
 a_{21} & a_{22} \end{array}  \right] =\left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2}\\
523
 \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] </math>
524
|}
525
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
526
|}
527
528
The component of <math display="inline">\boldsymbol g</math> can be used to compute any convenient membrane strain measures. For example the  Green-Lagrange strain tensor is obtained by
529
530
{| class="formulaSCP" style="width: 100%; text-align: left;" 
531
|-
532
| 
533
{| style="text-align: left; margin:auto;width: 100%;" 
534
|-
535
| style="text-align: center;" | <math>\mathbf{E}_{GL}\mathbf{=}\frac{1}{2}\left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1}-1 &  \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2}\\
536
 \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }1} &  \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2}-1  \end{array}  \right] =\frac{1}{2}\left[ \begin{array}{cc} a_{11}-1 & a_{12}\\
537
 a_{21} & a_{22}-1  \end{array}  \right] </math>
538
|}
539
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
540
|}
541
542
In the usual FEM matrix notation
543
544
{| class="formulaSCP" style="width: 100%; text-align: left;" 
545
|-
546
| 
547
{| style="text-align: left; margin:auto;width: 100%;" 
548
|-
549
| style="text-align: center;" | <math>\left[ \begin{array}{c} E_{11}\\
550
 E_{22}\\
551
 2E_{12} \end{array}  \right] =\frac{1}{2}\left[ \begin{array}{c} \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }1}-1\\
552
 \mathbf{\varphi }_{^{\prime }2}\cdot \mathbf{\varphi }_{^{\prime }2}-1\\
553
 2\mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] </math>
554
|}
555
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
556
|}
557
558
The virtual strains are obtained by the variation of above expression as
559
560
{| class="formulaSCP" style="width: 100%; text-align: left;" 
561
|-
562
| 
563
{| style="text-align: left; margin:auto;width: 100%;" 
564
|-
565
| style="text-align: center;" | <math>\delta \left[ \begin{array}{c} E_{11}\\
566
 E_{22}\\
567
 2E_{12} \end{array}  \right] =\left[ \begin{array}{c} \mathbf{\varphi }_{^{\prime }1}\cdot \delta \mathbf{\varphi }_{^{\prime }1}\\
568
 \mathbf{\varphi }_{^{\prime }2}\cdot \delta \mathbf{\varphi }_{^{\prime }2}\\
569
 \mathbf{\varphi }_{^{\prime }1}\cdot \delta \mathbf{\varphi }_{^{\prime }2} +\delta \mathbf{\varphi }_{^{\prime }1}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] </math>
570
|}
571
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
572
|}
573
574
In sides of elements where the adjacent element does not exist  (i.e. at the  boundaries of the shell), the gradient on that side is made equal to the  gradient computed over the 3 nodes of the central element of the patch. From  the definition of the gradients at the three mid-side points of the triangle,  the element formulation can be interpreted as an ''assumed strain'' approach, where  the metric tensor is interpolated over the element using the values computed  on the sides by
575
576
{| class="formulaSCP" style="width: 100%; text-align: left;" 
577
|-
578
| 
579
{| style="text-align: left; margin:auto;width: 100%;" 
580
|-
581
| style="text-align: center;" | <math>\mathbf{g}\left( \xi ,\eta \right) =\left( 1-2\zeta \right) \mathbf{g}^{1}  +\left( 1-2\xi \right) \mathbf{g}^{2} +\left( 1-2\eta \right) \mathbf{g}^{3} </math>
582
|}
583
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
584
|}
585
586
==5 Computation of curvatures==
587
588
Curvatures are assumed to be constant over the element. They are computed by performing  an average over the central element as
589
590
{| class="formulaSCP" style="width: 100%; text-align: left;" 
591
|-
592
| 
593
{| style="text-align: left; margin:auto;width: 100%;" 
594
|-
595
| style="text-align: center;" | <math>\kappa _{\alpha \beta }\mathbf{ =}\frac{-1}{A^{M}}\int _{A^{M}}\mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }\beta \alpha }\;dA^{M} </math>
596
|}
597
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
598
|}
599
600
where (<math display="inline">A^{M}</math> is the original area of the element).
601
602
Integrating by parts Eq.(43), the following integral over the element boundary <math display="inline">\Gamma ^{M}</math> is  obtained:
603
604
{| class="formulaSCP" style="width: 100%; text-align: left;" 
605
|-
606
| 
607
{| style="text-align: left; margin:auto;width: 100%;" 
608
|-
609
| style="text-align: center;" | <math>\kappa _{\alpha \beta }\mathbf{ =}\frac{-1}{A^{M}}{\displaystyle \oint  _{\Gamma ^{M}}}\mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }\alpha }\;n_{\beta  }\;d\Gamma ^{M} </math>
610
|}
611
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
612
|}
613
614
The three distinct components of the curvature tensor are
615
616
<span id="eq-44"></span>
617
{| class="formulaSCP" style="width: 100%; text-align: left;" 
618
|-
619
| 
620
{| style="text-align: left; margin:auto;width: 100%;" 
621
|-
622
| style="text-align: center;" | <math>\left[ \begin{array}{c} \kappa _{11}\\
623
 \kappa _{22}\\
624
 2\kappa _{12} \end{array}  \right] =\frac{-1}{A^{M}}{\displaystyle \oint _{\Gamma ^{M}}}\left[ \begin{array}{cc} n_{1} & 0\\
625
 0 & n_{2}\\
626
 n_{2} & n_{1} \end{array}  \right] \left[ \begin{array}{c} \mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }1}\\
627
 \mathbf{t}_{3}\cdot \mathbf{\varphi }_{^{\prime }2} \end{array}  \right] d\Gamma ^{M}  </math>
628
|}
629
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
630
|}
631
632
The numerical evaluation of the boundary integral in Eq.([[#eq-44|44]]) results in a sum over the integration  points on the element boundary which are, in fact, the same points used for the computation of  the gradients in the membrane formulations, i.e.
633
634
{| class="formulaSCP" style="width: 100%; text-align: left;" 
635
|-
636
| 
637
{| style="text-align: left; margin:auto;width: 100%;" 
638
|-
639
| style="text-align: center;" | <math>\left[ \begin{array}{c} \kappa _{11}\\
640
 \kappa _{22}\\
641
 2\kappa _{12} \end{array}  \right] =\frac{-1}{A^{M}}\sum _{G=1}^{3}l_{G}\left[ \begin{array}{cc} n_{1} & 0\\
642
 0 & n_{2}\\
643
 n_{2} & n_{1} \end{array}  \right] ^{G}\left[ \begin{array}{c} \mathbf{\varphi }_{1}^{G}\cdot \mathbf{t}_{3}\\
644
 \mathbf{\varphi }_{2}^{G}\cdot \mathbf{t}_{3} \end{array}  \right] </math>
645
|}
646
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
647
|}
648
649
For the definition of <math display="inline">\mathbf{t}_{3}</math>, the linear interpolation over the  central element is used. In this case the tangent plane components are (with <math display="inline">a_{i}</math>  and <math display="inline">b_{i}</math> being the usual cartesian projections of the sides when area coordinates  are used)
650
651
<span id="eq-46"></span>
652
{| class="formulaSCP" style="width: 100%; text-align: left;" 
653
|-
654
| 
655
{| style="text-align: left; margin:auto;width: 100%;" 
656
|-
657
| style="text-align: center;" | <math>\left[ \begin{array}{c} \mathbf{\varphi }_{1}\\
658
 \mathbf{\varphi }_{2} \end{array}  \right] ^{M}    =\frac{1}{2A^{M}}\sum _{i=1}^{3}\left[ \begin{array}{c} -b_{i}\\
659
 a_{i} \end{array}  \right] \mathbf{\varphi }^{i}=\left[ \begin{array}{ccc} \bar{N}_{^{\prime }1}^{1} & \bar{N}_{^{\prime }1}^{2} & \bar{N}_{^{\prime }1} ^{3}\\
660
 \bar{N}_{^{\prime }2}^{1} & \bar{N}_{^{\prime }2}^{2} & \bar{N}_{^{\prime }2}^{3} \end{array}  \right] \left[ \begin{array}{c} \mathbf{\varphi }^{1}\\
661
 \mathbf{\varphi }^{2}\\
662
 \mathbf{\varphi }^{3} \end{array}  \right]  </math>
663
|}
664
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
665
|}
666
667
<span id="eq-47"></span>
668
{| class="formulaSCP" style="width: 100%; text-align: left;" 
669
|-
670
| 
671
{| style="text-align: left; margin:auto;width: 100%;" 
672
|-
673
| style="text-align: center;" | <math>\mathbf{t}_{3}    =\frac{\mathbf{\varphi }_{1}^{M}\times \mathbf{\varphi }_{2}^{M}}{\left\vert \mathbf{\varphi }_{1}^{M}\times \mathbf{\varphi }_{2} ^{M}\right\vert }=\lambda \;\mathbf{\varphi }_{1}^{M}\times \mathbf{\varphi }_{2}^{M}  </math>
674
|}
675
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
676
|}
677
678
From these expressions it is also possible to compute in the original  configuration the element area <math display="inline">A^{M}</math>, the outer normals <math display="inline">\left( n_{1} ,n_{2}\right) ^{i}</math> at each side and its lengths <math display="inline">l_{i}</math>. Eq.([[#eq-47|47]]) 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>. In Eq.([[#eq-46|46]]), <math display="inline">\bar{N}^{I}</math> are the linear shape functions of the three node linear triangle <span id='citeF-3'></span>[[#cite-3|[3]]].
679
680
As one integration point is used over each side, it is not necessary to  distinguish  between sides (<math display="inline">I</math>) and integrations points (<math display="inline">G</math>) any more. In  this way the curvatures can be computed by:
681
682
{| class="formulaSCP" style="width: 100%; text-align: left;" 
683
|-
684
| 
685
{| style="text-align: left; margin:auto;width: 100%;" 
686
|-
687
| style="text-align: center;" | <math>\left[ \begin{array}{c} \kappa _{11}\\
688
 \kappa _{22}\\
689
 2\kappa _{12} \end{array}  \right]    =\frac{-1}{A^{M}}\sum _{I=1}^{3}l_{I}\left[ \begin{array}{cc}[ n_{1} & 0\\
690
 0 & n_{2}\\
691
 n_{2} & n_{1} \end{array}  \right] ^{G}\left[ \begin{array}{c} \mathbf{\varphi }_{1}^{I}\cdot \mathbf{t}_{3}\\
692
 \mathbf{\varphi }_{2}^{I}\cdot \mathbf{t}_{3} \end{array}  \right]</math>
693
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
694
|-
695
| style="text-align: center;" | <math>   =2\sum _{I=1}^{3}\left[ \begin{array}{cc} \bar{N}_{^{\prime }1}^{I} & 0\\
696
 0 & \bar{N}_{^{\prime }2}^{I}\\
697
 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left[ \begin{array}{c} \mathbf{\varphi }_{1}^{I}\cdot \mathbf{t}_{3}\\
698
 \mathbf{\varphi }_{2}^{I}\cdot \mathbf{t}_{3} \end{array}  \right] </math>
699
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
700
|}
701
|}
702
703
In the original BST element <span id='citeF-4'></span><span id='citeF-10'></span>[[#cite-4|[4,10]]] the gradient <math display="inline">\mathbf{\varphi  }_{\alpha }^{I}</math> was computed as the average of the linear approximations over  the two adjacent elements. In the new version, the gradient is  computed at each side <math display="inline">I</math> from the quadratic interpolation
704
705
{| class="formulaSCP" style="width: 100%; text-align: left;" 
706
|-
707
| 
708
{| style="text-align: left; margin:auto;width: 100%;" 
709
|-
710
| style="text-align: center;" | <math>\left[ \begin{array}{c} \mathbf{\varphi }_{1}\\
711
 \mathbf{\varphi }_{2} \end{array}  \right] ^{I}=\left[ \begin{array}{cccc} N_{^{\prime }1}^{1} & N_{^{\prime }1}^{2} & N_{^{\prime }1}^{3} & N_{^{\prime } 1}^{I+3}\\
712
 N_{^{\prime }2}^{1} & N_{^{\prime }2}^{2} & N_{^{\prime }2}^{3} & N_{^{\prime } 2}^{I+3} \end{array}  \right] ^{I}\left[ \begin{array}{c} \mathbf{\varphi }^{1}\\
713
 \mathbf{\varphi }^{2}\\
714
 \mathbf{\varphi }^{3}\\
715
 \mathbf{\varphi }^{I+3} \end{array}  \right] </math>
716
|}
717
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
718
|}
719
720
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">I</math>) where the gradient is computed.
721
722
An alternative form to express the curvatures, which is useful when their  variations are needed, is to define:
723
724
<span id="eq-51"></span>
725
{| class="formulaSCP" style="width: 100%; text-align: left;" 
726
|-
727
| 
728
{| style="text-align: left; margin:auto;width: 100%;" 
729
|-
730
| style="text-align: center;" | <math>\mathbf{h}_{ij}=\sum _{I=1}^{3}\left( \bar{N}_{^{\prime }i}^{I} \;\mathbf{\varphi }_{j}^{I}+\bar{N}_{^{\prime }j}^{I}\;\mathbf{\varphi }_{i} ^{I}\right)  </math>
731
|}
732
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
733
|}
734
735
This gives
736
737
<span id="eq-52"></span>
738
{| class="formulaSCP" style="width: 100%; text-align: left;" 
739
|-
740
| 
741
{| style="text-align: left; margin:auto;width: 100%;" 
742
|-
743
| style="text-align: center;" | <math>\kappa _{ij}=\mathbf{h}_{ij}\cdot \mathbf{t}_{3}  </math>
744
|}
745
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
746
|}
747
748
This 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. Direction  '''t'''<math display="inline">_{3}</math> can be seen as a reference direction.  If a different  direction  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-47|47]]) for the  definition of '''t'''<math display="inline">_{3}</math> as a function exclusively of the three nodes of  the central triangle,  instead of the 6-node  isoparametric interpolation.
749
750
===5.1 Variation of the curvatures===
751
752
From Eq.([[#eq-52|52]]) the variation of the curvatures is obtained as
753
754
<span id="eq-53"></span>
755
{| class="formulaSCP" style="width: 100%; text-align: left;" 
756
|-
757
| 
758
{| style="text-align: left; margin:auto;width: 100%;" 
759
|-
760
| style="text-align: center;" | <math>\delta \kappa _{ij}=\delta \mathbf{h}_{ij}\cdot \mathbf{t}_{3}+\mathbf{h}_{ij}\cdot \delta \mathbf{t}_{3}  </math>
761
|}
762
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
763
|}
764
765
The first term is crucial in the evaluation of the curvatures variation. The variations of <math display="inline">\mathbf{h}_{ij}</math> are:
766
767
<span id="eq-54"></span>
768
{| class="formulaSCP" style="width: 100%; text-align: left;" 
769
|-
770
| 
771
{| style="text-align: left; margin:auto;width: 100%;" 
772
|-
773
| style="text-align: center;" | <math>\delta \mathbf{h}_{ij}=\sum _{I=1}^{3}\left( \bar{N}_{^{\prime }i}^{I} \;\delta \mathbf{\varphi }_{j}^{I}+\bar{N}_{^{\prime }j}^{I}\;\delta  \mathbf{\varphi }_{i}^{I}\right)  </math>
774
|}
775
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
776
|}
777
778
where
779
780
{| class="formulaSCP" style="width: 100%; text-align: left;" 
781
|-
782
| 
783
{| style="text-align: left; margin:auto;width: 100%;" 
784
|-
785
| style="text-align: center;" | <math>\left[ \begin{array}{c} \delta \mathbf{\varphi }_{1}\\
786
 \delta \mathbf{\varphi }_{2} \end{array}  \right] ^{I}=\left[ \begin{array}{cccc} N_{^{\prime }1}^{1} & N_{^{\prime }1}^{2} & N_{^{\prime }1}^{3} & N_{^{\prime } 1}^{I+3}\\
787
 N_{^{\prime }2}^{1} & N_{^{\prime }2}^{2} & N_{^{\prime }2}^{3} & N_{^{\prime } 2}^{I+3} \end{array}  \right] ^{I}\left[ \begin{array}{c} \delta \mathbf{u}^{1}\\
788
 \delta \mathbf{u}^{2}\\
789
 \delta \mathbf{u}^{3}\\
790
 \delta \mathbf{u}^{I+3} \end{array}  \right] </math>
791
|}
792
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
793
|}
794
795
The variations of '''t'''<math display="inline">_{3}</math> can be computed as shown in the  original BST element <span id='citeF-10'></span>[[#cite-10|[10]]].
796
797
{| class="formulaSCP" style="width: 100%; text-align: left;" 
798
|-
799
| 
800
{| style="text-align: left; margin:auto;width: 100%;" 
801
|-
802
| style="text-align: center;" | <math>\begin{array}{l}\delta t_{31}    =-\mathbf{t}_{3}\cdot \delta \mathbf{\varphi }_{^{\prime }1} ^{M}\\
803
 \\
804
 \delta t_{32}    =-\mathbf{t}_{3}\cdot \delta \mathbf{\varphi }_{^{\prime }2}^{M} \end{array}  </math>
805
|}
806
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
807
|}
808
809
leading to
810
811
{| class="formulaSCP" style="width: 100%; text-align: left;" 
812
|-
813
| 
814
{| style="text-align: left; margin:auto;width: 100%;" 
815
|-
816
| style="text-align: center;" | <math>\delta \mathbf{t}_{3}=\delta t_{31}\mathbf{\tilde{\varphi }}_{^{\prime }1}+\delta  t_{32}\mathbf{\tilde{\varphi }}_{^{\prime }2} </math>
817
|}
818
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
819
|}
820
821
where <math display="inline">\mathbf{\tilde{\varphi }}_{^{\prime }\alpha }</math> are the  contravariant base vectors in the central triangle
822
823
{| class="formulaSCP" style="width: 100%; text-align: left;" 
824
|-
825
| 
826
{| style="text-align: left; margin:auto;width: 100%;" 
827
|-
828
| style="text-align: center;" | <math>\begin{array}{l}\mathbf{\tilde{\varphi }}_{^{\prime }1}    =\lambda \;\mathbf{\varphi }_{^{\prime }2}^{M}\times \mathbf{t}_{3}\\
829
 \\
830
 \mathbf{\tilde{\varphi }}_{^{\prime }2}    =-\lambda \;\mathbf{\varphi  }_{^{\prime }1}^{M}\times \mathbf{t}_{3} \end{array}  </math>
831
|}
832
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
833
|}
834
835
We then have
836
837
{| class="formulaSCP" style="width: 100%; text-align: left;" 
838
|-
839
| 
840
{| style="text-align: left; margin:auto;width: 100%;" 
841
|-
842
| style="text-align: center;" | <math>\delta \mathbf{t}_{3}    =\left( -\mathbf{t}_{3}\cdot \delta \mathbf{\varphi  }_{^{\prime }1}^{M}\right) \mathbf{\tilde{\varphi }}_{^{\prime }1}+\left( -\mathbf{t}_{3}\cdot \delta \mathbf{\varphi }_{^{\prime }2}^{M}\right) \mathbf{\tilde{\varphi }}_{^{\prime }2}</math>
843
|-
844
| style="text-align: center;" | <math>   =-\sum _{J=1}^{3}\left[ \bar{N}_{^{\prime }1}^{J}\mathbf{\tilde{\varphi }}_{^{\prime }1}+\bar{N}_{^{\prime }2}^{J}\mathbf{\tilde{\varphi }}_{^{\prime } 2}\right] \left( \mathbf{t}_{3}\cdot \delta \mathbf{u}^{J}\right) </math>
845
|}
846
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
847
|}
848
849
Substituting the last expression in Eq.([[#eq-53|53]]), results
850
851
<span id="eq-60"></span>
852
{| class="formulaSCP" style="width: 100%; text-align: left;" 
853
|-
854
| 
855
{| style="text-align: left; margin:auto;width: 100%;" 
856
|-
857
| style="text-align: center;" | <math>\left[ \begin{array}{c} \delta \chi _{11}\\
858
 \delta \chi _{22}\\
859
 2\delta \chi _{12} \end{array}  \right]    =2\sum \limits _{I=1}^{3}\left[ \begin{array}{cc} \bar{N}_{^{\prime }1}^{I} & 0\\
860
 0 & \bar{N}_{^{\prime }2}^{I}\\
861
 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left\{ \sum \limits _{J=1}^{3}\left[ \begin{array}{c} N_{^{\prime }1}^{J\left( I\right) }\left( \mathbf{t}_{3}\cdot \delta  \mathbf{u}^{J}\right)\\
862
 N_{^{\prime }2}^{J\left( I\right) }\left( \mathbf{t}_{3}\cdot \delta  \mathbf{u}^{J}\right) \end{array}  \right] +\left[ \begin{array}{c} N_{^{\prime }1}^{I+3\left( I\right) }\left( \mathbf{t}_{3}\cdot  \delta \mathbf{u}^{I+3}\right)\\
863
 N_{^{\prime }2}^{I+3\left( I\right) }\left( \mathbf{t}_{3}\cdot  \delta \mathbf{u}^{I+3}\right) \end{array}  \right] \right\}</math>
864
|-
865
| style="text-align: center;" | <math>   -2\sum \limits _{I=1}^{3}\left[ \begin{array}{c} \left( \bar{N}_{^{\prime }1}^{I}\rho _{11}^{1}+\bar{N}_{^{\prime }2}^{I} \rho _{11}^{2}\right)\\
866
 \left( \bar{N}_{^{\prime }1}^{I}\rho _{22}^{1}+\bar{N}_{^{\prime }2}^{I} \rho _{22}^{2}\right)\\
867
 \left( \bar{N}_{^{\prime }1}^{I}\rho _{12}^{1}+\bar{N}_{^{\prime }2}^{I} \rho _{12}^{2}\right) \end{array}  \right] \left( \mathbf{t}_{3}\cdot \delta \mathbf{u}^{I}\right)  </math>
868
|}
869
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
870
|}
871
872
where the projections of the vectors <math display="inline">\mathbf{h}_{ij}</math> over the  contravariant base vectors <math display="inline">\mathbf{\tilde{\varphi }}_{^{\prime }\alpha }</math> have  been included
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;width: 100%;" 
879
|-
880
| style="text-align: center;" | <math>\rho _{ij}^{\alpha }=\mathbf{h}_{ij}\cdot \mathbf{\tilde{\varphi }}_{^{\prime  }\alpha }  </math>
881
|}
882
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
883
|}
884
885
These projections together with Eq.([[#eq-52|52]]) allows to write
886
887
{| class="formulaSCP" style="width: 100%; text-align: left;" 
888
|-
889
| 
890
{| style="text-align: left; margin:auto;width: 100%;" 
891
|-
892
| style="text-align: center;" | <math>\mathbf{h}_{ij}=\sum _{\alpha=1}^{2}\rho _{ij}^{\alpha }\;\mathbf{\varphi  }_{^{\prime }\alpha }+\kappa _{ij}\mathbf{t}_{3} </math>
893
|}
894
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
895
|}
896
897
===5.2 Boundary conditions===
898
899
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 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">\mathbf{\varphi  }_{^{\prime }n}^{\left( 0\right) }</math> that does not change during the process.
900
901
<div id='img-2'></div>
902
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
903
|-
904
|[[Image:Draft_Samper_340530075-test-fig2.png|300px|Local cartesian system for the treatment of boundary conditions]]
905
|- style="text-align: center; font-size: 75%;"
906
| colspan="1" | '''Figure 2:''' Local cartesian system for the treatment of boundary conditions
907
|}
908
909
The tangent plane at the boundary is expressed in terms of two orthogonal  unit vectors referred to as local-to-the-boundary cartesian system (see Figure  [[#img-2|2]]) and defined as
910
911
{| class="formulaSCP" style="width: 100%; text-align: left;" 
912
|-
913
| 
914
{| style="text-align: left; margin:auto;width: 100%;" 
915
|-
916
| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }n}^{\left( 0\right) },\;\mathbf{\bar  {\varphi }}_{^{\prime }s}\right] </math>
917
|}
918
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
919
|}
920
921
where vector <math display="inline">\mathbf{\varphi }_{^{\prime }n}^{\left( 0\right) }</math> is fixed while  direction <math display="inline">\mathbf{\bar{\varphi }}_{^{\prime }s}</math> emerges as 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 is
922
923
{| class="formulaSCP" style="width: 100%; text-align: left;" 
924
|-
925
| 
926
{| style="text-align: left; margin:auto;width: 100%;" 
927
|-
928
| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1}^{M},\;\mathbf{\varphi }_{^{\prime }2} ^{M}\right] </math>
929
|}
930
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
931
|}
932
933
the intersection 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_{0}</math>, i.e.
934
935
{| class="formulaSCP" style="width: 100%; text-align: left;" 
936
|-
937
| 
938
{| style="text-align: left; margin:auto;width: 100%;" 
939
|-
940
| style="text-align: center;" | <math>\mathbf{\varphi }_{^{\prime }s}=\frac{1}{L_{0}}\left( \mathbf{\varphi }^{K}-\mathbf{\varphi }^{J}\right) </math>
941
|}
942
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
943
|}
944
945
That together with the outer normal to the side <math display="inline">\mathbf{n =}\left[ \nu  _{1},\nu _{2}\right]^T  </math>(resolved in the selected original convective cartesian  system) leads to
946
947
{| class="formulaSCP" style="width: 100%; text-align: left;" 
948
|-
949
| 
950
{| style="text-align: left; margin:auto;width: 100%;" 
951
|-
952
| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1}^{I},\;\mathbf{\varphi }_{^{\prime }2} ^{I}\right] =\left[ \mathbf{\varphi }_{^{\prime }n},\;\mathbf{\varphi  }_{^{\prime }s}\right] \left[ \begin{array}{cc} \nu _{1} & \nu _{2}\\
953
 -\nu _{2} & \nu _{1} \end{array}  \right] </math>
954
|}
955
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
956
|}
957
958
where the normal component of the gradient <math display="inline">\mathbf{\varphi }_{^{\prime }n}</math> is
959
960
{| class="formulaSCP" style="width: 100%; text-align: left;" 
961
|-
962
| 
963
{| style="text-align: left; margin:auto;width: 100%;" 
964
|-
965
| style="text-align: center;" | <math>\mathbf{\varphi }_{^{\prime }n}=\frac{\mathbf{\varphi }_{^{\prime }n}^{\left( 0\right) }}{\lambda |\mathbf{\varphi }_{^{\prime }s}|} </math>
966
|}
967
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
968
|}
969
970
In this way the contribution of the gradient at side <math display="inline">I</math> to vectors <math display="inline">\mathbf{h}_{ij}</math> results in
971
972
{| class="formulaSCP" style="width: 100%; text-align: left;" 
973
|-
974
| 
975
{| style="text-align: left; margin:auto;width: 100%;" 
976
|-
977
| style="text-align: center;" | <math>\left[ \begin{array}{c} \mathbf{h}_{11}^{T}\\
978
 \mathbf{h}_{22}^{T}\\
979
 2\mathbf{h}_{12}^{T} \end{array}  \right] ^{I}=2\left[ \begin{array}{cc} \bar{N}_{^{\prime }1}^{I} & 0\\
980
 0 & \bar{N}_{^{\prime }2}^{I}\\
981
 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left[ \begin{array}{cc} \nu _{1} & -\nu _{2}\\
982
 \nu _{2} & \nu _{1} \end{array}  \right] \left[ \begin{array}{c} \mathbf{\varphi }_{^{\prime }n}^{T}\\
983
 \mathbf{\varphi }_{^{\prime }s}^{T} \end{array}  \right] </math>
984
|}
985
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
986
|}
987
988
The only difference with Eq.([[#eq-54|54]]) necessary for the computation  of the curvature variations, is that now the contribution from the gradient at  side <math display="inline">I</math> is
989
990
{| class="formulaSCP" style="width: 100%; text-align: left;" 
991
|-
992
| 
993
{| style="text-align: left; margin:auto;width: 100%;" 
994
|-
995
| style="text-align: center;" | <math>\delta \left[ \begin{array}{c} \mathbf{h}_{11}^{T}\\
996
 \mathbf{h}_{22}^{T}\\
997
 2\mathbf{h}_{12}^{T} \end{array}  \right] ^{I}    =2\left[ \begin{array}{cc} \bar{N}_{^{\prime }1}^{I} & 0\\
998
 0 & \bar{N}_{^{\prime }2}^{I}\\
999
 \bar{N}_{^{\prime }2}^{I} & \bar{N}_{^{\prime }1}^{I} \end{array}  \right] \left[ \begin{array}{cc} \nu _{1} & -\nu _{2}\\
1000
 \nu _{2} & \nu _{1} \end{array}  \right] \left[ \begin{array}{c} \mathbf{0}\\
1001
 \frac{1}{L_{o}}\left[ \delta \mathbf{u}^{K}-\delta \mathbf{u}^{J}\right] ^{T} \end{array}  \right]</math>
1002
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
1003
|-
1004
| style="text-align: center;" | <math>   =\frac{2}{L_{0}}\left[ \begin{array}{c} -\bar{N}_{^{\prime }1}^{I}\nu _{2}\\
1005
 \bar{N}_{^{\prime }2}^{I}\nu _{1}\\
1006
 \bar{N}_{^{\prime }1}^{I}\nu _{1}-\bar{N}_{^{\prime }2}^{I}\nu _{2} \end{array}  \right] \left[ \delta \mathbf{u}^{K}-\delta \mathbf{u}^{J}\right] ^{T} </math>
1007
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
1008
|}
1009
|}
1010
1011
where the influence of variations in the length of vector <math display="inline">\mathbf{\varphi  }_{^{\prime }n}</math> has been neglected.
1012
1013
In the case of natural boundary conditions (either free or simple supported) the  gradient at the sides is supposed to be equal to the gradient in the central  element, similarly as in the membrane approach, i.e.
1014
1015
<span id="eq-71"></span>
1016
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1017
|-
1018
| 
1019
{| style="text-align: left; margin:auto;width: 100%;" 
1020
|-
1021
| style="text-align: center;" | <math>\left[ \mathbf{\varphi }_{^{\prime }1}^{I},\;\mathbf{\varphi }_{^{\prime }2} ^{I}\right] =\left[ \mathbf{\varphi }_{^{\prime }1}^{M},\;\mathbf{\varphi  }_{^{\prime }2}^{M}\right]   </math>
1022
|}
1023
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
1024
|}
1025
1026
Using Eq.([[#eq-71|71]]), vectors <math display="inline">\mathbf{h}_{ij}</math> and their variations can be easily  computed. A second possibility is to make use of the natural boundary  condition constraining the normal curvature to a zero  value. In simple supported  sides where the curvature along the side is zero, this leads to zero values for  both bending moments.
1027
1028
==6 Stiffness matrix evaluation==
1029
1030
When a predictor-corrector scheme is used to trace the movement of the shell,  the derivative (stiffness matrix) of the weak form of the momentum equations  ([[#eq-23|23]]) is needed. As usual, '''material''' and '''geometric'''  components are computed separately. The material part does not offer difficulties and  it is computed from the integral
1031
1032
<span id="eq-72"></span>
1033
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1034
|-
1035
| 
1036
{| style="text-align: left; margin:auto;width: 100%;" 
1037
|-
1038
| style="text-align: center;" | <math>\int _{A}\mathbf{B}^{T}\;\mathbf{C\;B}\;dA    </math>
1039
|}
1040
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
1041
|}
1042
1043
where matrix <math display="inline">\mathbf{B}={\boldsymbol B}_m+{\boldsymbol B}_\phi </math>  includes:
1044
1045
a) a '''membrane part''' <math display="inline">{\boldsymbol B}_m</math> computed at each mid side point <math display="inline">I</math> from
1046
1047
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1048
|-
1049
| 
1050
{| style="text-align: left; margin:auto;width: 100%;" 
1051
|-
1052
| style="text-align: center;" | <math>\delta \left[ \begin{array}{c} E_{11}\\
1053
 E_{22}\\
1054
 2E_{12} \end{array}  \right] ^{I}    =\left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }1}^{T} & \mathbf{0}_{3\times{1}}^{T}\\
1055
 \mathbf{0}_{3\times{1}}^{T} & \mathbf{\varphi }_{^{\prime }2}^{T}\\
1056
 \mathbf{\varphi }_{^{\prime }2}^{T} & \mathbf{\varphi }_{^{\prime }1}^{T} \end{array}  \right] ^{I}\left[ \begin{array}{cccc} N_{^{\prime }1}^{1} & N_{^{\prime }1}^{2} & N_{^{\prime }1}^{3} & N_{^{\prime } 1}^{I+3}\\
1057
 N_{^{\prime }2}^{1} & N_{^{\prime }2}^{2} & N_{^{\prime }2}^{3} & N_{^{\prime } 2}^{I+3} \end{array}  \right] \left[ \begin{array}{c} \delta \mathbf{u}^{1}\\
1058
 \delta \mathbf{u}^{2}\\
1059
 \delta \mathbf{u}^{3}\\
1060
 \delta \mathbf{u}^{I+3} \end{array}  \right]</math>
1061
|-
1062
| style="text-align: center;" | <math>   =\left[ \begin{array}{cc} \mathbf{\varphi }_{^{\prime }1}^{T} & \mathbf{0}_{3\times{1}}^{T}\\
1063
 \mathbf{0}_{3\times{1}}^{T} & \mathbf{\varphi }_{^{\prime }2}^{T}\\
1064
 \mathbf{\varphi }_{^{\prime }2}^{T} & \mathbf{\varphi }_{^{\prime }1}^{T} \end{array}  \right] ^{I}\left[ \begin{array}{cccc} N_{^{\prime }1}^{1\left( I\right) } & N_{^{\prime }1}^{2\left( I\right) } &  N_{^{\prime }1}^{3\left( I\right) } & N_{^{\prime }1}^{4\left( I\right) }\\
1065
 N_{^{\prime }2}^{1\left( I\right) } & N_{^{\prime }2}^{2\left( I\right) } &  N_{^{\prime }2}^{3\left( I\right) } & N_{^{\prime }2}^{4\left( I\right) } \end{array}  \right] \left[ \begin{array}{c} \delta \mathbf{u}^{1}\\
1066
 \delta \mathbf{u}^{2}\\
1067
 \delta \mathbf{u}^{3}\\
1068
 \delta \mathbf{u}^{4}={\boldsymbol B}\delta {\boldsymbol u} \end{array}  \right] ^{\left( I\right) } </math>
1069
|}
1070
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
1071
|}
1072
1073
Typically one integration point is used for computing the terms in <math display="inline">{\boldsymbol B}_m</math>.
1074
1075
b) a '''bending part''' <math display="inline">{\boldsymbol B}_\phi </math> which results from Eq.([[#eq-60|60]]) that is  constant over the element.
1076
1077
In Eq.([[#eq-72|72]]), '''C''' is the elasticity matrix integrated in the thickness. If material  non-linearity is considered, a layer wise numerical integration across the  thickness is required This leads to the tangent or algorithmic elastic-plastic  constitutive matrix <math display="inline">\mathbf{C}_{ep}</math>.
1078
1079
===6.1 Geometric stiffness (membrane part)===
1080
1081
The geometric stiffness due to membrane forces results from computing
1082
1083
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1084
|-
1085
| 
1086
{| style="text-align: left; margin:auto;width: 100%;" 
1087
|-
1088
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\;\mathbf{K}_{m}^{G}\;\Delta \mathbf{u=}\int _{A} \frac{\partial }{\partial \mathbf{u}}\left( \delta \mathbf{\varepsilon }^{T}\mathbf{N}\right) \Delta \mathbf{u}\;dA  </math>
1089
|}
1090
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
1091
|}
1092
1093
This can be written as the sum of the contributions of the three sides, i.e.
1094
1095
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1096
|-
1097
| 
1098
{| style="text-align: left; margin:auto;width: 100%;" 
1099
|-
1100
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{m}^{G}\mathbf{\;}\Delta \mathbf{u}    =\int _{A}\frac{\partial }{\partial \mathbf{u}}\left( \delta \mathbf{\varepsilon  }^{T}\right) \mathbf{N}\;\Delta \mathbf{u\;}dA</math>
1101
|-
1102
| style="text-align: center;" | <math>   =\frac{A^{M}}{3}\sum _{K=1}^{3}\sum _{I=1}^{4}\sum _{J=1}^{4}\left[ N_{^{\prime }1}^{I\left( K\right) }N_{^{\prime }1}^{J\left( K\right) }\;N_{11}^{\left( K\right) }+N_{^{\prime }2}^{I\left( K\right) } N_{^{\prime }2}^{J\left( K\right) }\;N_{22}^{\left( K\right) }\right. </math>
1103
|-
1104
| style="text-align: center;" | <math>   +\left. \left( N_{^{\prime }1}^{I\left( K\right) }N_{^{\prime } 2}^{J\left( K\right) }+N_{^{\prime }2}^{I\left( K\right) }N_{^{\prime } 1}^{J\left( K\right) }\right) N_{12}^{\left( K\right) }\right] \delta \mathbf{u}^{J\left( K\right) }\cdot \Delta \mathbf{u}^{I\left( K\right) }</math>
1105
|-
1106
| style="text-align: center;" | <math>   =\frac{A^{M}}{3}\sum _{K=1}^{3}\sum _{I=1}^{4}\sum _{J=1}^{4}\left\{ \delta \mathbf{u}^{I}\;\left[ \begin{array}{cc} N_{^{\prime }1}^{J} & N_{^{\prime }2}^{J} \end{array}  \right] \left[ \begin{array}{c}[c]{cc} N_{11} & N_{12}\\
1107
 N_{21} & N_{22} \end{array}  \right] \left[ \begin{array}{c} N_{^{\prime }1}^{J}\\
1108
 N_{^{\prime }2}^{J} \end{array}  \right] \Delta \mathbf{u}^{J}\right\} ^{\left( K\right) }</math>
1109
|-
1110
| style="text-align: center;" | <math>   .  </math>
1111
|}
1112
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
1113
|}
1114
1115
===6.2 Geometric stiffness (bending part)===
1116
1117
The geometric stiffness associated to bending moments is  more  involved. Its expression stems from
1118
1119
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1120
|-
1121
| 
1122
{| style="text-align: left; margin:auto;width: 100%;" 
1123
|-
1124
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{b}^{G}\mathbf{\;}\Delta \mathbf{u=}\int  _{A}\left( \frac{\partial }{\partial \mathbf{u}}\left( \delta \mathbf{\chi }^{T}\right) \mathbf{M}\right) \Delta \mathbf{u\;}dA  </math>
1125
|}
1126
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
1127
|}
1128
1129
Recalling expressions ([[#eq-51|51]]) and ([[#eq-52|52]]) it can be written
1130
1131
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1132
|-
1133
| 
1134
{| style="text-align: left; margin:auto;width: 100%;" 
1135
|-
1136
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{b}^{G}\Delta \mathbf{u}=A^{M}\;M_{ij} \;\Delta \left[ \delta \left( \mathbf{t}_{3}\cdot \mathbf{h}_{ij}\right) \right] </math>
1137
|}
1138
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
1139
|}
1140
1141
where
1142
1143
<span id="eq-78"></span>
1144
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1145
|-
1146
| 
1147
{| style="text-align: left; margin:auto;width: 100%;" 
1148
|-
1149
| style="text-align: center;" | <math>\Delta \left[ \delta \left( \mathbf{t}_{3}\cdot \mathbf{h}_{ij}\right) \right] =\Delta \mathbf{t}_{3}\cdot \delta \mathbf{h}_{ij}+\delta \mathbf{t}_{3}\cdot \Delta \mathbf{h}_{ij}+\Delta \left( \delta \mathbf{t}_{3}\right) \cdot \mathbf{h}_{ij}  </math>
1150
|}
1151
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
1152
|}
1153
1154
The first two terms lead to  symmetric components. The  second one <math display="inline">\left( \delta \mathbf{t}_{3}\cdot \Delta \mathbf{h}_{ij}\right) </math>  can be expressed as
1155
1156
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1157
|-
1158
| 
1159
{| style="text-align: left; margin:auto;width: 100%;" 
1160
|-
1161
| style="text-align: center;" | <math>=\left\{ -\sum _{J=1}^{3}\left[ \bar{N}_{^{\prime }1}^{J}\mathbf{\tilde  {\varphi }}_{^{\prime }1}+\bar{N}_{^{\prime }2}^{J}\mathbf{\tilde{\varphi }}_{^{\prime }2}\right] \left( \mathbf{t}_{3}\cdot \delta \mathbf{u}^{J}\right) \right\} \cdot \left\{ \sum _{I=1}^{3}\sum _{K=1}^{4}\left( \bar{N}_{^{\prime  }i}^{I}\;N_{^{\prime }j}^{K}+\bar{N}_{^{\prime }j}^{I}\;N_{i}^{K}\right) \Delta \mathbf{u}^{K\left( I\right) }\right\}</math>
1162
|-
1163
| style="text-align: center;" | <math>   =-\sum _{J=1}^{3}\sum _{I=1}^{3}\sum _{K=1}^{4}\left( \delta \mathbf{u}^{J}\right) ^{T}\left[ \bar{N}_{^{\prime }1}^{J}\left( \mathbf{t}_{3} \otimes \mathbf{\tilde{\varphi }}_{^{\prime }1}\right) +\bar{N}_{^{\prime }2} ^{J}\left( \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }2}\right) \right] \left( \bar{N}_{^{\prime }i}^{I}\;N_{^{\prime }j}^{K}+\bar  {N}_{^{\prime }j}^{I}\;N_{i}^{K}\right) \Delta \mathbf{u}^{K\left( I\right) }</math>
1164
|-
1165
| style="text-align: center;" | <math>   \;\;
1166
 </math>
1167
|}
1168
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
1169
|}
1170
1171
Then, a first component of the geometric bending stiffness  can be  written as
1172
1173
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1174
|-
1175
| 
1176
{| style="text-align: left; margin:auto;width: 100%;" 
1177
|-
1178
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{b1}^{G}\Delta \mathbf{u}    =A^{M}\left\{ -\sum _{J=1}^{3}\left( \delta \mathbf{u}^{J}\right) ^{T}\left[ \begin{array}{c}[c]{cc} \bar{N}_{^{\prime }1}^{J} \mathbf{1} & \bar{N}_{^{\prime }2}^{J} \mathbf{1}\end{array}  \right] \left[ \begin{array}{c} \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }1}\\
1179
 \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }2} \end{array}  \right] \right\}</math>
1180
|-
1181
| style="text-align: center;" | <math>   \sum _{I=1}^{3}\left[ \begin{array}{cc} \bar{N}_{^{\prime }1}^{I} & \bar{N}_{^{\prime }2}^{I} \end{array}  \right] \left[ \begin{array}{cc} M_{11} & M_{12}\\
1182
 M_{12} & M_{22} \end{array}  \right] \sum _{K=1}^{4}\left[ \begin{array}{c} N_{^{\prime }1}^{K\left( I\right) }\\
1183
 N_{^{\prime }2}^{K\left( I\right) } \end{array}  \right] \Delta \mathbf{u}^{K\left( I\right) }</math>
1184
|-
1185
| style="text-align: center;" | 
1186
|}
1187
| style="width: 5px;text-align: right;white-space: nowrap;" | (80)
1188
|}
1189
1190
The last term in ([[#eq-78|78]]), can be obtained observing that (sum in <math display="inline">\alpha </math> and <math display="inline">\beta </math>)
1191
1192
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1193
|-
1194
| 
1195
{| style="text-align: left; margin:auto;width: 100%;" 
1196
|-
1197
| style="text-align: center;" | <math>\Delta \left( \delta \mathbf{t}_{3}\right)    =\left[ \Delta \left( \delta \mathbf{t}_{3}\right) \right] _{1}\mathbf{\tilde{\varphi }}_{^{\prime  }1}+\left[ \Delta \left( \delta \mathbf{t}_{3}\right) \right] _{2} \mathbf{\tilde{\varphi }}_{^{\prime }2}+\left[ \Delta \left( \delta  \mathbf{t}_{3}\right) \right] _{3}\mathbf{t}_{3}</math>
1198
|-
1199
| style="text-align: center;" | <math>   =-\left( \mathbf{\tilde{\varphi }}_{^{\prime }\alpha }\cdot \mathbf{\tilde  {\varphi }}_{^{\prime }\beta }\right) \left[ \delta \mathbf{\varphi }_{^{\prime  }\alpha }^{T}\left( \mathbf{t}_{3}\otimes \mathbf{t}_{3}\right) \Delta  \mathbf{\varphi }_{^{\prime }\beta }\right] \mathbf{t}_{3}+</math>
1200
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
1201
|-
1202
| style="text-align: center;" | <math>   +\left[ \delta \mathbf{\varphi }_{^{\prime }\alpha }^{T}\left( \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }\alpha }\right) \Delta  \mathbf{\varphi }_{^{\prime }\beta }\right] \mathbf{\tilde{\varphi }}_{^{\prime  }\beta }+\left[ \delta \mathbf{\varphi }_{^{\prime }\alpha }^{T}\left( \mathbf{\tilde{\varphi }}_{^{\prime }\beta }\otimes \mathbf{t}_{3}\right) \Delta \mathbf{\varphi }_{^{\prime }\beta }\right] \mathbf{\tilde{\varphi }}_{^{\prime }\alpha } </math>
1203
|}
1204
|}
1205
1206
then
1207
1208
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1209
|-
1210
| 
1211
{| style="text-align: left; margin:auto;width: 100%;" 
1212
|-
1213
| style="text-align: center;" | <math>\Delta \left( \delta \mathbf{t}_{3}\right) \cdot \mathbf{h}_{ij}    =-\left( \mathbf{\tilde{\varphi }}_{^{\prime }\alpha }\cdot \mathbf{\tilde{\varphi }}_{^{\prime }\beta }\right) \left[ \delta \mathbf{\varphi }_{^{\prime }\alpha  }^{T}\left( \mathbf{t}_{3}\otimes \mathbf{t}_{3}\right) \Delta \mathbf{\varphi  }_{^{\prime }\beta }\right] \left( \mathbf{t}_{3}\cdot \mathbf{h}_{ij}\right) +</math>
1214
|-
1215
| style="text-align: center;" | <math>   +\left[ \delta \mathbf{\varphi }_{^{\prime }\alpha }^{T}\left( \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }\alpha }\right) \Delta  \mathbf{\varphi }_{^{\prime }\beta }\right] \left( \mathbf{\tilde{\varphi }}_{^{\prime }\beta }\cdot \mathbf{h}_{ij}\right)</math>
1216
|-
1217
| style="text-align: center;" | <math>   +\left[ \delta \mathbf{\varphi }_{^{\prime }\alpha }^{T}\left( \mathbf{\tilde  {\varphi }}_{^{\prime }\beta }\otimes \mathbf{t}_{3}\right) \Delta \mathbf{\varphi  }_{^{\prime }\beta }\right] \left( \mathbf{\tilde{\varphi }}_{^{\prime }\alpha  }\cdot \mathbf{h}_{ij}\right)</math>
1218
|-
1219
| style="text-align: center;" | <math>   =\delta \mathbf{\varphi }_{^{\prime }\alpha }^{T}\left\{ \rho _{ij}^{\beta  }\left( \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }\alpha  }\right) +\rho _{ij}^{\alpha }\left( \mathbf{\tilde{\varphi }}_{^{\prime }\beta  }\otimes \mathbf{t}_{3}\right) -g^{\alpha \beta }\kappa _{ij}\left( \mathbf{t}_{3}\otimes \mathbf{t}_{3}\right) \right\} \Delta \mathbf{\varphi  }_{^{\prime }\beta } </math>
1220
|}
1221
| style="width: 5px;text-align: right;white-space: nowrap;" | (82)
1222
|}
1223
1224
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1225
|-
1226
| 
1227
{| style="text-align: left; margin:auto;width: 100%;" 
1228
|-
1229
| style="text-align: center;" | <math>\Delta \left( \delta \mathbf{t}_{3}\right) \cdot \mathbf{h}_{ij}=\sum _{J=1} ^{3}\sum _{K=1}^{3}\sum _{\alpha=1}^{2}\sum _{\beta=1}^{2}\bar{N}_{^{\prime  }\alpha }^{J}\bar{N}_{^{\prime }\beta }^{K}\left[ \delta \mathbf{u}^{J}\right] ^{T}\;\mathbf{P}_{ij}^{a\beta }\mathbf{\;}\Delta \mathbf{u}^{K} </math>
1230
|}
1231
| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
1232
|}
1233
1234
with
1235
1236
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1237
|-
1238
| 
1239
{| style="text-align: left; margin:auto;width: 100%;" 
1240
|-
1241
| style="text-align: center;" | <math>\mathbf{P}_{ij}^{\alpha \beta }=\left[ -a^{\alpha \beta }\;\kappa _{ij}\;\left( \mathbf{t}_{3}\otimes \mathbf{t}_{3}\right) +\rho _{ij}^{\beta }\left( \mathbf{t}_{3}\otimes \mathbf{\tilde{\varphi }}_{^{\prime }\alpha }\right) +\rho _{ij}^{\alpha }\left( \mathbf{\tilde{\varphi }}_{^{\prime }\beta } \otimes \mathbf{t}_{3}\right) \right] </math>
1242
|}
1243
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
1244
|}
1245
1246
where the covariant metric tensor at the middle surface <math display="inline">a^{\alpha \beta }</math> has  been used.
1247
1248
Denoting the sum
1249
1250
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1251
|-
1252
| 
1253
{| style="text-align: left; margin:auto;width: 100%;" 
1254
|-
1255
| style="text-align: center;" | <math>\mathbf{Q}^{\alpha \beta }=A^{M}\;\sum _{i=1}^{2}\sum _{j=1}^{2}\mathbf{P}_{ij}^{a\beta }M_{\iota j} </math>
1256
|}
1257
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
1258
|}
1259
1260
the term
1261
1262
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1263
|-
1264
| 
1265
{| style="text-align: left; margin:auto;width: 100%;" 
1266
|-
1267
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{b2}^{G}\Delta \mathbf{u}=A^{M}\;M_{ij}\left( \Delta \left( \delta \mathbf{t}_{3}\right) \cdot \mathbf{h}_{ij}\right) </math>
1268
|}
1269
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
1270
|}
1271
1272
results in
1273
1274
<span id="eq-87"></span>
1275
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1276
|-
1277
| 
1278
{| style="text-align: left; margin:auto;width: 100%;" 
1279
|-
1280
| style="text-align: center;" | <math>\delta \mathbf{u}^{T}\mathbf{K}_{b2}^{G}\Delta \mathbf{u}=\sum _{J=1}^{3} \sum _{K=1}^{3}\sum _{\alpha=1}^{2}\sum _{\beta=1}^{2}\bar{N}_{^{\prime }\alpha  }^{J}\bar{N}_{^{\prime }\beta }^{K}\left[ \delta \mathbf{u}^{J}\right] ^{T}\;\mathbf{Q}^{a\beta }\mathbf{\;}\Delta \mathbf{u}^{K}  </math>
1281
|}
1282
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
1283
|}
1284
1285
The last expression ([[#eq-87|87]]) has components only in the nodes of the  central element, which stems from the definition used for '''t'''<math display="inline">_{3}</math>.
1286
1287
Numerical experiments have shown that the bending part of the geometric stiffness is not so important and can be disregarded in the iterative process.
1288
1289
==7 Numerical examples in the linear range==
1290
1291
In this section, a summary of the results obtained in the assessment of the  proposed element in the linear range are presented. In this first part of the  numerical experiments, results have been obtained with a static/dynamic code  that uses a implicit solution of the discretized equilibrium equations. In the examples,  the original BST element <span id='citeF-10'></span>[[#cite-10|[10]]] is compared with the enhanced version here proposed, denoted by CBST when a numerical integration is  performed with three points, and  CBST1 when only one integration point is  used. Note that one integration quadrature is equivalent to averaging  the metric tensors computed  at each side.
1292
1293
===7.1 Membrane patch test===
1294
1295
The main objective of the present formulation is to obtain a (non conforming)  membrane element  with a similar performance than the linear strain triangle,  that satisfies the patch test. To assess this, a square domain subjected to  nodal forces associated to a constant stress state (in both directions and in  shear) is considered. Figure [[#img-3|3]] shows the patch of elements and the  necessary nodal forces to obtain a uniform tensile stress in direction <math display="inline">x_{1} </math>. Note that the nodal forces used correspond to the constant strain  triangle, not to the linear one. For the distorted mesh shown, correct results  are obtained using either  1 or 3 integration points.
1296
1297
<div id='img-3'></div>
1298
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1299
|-
1300
|[[Image:Draft_Samper_340530075-test-fig3.png|600px|Membrane patch  test]]
1301
|- style="text-align: center; font-size: 75%;"
1302
| colspan="1" | '''Figure 3:''' Membrane patch  test
1303
|}
1304
1305
===7.2 Bending patch test (torsion)===
1306
1307
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-4|4]] 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 displaced from the center, the results  obtained with the CBST element are not correct. The original element BST gives correct  results in all cases, if natural  boundary conditions are imposed in the formulation. If this is  not the case, incorrent  results are obtained even with structured meshes.
1308
1309
<div id='img-4'></div>
1310
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1311
|-
1312
|[[Image:Draft_Samper_340530075-test-fig4.png|600px|Patch test for uniform  torsion]]
1313
|- style="text-align: center; font-size: 75%;"
1314
| colspan="1" | '''Figure 4:''' Patch test for uniform  torsion
1315
|}
1316
1317
===7.3 Cook's membrane problem===
1318
1319
This example is used to assess the membrane performance of the CBST element  and to compare it with the 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-5|5]].a shows the geometry and the applied load. Figure [[#img-5|5]].b  plots the vertical displacement of the upper vertex as a function of the  number of nodes in the mesh. In this figure results obtained with other  isoparametric elements have been also plotted for comparison. They include  the constain strain triangle (CST), the bilinear  quadrilateral (QUAD4) and the linear strain triangle (LST). Note that in this membrane problem both the BST and the CST elements give identical results.
1320
1321
<div id='img-5'></div>
1322
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1323
|-
1324
|[[Image:Draft_Samper_340530075-test-fig5.png|600px|Cook membrane problem (a)  Geometry (b) Results]]
1325
|- style="text-align: center; font-size: 75%;"
1326
| colspan="1" | '''Figure 5:''' Cook membrane problem (a)  Geometry (b) Results
1327
|}
1328
1329
From the plot shown it can be seen that the new element with three integration  points (CBST)  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 (CBST1) the convergence in the reported displacement is notably  better that for the rest of the elements.
1330
1331
===7.4 Cylindrical roof===
1332
1333
In this example an effective membrane interpolation is of primary importance. Hence this is good test to assess the new element. The geometry is a  cylindrical roof supported by a rigid diaphragm at both ends and it is loaded by  a uniform dead weight (see [[#img-6|6]].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-6|6]] shows orientation B).
1334
1335
Tables [[#table-1|1]], [[#table-2|2]] and [[#table-3|3]] present the normalized  vertical displacements at the crown (point A) and in the midpoint of the free  side (point B) for the two orientations of structured meshes and  non-structured meshes. 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 <span id='citeF-13'></span>[[#cite-13|[13]]].
1336
1337
<div id='img-6'></div>
1338
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1339
|-
1340
|[[Image:Draft_Samper_340530075-test-fig6.png|600px|Cylindrical roof under dead  weight. E=3 E⁶, ν=0.0, Thickness =3.0, shell weight =0.625 per unit  area.]]
1341
|- style="text-align: center; font-size: 75%;"
1342
| colspan="1" | '''Figure 6:''' Cylindrical roof under dead  weight. E=3 E<math>^{6}</math>, <math>\nu=0.0</math>, Thickness =3.0, shell weight =0.625 per unit  area.
1343
|}
1344
1345
1346
{|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
1347
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Cylindrical roof under dead weigth. Normalized displacements for mesh  orientation A
1348
|- style="border-top: 2px solid;"
1349
| style="border-left: 2px solid;border-right: 2px solid;" |  
1350
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Point-A
1351
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Point-B
1352
|- style="border-top: 2px solid;"
1353
| style="border-left: 2px solid;border-right: 2px solid;" |  NDOFs 
1354
| style="border-left: 2px solid;border-right: 2px solid;" | CBST 
1355
| style="border-left: 2px solid;border-right: 2px solid;" | CBST1 
1356
| style="border-left: 2px solid;border-right: 2px solid;" | BST 
1357
| style="border-left: 2px solid;border-right: 2px solid;" | CBST 
1358
| style="border-left: 2px solid;border-right: 2px solid;" | CBST1 
1359
| style="border-left: 2px solid;border-right: 2px solid;" | BST
1360
|-
1361
| style="border-left: 2px solid;border-right: 2px solid;" | 16 
1362
| style="border-left: 2px solid;border-right: 2px solid;" | 0.65724 
1363
| style="border-left: 2px solid;border-right: 2px solid;" | 0.91855 
1364
| style="border-left: 2px solid;border-right: 2px solid;" | 0.74161 
1365
| style="border-left: 2px solid;border-right: 2px solid;" | 0.40950 
1366
| style="border-left: 2px solid;border-right: 2px solid;" | 0.70100 
1367
| style="border-left: 2px solid;border-right: 2px solid;" | 1.35230
1368
|-
1369
| style="border-left: 2px solid;border-right: 2px solid;" | 56 
1370
| style="border-left: 2px solid;border-right: 2px solid;" | 0.53790 
1371
| style="border-left: 2px solid;border-right: 2px solid;" | 1.03331 
1372
| style="border-left: 2px solid;border-right: 2px solid;" | 0.74006 
1373
| style="border-left: 2px solid;border-right: 2px solid;" | 0.54859 
1374
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00759 
1375
| style="border-left: 2px solid;border-right: 2px solid;" | 0.75590
1376
|-
1377
| style="border-left: 2px solid;border-right: 2px solid;" | 208 
1378
| style="border-left: 2px solid;border-right: 2px solid;" | 0.89588 
1379
| style="border-left: 2px solid;border-right: 2px solid;" | 1.04374 
1380
| style="border-left: 2px solid;border-right: 2px solid;" | 0.88491 
1381
| style="border-left: 2px solid;border-right: 2px solid;" | 0.91612 
1382
| style="border-left: 2px solid;border-right: 2px solid;" | 1.02155 
1383
| style="border-left: 2px solid;border-right: 2px solid;" | 0.88269
1384
|-
1385
| style="border-left: 2px solid;border-right: 2px solid;" | 800 
1386
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99658 
1387
| style="border-left: 2px solid;border-right: 2px solid;" | 1.01391 
1388
| style="border-left: 2px solid;border-right: 2px solid;" | 0.96521 
1389
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99263 
1390
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00607 
1391
| style="border-left: 2px solid;border-right: 2px solid;" | 0.96393
1392
|- style="border-bottom: 2px solid;"
1393
| style="border-left: 2px solid;border-right: 2px solid;" | 3136 
1394
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00142 
1395
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00385 
1396
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99105 
1397
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99881 
1398
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00102 
1399
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98992
1400
1401
|}
1402
1403
1404
{|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
1405
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Cylindrical roof under dead weigth. Normalized displacements for mesh  orientation A
1406
|- style="border-top: 2px solid;"
1407
| style="border-left: 2px solid;border-right: 2px solid;" |  
1408
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Point-A
1409
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Point-B
1410
|- style="border-top: 2px solid;"
1411
| style="border-left: 2px solid;border-right: 2px solid;" |  NDOFs 
1412
| style="border-left: 2px solid;border-right: 2px solid;" | CBST 
1413
| style="border-left: 2px solid;border-right: 2px solid;" | CBST1 
1414
| style="border-left: 2px solid;border-right: 2px solid;" | BST 
1415
| style="border-left: 2px solid;border-right: 2px solid;" | CBST 
1416
| style="border-left: 2px solid;border-right: 2px solid;" | CBST1 
1417
| style="border-left: 2px solid;border-right: 2px solid;" | BST
1418
|-
1419
| style="border-left: 2px solid;border-right: 2px solid;" | 16 
1420
| style="border-left: 2px solid;border-right: 2px solid;" | 0.26029 
1421
| style="border-left: 2px solid;border-right: 2px solid;" | 0.83917 
1422
| style="border-left: 2px solid;border-right: 2px solid;" | 0.40416 
1423
| style="border-left: 2px solid;border-right: 2px solid;" | 0.52601 
1424
| style="border-left: 2px solid;border-right: 2px solid;" | 0.86133 
1425
| style="border-left: 2px solid;border-right: 2px solid;" | 0.89778
1426
|-
1427
| style="border-left: 2px solid;border-right: 2px solid;" | 56 
1428
| style="border-left: 2px solid;border-right: 2px solid;" | 0.81274 
1429
| style="border-left: 2px solid;border-right: 2px solid;" | 1.10368 
1430
| style="border-left: 2px solid;border-right: 2px solid;" | 0.61642 
1431
| style="border-left: 2px solid;border-right: 2px solid;" | 0.67898 
1432
| style="border-left: 2px solid;border-right: 2px solid;" | 0.93931 
1433
| style="border-left: 2px solid;border-right: 2px solid;" | 0.68238
1434
|-
1435
| style="border-left: 2px solid;border-right: 2px solid;" | 208 
1436
| style="border-left: 2px solid;border-right: 2px solid;" | 0.97651 
1437
| style="border-left: 2px solid;border-right: 2px solid;" | 1.04256 
1438
| style="border-left: 2px solid;border-right: 2px solid;" | 0.85010 
1439
| style="border-left: 2px solid;border-right: 2px solid;" | 0.93704 
1440
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00255 
1441
| style="border-left: 2px solid;border-right: 2px solid;" | 0.86366
1442
|-
1443
| style="border-left: 2px solid;border-right: 2px solid;" | 800 
1444
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00085 
1445
| style="border-left: 2px solid;border-right: 2px solid;" | 1.01195 
1446
| style="border-left: 2px solid;border-right: 2px solid;" | 0.95626 
1447
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99194 
1448
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00211 
1449
| style="border-left: 2px solid;border-right: 2px solid;" | 0.95864
1450
|- style="border-bottom: 2px solid;"
1451
| style="border-left: 2px solid;border-right: 2px solid;" | 3136 
1452
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00129 
1453
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00337 
1454
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98879 
1455
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99828 
1456
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00017 
1457
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98848
1458
1459
|}
1460
1461
1462
{|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
1463
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Cylindrical roof under dead weigth. Normalized displacements for  non-structured mesh
1464
|- style="border-top: 2px solid;"
1465
| style="border-left: 2px solid;border-right: 2px solid;" |  
1466
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Point-A
1467
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Point-B
1468
|- style="border-top: 2px solid;"
1469
| style="border-left: 2px solid;border-right: 2px solid;" |  NDOFs 
1470
| style="border-left: 2px solid;border-right: 2px solid;" | CBST 
1471
| style="border-left: 2px solid;border-right: 2px solid;" | CBST1 
1472
| style="border-left: 2px solid;border-right: 2px solid;" | BST 
1473
| style="border-left: 2px solid;border-right: 2px solid;" | CBST 
1474
| style="border-left: 2px solid;border-right: 2px solid;" | CBST1 
1475
| style="border-left: 2px solid;border-right: 2px solid;" | BST
1476
|-
1477
| style="border-left: 2px solid;border-right: 2px solid;" | 851 
1478
| style="border-left: 2px solid;border-right: 2px solid;" | 0.97546 
1479
| style="border-left: 2px solid;border-right: 2px solid;" | 0.8581 
1480
| style="border-left: 2px solid;border-right: 2px solid;" | 0.97598 
1481
| style="border-left: 2px solid;border-right: 2px solid;" | 0.97662 
1482
| style="border-left: 2px solid;border-right: 2px solid;" | 1.0027 
1483
| style="border-left: 2px solid;border-right: 2px solid;" | 0.97194
1484
|-
1485
| style="border-left: 2px solid;border-right: 2px solid;" | 3311 
1486
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98729 
1487
| style="border-left: 2px solid;border-right: 2px solid;" | 0.9682 
1488
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98968 
1489
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98476 
1490
| style="border-left: 2px solid;border-right: 2px solid;" | 1.0083 
1491
| style="border-left: 2px solid;border-right: 2px solid;" | 0.98598
1492
|- style="border-bottom: 2px solid;"
1493
| style="border-left: 2px solid;border-right: 2px solid;" | 13536 
1494
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99582 
1495
| style="border-left: 2px solid;border-right: 2px solid;" | 0.9992 
1496
| style="border-left: 2px solid;border-right: 2px solid;" | 1.00057 
1497
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99316 
1498
| style="border-left: 2px solid;border-right: 2px solid;" | 0.9973 
1499
| style="border-left: 2px solid;border-right: 2px solid;" | 0.99596
1500
1501
|}
1502
1503
Figure [[#img-6|6]].b shows the normalized displacement of point-B  over the structured meshes as a function of the number of degrees of freedom. An  excelent convergence the CBST element can be seen. The version with only one  integration point (CBST1) presents a behavior a little more flexible and  converges from above. Table [[#table-3|3]] shows that for non-structured meshes  the element converges to the correct value but more slowly.
1504
1505
===7.5 Open semi-espherical dome with point loads===
1506
1507
The main problem of finite elements with initially curved geometry is the so  called membrane locking. The CBST element  proposed  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-7|7]].a shows the discretized geometry (only one quarter of the geometry  is considered due to symmetry).
1508
1509
<div id='img-7'></div>
1510
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1511
|-
1512
|[[Image:Draft_Samper_340530075-test-fig7.png|600px|Pinched hemispherical shell with a hole, (a)geometry, (b)normalized  displacement]]
1513
|- style="text-align: center; font-size: 75%;"
1514
| colspan="1" | '''Figure 7:''' Pinched hemispherical shell with a hole, (a)geometry, (b)normalized  displacement
1515
|}
1516
1517
In Figure [[#img-7|7]].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 and a transverse  shear-deformable quadrilateral (QUAD); a transverse shear deformable quadratic  triangle (TRIA) <span id='citeF-15'></span>[[#cite-1|[1]]] that is vulnerable to locking and an  assumed strain version of this triangle (TRIA-1)<sup><span id='citeF-1'></span>[[#cite-1|[1]]]</sup> that  alleviates membrane locking but does not eliminate it.
1518
1519
From the plotted results it can be seen that the  CBST element presents  membrane locking in bending dominated problems with initially curved  geometries. This locking is much less severe than in a standard quadratic  triangle and even than for the assumed strain version used for comparison. When only  one integration point is used (CBST1 element) membrane locking disappears.
1520
1521
==8 Non linear numerical examples==
1522
1523
Results for examples with  geometric and material non-linearities are presented next using the CBST1 element. Due to the features of  the modelized problems, with strong nonlinearities associated to instabilities  and frictional contact conditions, a code with explicit integration of the dynamic  equilibrium equations has been used <span id='citeF-14'></span>[[#cite-14|[14]]]. This code allows to  obtain pseudo-static solutions throught dynamic relaxation.
1524
1525
===8.1 Inflation of a sphere===
1526
1527
The example is the inflation of a spherical shell under internal pressure. An  incompressible Mooney-Rivlin constitutive material have 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 <span id='citeF-15'></span>[[#cite-15|[15]]] (with <math display="inline">\gamma =R/R^{\left( 0\right) }</math>):
1528
1529
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1530
|-
1531
| 
1532
{| style="text-align: left; margin:auto;width: 100%;" 
1533
|-
1534
| style="text-align: center;" | <math>  p=\frac{h^{\left( 0\right) }}{R^{\left( 0\right) }\gamma ^{2}}\frac{dW}{d\gamma }=\frac{8h^{\left( 0\right) }}{R^{\left( 0\right) }\gamma ^{2} }\left( \gamma ^{6}-1\right) \left( \mu _{1}-\mu _{2}\gamma ^{2}\right) </math>
1535
|}
1536
|}
1537
1538
In this numerical simulation the same geometric and material parameters used  in Ref. <span id='citeF-16'></span>[[#cite-16|[16]]] have been adopted: <math display="inline">R^{\left( 0\right) }=1</math> and <math display="inline">h^{\left( 0\right) }=0.02</math>. The three meshes considered to evaluate  convergence are shown in Figure [[#img-8|8]].a. The value of the actual radius   as a function of the internal pressure is plotted in Figure  [[#img-8|8]].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 thickness radius ratio of <math display="inline">h/R=0.00024</math>.
1539
1540
<div id='img-8'></div>
1541
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1542
|-
1543
|[[Image:Draft_Samper_340530075-test-fig8a.png|600px|]]
1544
|[[Image:Draft_Samper_340530075-test-fig8b.png|300px|Inflation  of sphere of Mooney-Rivlin material. (a) Meshes used in the analysis (b) Change of radius as a function of the internal pressure.]]
1545
|- style="text-align: center; font-size: 75%;"
1546
| colspan="2" | '''Figure 8:''' Inflation  of sphere of Mooney-Rivlin material. (a) Meshes used in the analysis (b) Change of radius as a function of the internal pressure.
1547
|}
1548
1549
===8.2 Inflation of an air-bag===
1550
1551
This example has also been taken from Ref.<span id='citeF-16'></span>[[#cite-16|[16]]] 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 <span id='citeF-16'></span>[[#cite-16|[16]]] it is also  discussed the important regularizing properties of the bending energy, that when  disregarded leads to massive wrinkling in the compressed zones.
1552
1553
The air bag geometry is initially square with an undeformed diagonal of <math display="inline">1.20</math>. The constitutive material is a linear isotropic elastic one with  modulus of elasticity <math display="inline">E=5.88\times{10}^{8}</math> and Poisson's ratio <math display="inline">\nu=0.4</math>. Only  one quarter of the geometry has been modelled due to symmetry. Only the normal  to the original plane is constrained along the boundaries. The thickness  considered is <math display="inline">h=0.001</math> and the inflation pressure is <math display="inline">5000</math>. Using a density <math display="inline">\delta=1000</math>, pressure is linearly incremented from 0 to the final value in <math display="inline">t=0.1</math>.
1554
1555
With comparative purposes and also to backup the comments in Ref. <span id='citeF-16'></span>[[#cite-16|[16]]]  two analyses have been performed, a purely membranal one and one including bending  effects. Figure [[#img-9|9]] shows the final deformed configurations for three  meshes with 289, 1089 and 4225 nodes. The top row corresponds to a full  analysis including bending and the central row is a pure  membrane  analysis.  The bottom row is also an analysis including bending where the mesh  orientation has been changed.
1556
1557
<div id='img-9'></div>
1558
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1559
|-
1560
|[[Image:Draft_Samper_340530075-test-fig9.png|500px|Inflation of a square  air-bag . Deformed configurations for three different meshes with 800, 3136  and 12416 degrees of freedom.]]
1561
|- style="text-align: center; font-size: 75%;"
1562
| colspan="1" | '''Figure 9:''' Inflation of a square  air-bag . Deformed configurations for three different meshes with 800, 3136  and 12416 degrees of freedom.
1563
|}
1564
1565
The top and bottom lines  show the final shapes change according  to the degree of discretization and mesh orientation due to instabilities and  non uniqueness of the solution. The central row shows the noteworthy increment  of wrinkles as the mesh is refined. Finally in Figure [[#img-10|10]] presents  the convergence of the maximum displacement of the central point as the mesh  is refined.
1566
1567
<div id='img-10'></div>
1568
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1569
|-
1570
|[[Image:Draft_Samper_340530075-test-fig10.png|500px|Inflation of a square  air-bag. Convergence of the maximum displacement versus the number of degrees  of freedom in the model.]]
1571
|- style="text-align: center; font-size: 75%;"
1572
| colspan="1" | '''Figure 10:''' Inflation of a square  air-bag. Convergence of the maximum displacement versus the number of degrees  of freedom in the model.
1573
|}
1574
1575
===8.3 Deep drawing of a rolled sheet===
1576
1577
The element performance in problems involving large strains and  anisotropic plastic behaviour is assesed in a benchmark proposed in a recent  NUMISHEET <span id='citeF-175'></span>[[#cite-17|[17]]] meeting. This is the deep  drawing of a circular mild steel sheet with a spherical punch of radius <math display="inline">50 mm</math> (<math display="inline">R_{p}</math>). The original radius of the sheet (<math display="inline">R_{s}</math>) is <math display="inline">100 mm</math>  (<math display="inline">\frac{R_{L}}{R_{P}}=\beta=2</math>). The drawing depth is (<math display="inline">DD</math>) is 85 mm  (<math display="inline">\frac{DD}{R_{P}}=1.7</math>) and the force over the blank holder is <math display="inline">80 kN</math>. The  constitutive material of the rolled sheet is characterized by three tensile  tests along three different directions respect to rolling direction (<math display="inline">RD</math>) (Table  [[#table-4|4]]). Note that the modellization requires to treat  the material as transversally anisotropic, as this improves the deep  drawability and if not considered leads to an early localized necking of the  material. The yield function originally proposed by Hill<sup><span id='citeF-11'></span>[[#cite-11|[11]]]</sup> has been chosen for both associative and non-associative  potential functions. For elastic-plastic problems a numerical integration along  the thickness direction is necessary. This is accomplished here using four integration points.
1578
1579
1580
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1581
|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 Plastic characterization of the constitutive material
1582
|- style="border-top: 2px solid;"
1583
| style="border-left: 2px solid;border-right: 2px solid;" |  Thickness 
1584
| style="border-left: 2px solid;border-right: 2px solid;" | Orientation 
1585
| style="border-left: 2px solid;border-right: 2px solid;" | Yield 
1586
| style="border-left: 2px solid;border-right: 2px solid;" | Tensile 
1587
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\varepsilon _{u}</math>
1588
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\varepsilon  _{t}</math>
1589
| style="border-left: 2px solid;border-right: 2px solid;" | n 
1590
| style="border-left: 2px solid;border-right: 2px solid;" | r
1591
|-
1592
| style="border-left: 2px solid;border-right: 2px solid;" | 
1593
| style="border-left: 2px solid;border-right: 2px solid;" | respect to RD 
1594
| style="border-left: 2px solid;border-right: 2px solid;" | stres 
1595
| style="border-left: 2px solid;border-right: 2px solid;" | strength 
1596
| style="border-left: 2px solid;border-right: 2px solid;" | (uniform) 
1597
| style="border-left: 2px solid;border-right: 2px solid;" | (total) 
1598
| style="border-left: 2px solid;border-right: 2px solid;" | 
1599
| style="border-left: 2px solid;border-right: 2px solid;" | 
1600
|- style="border-top: 2px solid;"
1601
| style="border-left: 2px solid;border-right: 2px solid;" |  [mm] 
1602
| style="border-left: 2px solid;border-right: 2px solid;" | [<math display="inline">^{o}</math>] 
1603
| style="border-left: 2px solid;border-right: 2px solid;" | [N/mm<math display="inline">^{2}</math>] 
1604
| style="border-left: 2px solid;border-right: 2px solid;" | [N/mm<math display="inline">^{2}</math>] 
1605
| style="border-left: 2px solid;border-right: 2px solid;" | [%] 
1606
| style="border-left: 2px solid;border-right: 2px solid;" | [%] 
1607
| style="border-left: 2px solid;border-right: 2px solid;" | - 
1608
| style="border-left: 2px solid;border-right: 2px solid;" | -
1609
|- style="border-top: 2px solid;"
1610
| style="border-left: 2px solid;border-right: 2px solid;" |  
1611
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
1612
| style="border-left: 2px solid;border-right: 2px solid;" | 176 
1613
| style="border-left: 2px solid;border-right: 2px solid;" | 322 
1614
| style="border-left: 2px solid;border-right: 2px solid;" | 24 
1615
| style="border-left: 2px solid;border-right: 2px solid;" | 40 
1616
| style="border-left: 2px solid;border-right: 2px solid;" | 0.214 
1617
| style="border-left: 2px solid;border-right: 2px solid;" | 1.73
1618
|- style="border-top: 2px solid;"
1619
| style="border-left: 2px solid;border-right: 2px solid;" |  0.98 
1620
| style="border-left: 2px solid;border-right: 2px solid;" | 45 
1621
| style="border-left: 2px solid;border-right: 2px solid;" | 185 
1622
| style="border-left: 2px solid;border-right: 2px solid;" | 333 
1623
| style="border-left: 2px solid;border-right: 2px solid;" | 22 
1624
| style="border-left: 2px solid;border-right: 2px solid;" | 39 
1625
| style="border-left: 2px solid;border-right: 2px solid;" | 0.203 
1626
| style="border-left: 2px solid;border-right: 2px solid;" | 1.23
1627
|- style="border-top: 2px solid;border-bottom: 2px solid;"
1628
| style="border-left: 2px solid;border-right: 2px solid;" |  
1629
| style="border-left: 2px solid;border-right: 2px solid;" | 90 
1630
| style="border-left: 2px solid;border-right: 2px solid;" | 180 
1631
| style="border-left: 2px solid;border-right: 2px solid;" | 319 
1632
| style="border-left: 2px solid;border-right: 2px solid;" | 23 
1633
| style="border-left: 2px solid;border-right: 2px solid;" | 44 
1634
| style="border-left: 2px solid;border-right: 2px solid;" | 0.206 
1635
| style="border-left: 2px solid;border-right: 2px solid;" | 2.02
1636
1637
|}
1638
1639
Only one quarter of the geometry has been considered due to symmetry. Tools  are modelized as rigid surfaces, the punch has been defined by 1439  points and 2730 triangles. For the die 744 points and 490 quadrilaterals have  been used. Finally the blank holder has been discretized with 155 points and  120 quadrilaterals. Figure [[#img-11|11]].a shows a perspective of the tools and  Figure [[#img-11|11]].a the final deformed sheet. The sheet has been modelled  with 6370 3-node triangles and 3284 nodes (9274 DOFs).
1640
1641
<div id='img-11'></div>
1642
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1643
|-
1644
|[[Image:Draft_Samper_340530075-test-fig11.png|300px|Deep drawing of a  circular sheet. (a)geometry of the tools, (b)final deformed shape of the  sheet.]]
1645
|- style="text-align: center; font-size: 75%;"
1646
| colspan="1" | '''Figure 11:''' Deep drawing of a  circular sheet. (a)geometry of the tools, (b)final deformed shape of the  sheet.
1647
|}
1648
1649
The reported results are associated to three different meridiands: '''A'''  at 90<math display="inline">^{0}</math> of the rolling direction, '''B''': at 45<math display="inline">^{0}</math> of the rolling  direction and '''C''': in the rolling direction. The numerical results are  compared with a set of experimental values reported by Thyssen Krupp Stahl AG  (who proposed the benchmark and supplied the sheet samples) <span id='citeF-17'></span>[[#cite-17|[17]]]. It must be noted  that there was a great dispersion in the experimental data send to NUMISHEET,  so it is difficult to extract conclusions from a unique comparison. In any case the  experimental data supplied  seems to be consistent enought to have an  idea if the numerical model gives reasonable results.
1650
1651
The first element of comparison associated to the planar anisotropy, are the  different displacements (draw in) of the points in the boundary of the sheet. Table [[#table-5|5]] shows the values obtained for the three reference meridians chosen.
1652
1653
1654
{|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
1655
|+ style="font-size: 75%;" |<span id='table-5'></span>Table. 5 Flange draw in at three reference meridians
1656
|- style="border-top: 2px solid;"
1657
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  &nbsp; 
1658
| colspan='3' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" |  Draw in [mm]
1659
|- style="border-top: 2px solid;"
1660
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Model 
1661
| style="border-left: 2px solid;border-right: 2px solid;" | Section A 
1662
| style="border-left: 2px solid;border-right: 2px solid;" | Section B 
1663
| style="border-left: 2px solid;border-right: 2px solid;" | Section C
1664
|- style="border-top: 2px solid;"
1665
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  CBST1-Non Associative 
1666
| style="border-left: 2px solid;border-right: 2px solid;" | 29.06 
1667
| style="border-left: 2px solid;border-right: 2px solid;" | 31.91 
1668
| style="border-left: 2px solid;border-right: 2px solid;" | 30.77
1669
|- style="border-top: 2px solid;"
1670
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  CBST1-Associative 
1671
| style="border-left: 2px solid;border-right: 2px solid;" | 27.08 
1672
| style="border-left: 2px solid;border-right: 2px solid;" | 31.36 
1673
| style="border-left: 2px solid;border-right: 2px solid;" | 28.95
1674
|- style="border-top: 2px solid;border-bottom: 2px solid;"
1675
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Experimental 
1676
| style="border-left: 2px solid;border-right: 2px solid;" | 30.75 
1677
| style="border-left: 2px solid;border-right: 2px solid;" | 32.30 
1678
| style="border-left: 2px solid;border-right: 2px solid;" | 30.00
1679
1680
|}
1681
1682
The values for the non-associative model are much like the experimental ones,  those obtained with the associative model are a little lower. The largest  draw-in is at the meridian at 45<math display="inline">^{o}</math> of rolling direction ( B). In the numerial simulation the draw-in at meridian  C is larger than in meridian  A, while in the experimental results   the latter are larger than the former. It should be reminded that other  numerical results presented at NUMISHEET and most of the numerical simulations  showed the same tendency of present values.
1683
1684
The experimental data supplied are the three logarithmic principal strains  associated to the three meridians defined above. Here one strain at each  meridian has been choosen for comparison. Figure [[#img-12|12]].a shows  the meridian strain (<math display="inline">E_{1}</math>) along meridian A (direction transversal to rolling); Figure [[#img-12|12]].b shows the  circunferencial strain (<math display="inline">E_{2}</math>) along meridian  B and in Figure [[#img-12|12]].c the thickness strain (<math display="inline">E_{3}</math>)  along meridian C (rolling direction) has  been plotted.
1685
1686
The experimental data show a saw-tooth profile (specially for the thickness strain)  that is difficult to accept. These values are therefore not reliable point-wise but  as a whole. Based on this fact and on the above mentioned differences between  experimental values, it can be said that the present simulations agree quite  well with the experimental values, specially when the non-associative  plasticity model was used.
1687
1688
In the central part of the sheet the experimental data gives strains lower  than the numerical simulation. It can be said that experimentally the sheet  has better drawability. In the external zone there are less differences for in  plane strains but for thickness strains they are larger. Experimentally  thickness strains are more or less uniform between <math display="inline">-0.1</math> and <math display="inline">-0.2</math> what  differs from present numerical simulations and also from other numerical  simulations and experimental data presented at NUMISHEET <span id='citeF-17'></span>[[#cite-17|[17]]].
1689
1690
Finally Figure [[#img-12|12]] shows resistance to punch as a funtion of punch  travel. In the final part simulation forces are smaller than experimental  values. This may be due to an incorrect definition of friction with the tools.
1691
1692
<div id='img-12'></div>
1693
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1694
|-
1695
|[[Image:Draft_Samper_340530075-test-fig12.png|500px|Results from NUMISHEET  benchmark. (a) Meridional strain E₁ along transversal direction to  rolling (meridian A); (b) hoop strain E₂ along direction at 45⁰ of rolling (meridian  B); (c) Thickness strain E₃ along rolling direction  (meridian C); (d) Punch force versus Punch  travel; ]]
1696
|- style="text-align: center; font-size: 75%;"
1697
| colspan="1" | '''Figure 12:''' Results from NUMISHEET  benchmark. (a) Meridional strain <math>E_{1}</math> along transversal direction to  rolling (meridian A); (b) hoop strain <math>E_{2}</math> along direction at 45<math>^{0}</math> of rolling (meridian  B); (c) Thickness strain <math>E_{3}</math> along rolling direction  (meridian C); (d) Punch force versus Punch  travel; 
1698
|}
1699
1700
From the present results it is difficult to draw definite conclusions about what  should be improved in the numerical model due to all the factors involved,  including the constitutive model, the thin shell theory, the contact  treatment, the reliability of the experimental data of both the material  characterization and the deep drawing, etc..
1701
1702
==9 Conclusions==
1703
1704
A membrane and bending non-conforming rotation-free shell finite element has been presented.  The element pass the membrane patch test and the numerical  experiments performed do not show problems. From the bending point of view the  element is a little stiffer than the original BST element, but presents a smoother  continuity of displacement gradients. Convergence rate in membrane dominated  problems is similar to the linear strain triangle. Using three integration  points the elements is vulnerable to membrane locking for initially curved  surfaces. This locking effect dissapears when only one integration point is used, which the  usual case for elastic-plastic problems and for any large scale simulation.
1705
1706
==Acknowledgments==
1707
1708
The first author is a member of the scientific staff of the Science Research  Council of Argentina (CONICET). The support provided by grants of CONICET and  Agencia Córdoba Ciencia S.E. is gratefully acknowledged.
1709
1710
==BIBLIOGRAPHY==
1711
1712
<div id="cite-1"></div>
1713
'''[[#citeF-1|[1]]]''' F.G. Flores, E. Oñate and F. Zárate “New assumed  strain triangles for non-linear shell analysis”, ''Computational  Mechanics'', '''17''', 107-114, 1995.
1714
1715
<div id="cite-2"></div>
1716
'''[[#citeF-2|[2]]]''' 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.
1717
1718
<div id="cite-3"></div>
1719
'''[[#citeF-3|[3]]]'''  O.C. Zienkiewicz and R.L. Taylor. ''The finite element method. Solid Mechanics''. Vol II, Butterworth-Heinemann, 2000.
1720
1721
<div id="cite-4"></div>
1722
'''[[#citeF-4|[4]]]''' E. Oñate and F. Zárate, “Rotation-free plate and  shell triangles”, ''Int. J. Num. Meth. Engng'', '''47''', pp.  557&#8211;603, 2000.
1723
1724
<div id="cite-5"></div>
1725
'''[[#citeF-5|[5]]]''' 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, págs. 2039-2072.
1726
1727
<div id="cite-6"></div>
1728
'''[[#citeF-6|[6]]]''' D.Y. Yang, D.W. Jung, L.S. Song, D.J. Yoo and J.H. Lee,  “Comparative investigation into implicit, explicit and iterative  implic/explicit schemes for simulation of sheet metal forming processes”  ''NUMISHEET'93'', A. Makinouchi, E. Nakamachi, E. Oñate and R.H.  Wagoner (Eds.), RIKEN, pp. 35&#8211;42, Tokyo, 1993.
1729
1730
<div id="cite-7"></div>
1731
'''[[#citeF-7|[7]]]''' 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.
1732
1733
<div id="cite-8"></div>
1734
'''[[#citeF-8|[8]]]''' 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 Technoloy'', Ed. M.J.M. Barata Marques, 18th IDDRG  Biennial Congress, Lisbon, 1994.
1735
1736
<div id="cite-9"></div>
1737
'''[[#citeF-9|[9]]]''' Rojek J., Oñate E. and Postek E '''.''',  Application of explicit FE codes to simulation of sheet and  bulk metal forming processes ''Journal of the  Materials Processing Thecnology'', '''41''', págs. 620-627, 1998.
1738
1739
<div id="cite-10"></div>
1740
'''[[#citeF-10|[10]]]''' F.G. Flores F.G. 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.
1741
1742
<div id="cite-11"></div>
1743
'''[[#citeF-11|[11]]]''' R. Hill, “A Theory of the Yielding and Plastic Flow of  Anisotropic Metals”, ''Proc. Royal Society London'', vol. '''A193''',  pp. 281, 1948.
1744
1745
<div id="cite-12"></div>
1746
'''[[#citeF-12|[12]]]''' 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-584, 1972.
1747
1748
<div id="cite-13"></div>
1749
'''[[#citeF-13|[13]]]''' H.C. Huang, ''Static and Dynamic Analysis of Plates and  Shells'', page 40, Springer-Verlag, Berlin, 1989.
1750
1751
<div id="cite-14"></div>
1752
'''[[#citeF-14|[14]]]''' STAMPACK, ''A General Finite Element System for  Sheet Stamping and Forming Problems'', Quantech ATZ, Barcelona, Spain,  2003 (www.quantech.es).
1753
1754
<div id="cite-15"></div>
1755
'''[[#citeF-15|[15]]]''' A. Needleman, Inflation of spherical rubber ballons,  ''Int. J. of Solids and Structures'', 13, 409&#8211;421, 1977.
1756
1757
<div id="cite-16"></div>
1758
'''[[#citeF-16|[16]]]''' 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.
1759
1760
<div id="cite-17"></div>
1761
'''[[#citeF-17|[17]]]''' NUMISHEET'99, ''Fourth International Conference and  Workshop on Numerical Simulation of 3D Sheet Forming Processes, NUMISHEET'99'',  J.C. Gelin and P. Picart (Eds.), BURS Edition, Besancon, France, 1999.
1762

Return to Flores et al 2018c.

Back to Top

Document information

Published on 01/01/2005

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

Document Score

0

Times cited: 42
Views 126
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?