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

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


You can view and copy the source of this page.

x
 
1
Published in ''Computational Plasticity'', E. Oñate and R. Owen (Eds.), pp. 239-265, 2007<br />
2
DOI: 10.1007/978-1-4020-6577-4_14
3
==Abstract==
4
5
An enhanced rotation-free three node triangular shell element (termed EBST) is presented. The element formulation is based on a quadratic interpolation of the geometry in terms of the six  nodes of a patch of four triangles associated to each triangular element. This allows to compute an assumed constant curvature field and an assumed linear membrane strain field which improves the in-plane behaviour of the  element. A simple and economic version of the element using a single integration point is presented. The implementation of the element into an explicit dynamic scheme is described. The efficiency and accuracy of the EBST element  and the explicit dynamic scheme are demonstrated in many examples of application including the analysis of a cylindrical panel under impulse loading  and  sheet metal stamping problems.
6
7
==1 Introduction==
8
9
Triangular shell elements are very useful for the solution of large scale shell problems occurring in many practical engineering situations. Typical examples are the analysis of shell roofs under static and dynamic loads, sheet stamping processes, vehicle dynamics and crash-worthiness situations. Many of these problems involve high geometrical and material non linearities and time changing frictional contact conditions. These difficulties are usually increased by the need of discretizing complex geometrical shapes. Here the use of shell triangles and non-structured meshes becomes a critical necessity. Despite recent advances in the field <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-1|[1]]-<span id='citeF-6'></span>[[#cite-6|6]]] there are not so many simple shell triangles which are capable of accurately modelling the deformation of a shell structure under arbitrary loading conditions.
10
11
A promising line to derive simple shell triangles is to use the nodal displacements as the only unknowns for describing the shell kinematics. This idea goes back to the original attempts to solve thin plate bending problems using finite difference schemes with the deflection as the only nodal variable <span id='citeF-7'></span>[[#cite-7|[7]]-<span id='citeF-9'></span>[[#cite-9|9]]].
12
13
In past years some authors have derived a number of thin plate and shell triangular elements free of rotational degrees of freedom (d.o.f.) based on Kirchhoff's theory <span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-25'></span>[[#cite-10|[10]]-<span id='citeF-26'></span>[[#cite-26|26]]]. In essence all methods attempt to express the curvatures field over an element in terms of the displacements of a collection of nodes belonging to a patch of adjacent elements. Oñate and Cervera <span id='citeF-14'></span>[[#cite-14|[14]]] proposed a general procedure of this kind combining finite element and finite volume concepts for deriving thin plate triangles and quadrilaterals with the deflection as the only nodal variable and presented a simple and competitive rotation-free three d.o.f. triangular element termed BPT (for Basic Plate Triangle). These ideas were extended in <span id='citeF-20'></span>[[#cite-20|[20]]] to derive a number of rotation-free thin plate and shell triangles. The basic ingredients of the method are a mixed Hu-Washizu formulation, a standard discretization into three node triangles, a linear finite element interpolation of the displacement field within each triangle and a finite volume type approach for computing constant curvature and bending moment fields within appropriate non-overlapping control domains. The so called cell-centered and cell-vertex triangular domains yield different families of rotation-free plate and shell triangles. Both the BPT plate element and its extension to shell analysis (termed BST for Basic Shell Triangle) can be derived from the cell-centered formulation. Here the control domain is an individual triangle. The constant curvatures field within a triangle is computed in terms of the displacements of the six nodes belonging to the four elements patch formed by the chosen triangle and the three adjacent triangles. The cell-vertex approach yields a different family of rotation-free plate and shell triangles. Details of the derivation of both rotation-free triangular shell element families can be found in <span id='citeF-20'></span>[[#cite-20|[20]]].
14
15
An extension of the BST element to the non linear analysis of shells was implemented in an explicit dynamic code by Oñate ''et al.'' <span id='citeF-25'></span>[[#cite-25|[25]]] using an updated Lagrangian formulation and a hypo-elastic constitutive model. Excellent numerical results were obtained for non linear dynamics of shells involving frictional contact situations and sheet stamping problems <span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-25'></span>[[#cite-17|[17]],[[#cite-18|18]],[[#cite-19|19]],[[#cite-25|25]]].
16
17
A large strain formulation for the BST element using a total Lagrangian description was presented by Flores and Oñate <span id='citeF-23'></span>[[#cite-23|[23]]]. A recent extension of this formulation is based on a quadratic interpolation of the geometry of the patch formed by the BST element and the three adjacent triangles <span id='citeF-26'></span>[[#cite-26|[26]]]. This yields a linear displacement gradient field over the element from which linear membrane strains and  constant curvatures  can be computed within the BST element.
18
19
In this chapter an enhanced version of the BST element (termed EBST element) is derived using an assumed linear field for the membrane strains and an assumed constant curvature field. Both assumed fields are obtained from the quadratic interpolation of the patch geometry following the ideas presented in <span id='citeF-26'></span>[[#cite-26|[26]]]. Details of the element formulation are given. An efficient version of the  EBST element using one single quadrature point for integration of the tangent matrix is  presented. An explicit  scheme adequate for dynamic analysis is  briefly described.
20
21
The efficiency and accuracy of the EBST element is validated in a number of examples of application including the non linear analysis of a cylindrical shell under an impulse loading and practical sheet stamping problems.
22
23
==2 Basic Thin Shell Equations Using a Total Lagrangian Formulation==
24
25
===2.1 Shell Kinematics===
26
27
A summary of the most relevant hypothesis related to the kinematic behaviour of a thin shell are presented. Further details may be found in the wide literature dedicated to this field <span id='citeF-8'></span>[[#cite-8|[8]],<span id='citeF-9'></span>[[#cite-9|9]]].
28
29
Consider a shell with undeformed middle surface occupying the domain <math display="inline">\Omega ^{0}</math> in <math display="inline">R^{3}</math> with a boundary <math display="inline">\Gamma ^{0}</math>. At each point of the middle surface a thickness <math display="inline">h^{0}</math> is defined. The positions <math display="inline">\mathbf{x}^{0}</math> and <math display="inline">\mathbf{x}</math> of a point in the undeformed and the deformed configurations can be respectively written as a function of the coordinates of the middle surface <math display="inline">{\boldsymbol \varphi }</math> and the normal <math display="inline">\mathbf{t}_{3}</math> at the point as
30
31
{| class="formulaSCP" style="width: 100%; text-align: left;" 
32
|-
33
| 
34
{| style="text-align: left; margin:auto;width: 100%;" 
35
|-
36
| style="text-align: center;" | <math>\mathbf{x}^{0}\left( \xi _{1},\xi _{2},\zeta \right)    ={\boldsymbol \varphi }^{0}\left( \xi _{1},\xi _{2}\right) +\lambda \mathbf{t}_{3}^{0}</math>
37
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
38
|-
39
| style="text-align: center;" | <math> \mathbf{x}\left( \xi _{1},\xi _{2},\zeta \right)    ={\boldsymbol \varphi }\left( \xi  _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}</math>
40
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
41
|}
42
|}
43
44
where <math display="inline">\xi _{1},\xi _{2}</math> are arc-length curvilinear principal coordinates defined over the middle surface of the shell and <math display="inline">\zeta </math> is the distance from the point to the middle surface in the undeformed configuration. The product <math display="inline">\zeta \lambda </math> is the distance from the point to the middle surface measured on the deformed configuration. The parameter <math display="inline">\lambda </math> relates the thickness at the present and initial configurations as:
45
46
{| class="formulaSCP" style="width: 100%; text-align: left;" 
47
|-
48
| 
49
{| style="text-align: left; margin:auto;width: 100%;" 
50
|-
51
| style="text-align: center;" | <math>\lambda =\frac{h}{h^{0}}</math>
52
|}
53
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
54
|}
55
56
This approach implies a constant strain in the normal direction. Parameter <math display="inline">\lambda </math> will not be considered as an independent variable  and will be computed from purely geometrical considerations (''isochoric'' behaviour) via a staggered iterative update. Besides this, the usual plane stress condition of thin shell theory will be adopted.
57
58
A convective system is computed at each point as
59
60
{| class="formulaSCP" style="width: 100%; text-align: left;" 
61
|-
62
| 
63
{| style="text-align: left; margin:auto;width: 100%;" 
64
|-
65
| style="text-align: center;" | <math>\mathbf{g}_{i}\left( \mathbf{\xi }\right) =\frac{\partial \mathbf{x}}{\partial \xi _{i}}\qquad i=1,2,3 </math>
66
|}
67
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
68
|}
69
70
{| class="formulaSCP" style="width: 100%; text-align: left;" 
71
|-
72
| 
73
{| style="text-align: left; margin:auto;width: 100%;" 
74
|-
75
| style="text-align: center;" | <math>\mathbf{g}_{\alpha }\left( \mathbf{\xi }\right)    =\frac{\partial \left( \mathbf{\boldsymbol \varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \xi _{\alpha }}={\boldsymbol \varphi }_{^{\prime }\alpha }+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }\alpha }\quad \alpha=1,2</math>
76
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
77
|-
78
| style="text-align: center;" | <math> \mathbf{g}_{3}\left( \mathbf{\xi }\right)    =\frac{\partial \left( \mathbf{\boldsymbol \varphi }\left( \xi _{1},\xi _{2}\right) +\zeta \lambda \mathbf{t}_{3}\right) }{\partial \zeta }=\lambda \mathbf{t}_{3}</math>
79
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
80
|}
81
|}
82
83
This can be particularized for the points on the middle surface as
84
85
{| class="formulaSCP" style="width: 100%; text-align: left;" 
86
|-
87
| 
88
{| style="text-align: left; margin:auto;width: 100%;" 
89
|-
90
| style="text-align: center;" | <math>\mathbf{a}_{\alpha }    =\mathbf{g}_{\alpha }\left( \zeta=0\right) ={\boldsymbol \varphi  }_{^{\prime }\alpha }</math>
91
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
92
|-
93
| style="text-align: center;" | <math> \mathbf{a}_{3}    =\mathbf{g}_{3}\left( \zeta=0\right) =\lambda  \mathbf{t}_{3}</math>
94
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
95
|}
96
|}
97
98
The covariant (first fundamental form) metric tensor of the middle surface is
99
100
<span id="eq-9"></span>
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 } = {\boldsymbol \varphi }_{^{\prime }\alpha } \cdot  {\boldsymbol \varphi }_{^{\prime }\beta }  </math>
107
|}
108
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
109
|}
110
111
The Green-Lagrange strain vector of the middle surface points (membrane strains) is defined 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>{\boldsymbol \varepsilon }_{m}=[\varepsilon _{m_{11}},\varepsilon _{m_{12}},\varepsilon _{m_{12}}]^{T}</math>
119
|}
120
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
121
|}
122
123
with
124
125
<span id="eq-11"></span>
126
{| class="formulaSCP" style="width: 100%; text-align: left;" 
127
|-
128
| 
129
{| style="text-align: left; margin:auto;width: 100%;" 
130
|-
131
| style="text-align: center;" | <math>\varepsilon _{m_{ij}}=\frac{1}{2}(a_{ij}-a_{ij}^{0}) </math>
132
|}
133
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
134
|}
135
136
The curvatures (second fundamental form) of the middle surface are obtained by
137
138
{| class="formulaSCP" style="width: 100%; text-align: left;" 
139
|-
140
| 
141
{| style="text-align: left; margin:auto;width: 100%;" 
142
|-
143
| style="text-align: center;" | <math>\kappa _{\alpha \beta }=\frac{1}{2}\left( {\boldsymbol \varphi }_{^{\prime }\alpha }\cdot \mathbf{t}_{3^{\prime }\beta }+{\boldsymbol \varphi }_{^{\prime }\beta }\cdot  \mathbf{t}_{3^{\prime }\alpha }\right) =- \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{{\prime }\alpha \beta }\quad , \quad \alpha ,\beta=1,2 </math>
144
|}
145
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
146
|}
147
148
The deformation gradient tensor is
149
150
{| class="formulaSCP" style="width: 100%; text-align: left;" 
151
|-
152
| 
153
{| style="text-align: left; margin:auto;width: 100%;" 
154
|-
155
| style="text-align: center;" | <math>\mathbf{F=} [{\boldsymbol x}_{{\prime }1},{\boldsymbol x}_{{\prime }2},{\boldsymbol x}_{{\prime }3}]=\left[ \begin{array}{ccc}{\boldsymbol \varphi }_{^{\prime }1}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime  }1} & {\boldsymbol \varphi }_{^{\prime }2}+\zeta \left( \lambda \mathbf{t}_{3}\right) _{^{\prime }2} & \lambda \mathbf{t}_{3}\end{array} \right] </math>
156
|}
157
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
158
|}
159
160
The product <math display="inline">\mathbf{F}^{T}\mathbf{F=U}^{2}=\mathbf{C}</math> (where <math display="inline">\mathbf{U}</math> is the right stretch tensor, and <math display="inline">\mathbf{C}</math> the right Cauchy-Green deformation tensor) can be written as
161
162
<span id="eq-14"></span>
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}+2\kappa _{11}\zeta \lambda & a_{12}+2\kappa _{12}\zeta \lambda & 0\\ a_{12}+2\kappa _{12}\zeta \lambda & a_{22}+2\kappa _{22}\zeta \lambda & 0\\ 0 & 0 & \lambda ^{2}\end{array} \right] </math>
169
|}
170
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
171
|}
172
173
In the derivation of expression ([[#eq-14|14]]) the derivatives of the thickness ratio <math display="inline">\lambda _{^{\prime }a}</math> and the terms associated to <math display="inline">\zeta ^{2}</math> have been neglected.
174
175
Equation ([[#eq-14|14]]) shows that <math display="inline">\mathbf{U}^{2}</math> is not a unit tensor at the original configuration for curved surfaces (<math display="inline">\kappa _{ij}^{0}\neq{0}</math>). The changes of curvature of the middle surface are computed by
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;" | (15)
185
|}
186
187
Note that <math display="inline">\delta \chi _{ij}=\delta \kappa _{ij}</math>.
188
189
For computational convenience the following approximate expression (which is exact for initially flat surfaces) will be adopted
190
191
<span id="eq-16"></span>
192
{| class="formulaSCP" style="width: 100%; text-align: left;" 
193
|-
194
| 
195
{| style="text-align: left; margin:auto;width: 100%;" 
196
|-
197
| style="text-align: center;" | <math>\mathbf{U}^{2}=\left[ \begin{array}{ccc}a_{11}+2\chi _{11}\zeta \lambda & a_{12}+2\chi _{12}\zeta \lambda & 0\\ a_{12}+2\chi _{12}\zeta \lambda & a_{22}+2\chi _{22}\zeta \lambda & 0\\ 0 & 0 & \lambda ^{2}\end{array} \right]  </math>
198
|}
199
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
200
|}
201
202
This expression is useful to compute different Lagrangian strain measures. An advantage of these measures is that they are associated to material fibres, what makes it easy to take into account material anisotropy. It is also useful to compute the eigen decomposition of <math display="inline">\mathbf{U}</math> as
203
204
{| class="formulaSCP" style="width: 100%; text-align: left;" 
205
|-
206
| 
207
{| style="text-align: left; margin:auto;width: 100%;" 
208
|-
209
| style="text-align: center;" | <math>\mathbf{U=}\sum _{\alpha=1}^{3}\lambda _{\alpha } \mathbf{r}_{\alpha }\otimes \mathbf{r}_{\alpha }</math>
210
|}
211
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
212
|}
213
214
where <math display="inline">\lambda _{\alpha }</math> and <math display="inline">\mathbf{r}_{\alpha }</math> are the eigenvalues and eigenvectors of <math display="inline">\mathbf{U}</math>.
215
216
The resultant stresses  (axial forces and moments) are obtained by integrating across the original thickness the second Piola-Kirchhoff stress vector <math display="inline">{ \boldsymbol \sigma }</math> using the actual distance to the middle surface for  evaluating the bending moments. This gives
217
218
<span id="eq-18"></span>
219
{| class="formulaSCP" style="width: 100%; text-align: left;" 
220
|-
221
| 
222
{| style="text-align: left; margin:auto;width: 100%;" 
223
|-
224
| style="text-align: center;" | <math>{\boldsymbol \sigma }_{m}\equiv \lbrack N_{11},N_{22},N_{12}]^{T}=\int _{h^{0}}{\boldsymbol \sigma }d\zeta </math>
225
|}
226
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
227
|}
228
229
<span id="eq-19"></span>
230
{| class="formulaSCP" style="width: 100%; text-align: left;" 
231
|-
232
| 
233
{| style="text-align: left; margin:auto;width: 100%;" 
234
|-
235
| style="text-align: center;" | <math>{\boldsymbol \sigma }_{b}\equiv \lbrack M_{11},M_{22},M_{12}]^{T}=\int _{h^{0}}{\boldsymbol \sigma  }\lambda \zeta  d\zeta </math>
236
|}
237
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
238
|}
239
240
With these values the virtual work can be written as
241
242
<span id="eq-20"></span>
243
{| class="formulaSCP" style="width: 100%; text-align: left;" 
244
|-
245
| 
246
{| style="text-align: left; margin:auto;width: 100%;" 
247
|-
248
| style="text-align: center;" | <math>\int \int _{A^{0}}\left[ \delta{\boldsymbol \varepsilon }_{m}^{T}{\boldsymbol \sigma }_{m}+\delta{\boldsymbol \kappa  }^{T}{\boldsymbol \sigma }_{b}\right] dA=\int \int _{A^{0}}\delta \mathbf{u}^{T}\mathbf{t}dA </math>
249
|}
250
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
251
|}
252
253
where <math display="inline">\delta \mathbf{u}</math> are virtual displacements, <math display="inline">\delta{\boldsymbol \varepsilon }_{m}</math> is the virtual Green-Lagrange membrane strain vector, <math display="inline">\delta{\boldsymbol \kappa }</math> are the virtual curvatures and <math display="inline">\mathbf{t}</math> are the surface loads. Other load types can be easily included into ([[#eq-20|20]]).
254
255
===2.2 Constitutive Models===
256
257
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
258
259
<span id="eq-21"></span>
260
{| class="formulaSCP" style="width: 100%; text-align: left;" 
261
|-
262
| 
263
{| style="text-align: left; margin:auto;width: 100%;" 
264
|-
265
| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=}\left[ \begin{array}{ccc}\varepsilon _{11} & \varepsilon _{21} & 0\\ \varepsilon _{12} & \varepsilon _{22} & 0\\ 0 & 0 & \varepsilon _{33}\end{array} \right] =\sum _{\alpha=1}^{3}\ln \left( \lambda _{\alpha }\right) \mathbf{r}_{\alpha }\otimes \mathbf{r}_{\alpha } </math>
266
|}
267
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
268
|}
269
270
For the type of problems dealt within the paper we use an elastic-plastic material associated to thin rolled metal sheets.
271
272
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
273
274
<span id="eq-22"></span>
275
{| class="formulaSCP" style="width: 100%; text-align: left;" 
276
|-
277
| 
278
{| style="text-align: left; margin:auto;width: 100%;" 
279
|-
280
| style="text-align: center;" | <math>\mathbf{E}_{\ln }\mathbf{=E}_{\ln }^{e}+\mathbf{E}_{\ln }^{p} </math>
281
|}
282
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
283
|}
284
285
A  linear relationship between the (plane) Hencky stresses and the logarithmic elastic strains is  chosen giving
286
287
<span id="eq-23"></span>
288
{| class="formulaSCP" style="width: 100%; text-align: left;" 
289
|-
290
| 
291
{| style="text-align: left; margin:auto;width: 100%;" 
292
|-
293
| style="text-align: center;" | <math>\mathbf{T}=\mathbf{H} \mathbf{E}_{\ln }^{e} </math>
294
|}
295
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
296
|}
297
298
where <math display="inline">\boldsymbol H</math> is the constitutive matrix. The constitutive equations are integrated using a standard return algorithm. The following Mises-Hill <span id='citeF-30'></span>[[#cite-30|[30]]] yield function with non-linear isotropic hardening is chosen
299
300
{| class="formulaSCP" style="width: 100%; text-align: left;" 
301
|-
302
| 
303
{| style="text-align: left; margin:auto;width: 100%;" 
304
|-
305
| style="text-align: center;" | <math>\left( G+H\right) \;T_{11}^{2}+\left( F+H\right) \;T_{22}^{2}-2H\;T_{11}T_{22}+2N\;T_{12}^{2}=\sigma _0\left(e_{0}+e^{p}\right) ^{n}</math>
306
|}
307
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
308
|}
309
310
where <math display="inline">F, G, H</math> and <math display="inline">N</math> define the non-isotropic shape of the yield surface and the parameters <math display="inline">\sigma _{0}</math>, <math display="inline">e_{0}</math> and <math display="inline">n</math> define its size as a function of the effective plastic strain <math display="inline">e^{p}</math>.
311
312
The simple Mises-Hill yield function  allows, as a first approximation, to treat rolled thin metal sheets with planar and transversal anisotropy.
313
314
The Hencky stress tensor <math display="inline">\mathbf{T}</math> can be easily particularized for the plane stress case.
315
316
We define the rotated Hencky and second Piola-Kirchhoff stress tensors as
317
318
<span id="eq-25"></span>
319
<span id="eq-26"></span>
320
{| class="formulaSCP" style="width: 100%; text-align: left;" 
321
|-
322
| 
323
{| style="text-align: left; margin:auto;width: 100%;" 
324
|-
325
| style="text-align: center;" | <math>\mathbf{T}_{L}    =\mathbf{R}_{L}^{T}\;\mathbf{T\;R}_{L}</math>
326
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
327
|-
328
| style="text-align: center;" | <math> \mathbf{S}_{L}    =\mathbf{R}_{L}^{T}\;\mathbf{S\;R}_{L}</math>
329
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
330
|}
331
|}
332
333
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
334
335
{| class="formulaSCP" style="width: 100%; text-align: left;" 
336
|-
337
| 
338
{| style="text-align: left; margin:auto;width: 100%;" 
339
|-
340
| style="text-align: center;" | <math>\mathbf{R}_{L}=\left[ \begin{array}{ccc}\mathbf{r}_{1}\quad ,& \mathbf{r}_{2} \quad ,& \mathbf{r}_{3}\end{array} \right] </math>
341
|}
342
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
343
|}
344
345
The relationship between the rotated Hencky and Piola-Kirchhoff stresses is <math display="inline">\left(\alpha , \beta=1,2 \right)</math>
346
347
{| class="formulaSCP" style="width: 100%; text-align: left;" 
348
|-
349
| 
350
{| style="text-align: left; margin:auto;width: 100%;" 
351
|-
352
| style="text-align: center;" | <math>\left[ S_{L}\right] _{\alpha \alpha }    =\frac{1}{\lambda _{\alpha }^{2}}\left[ T_{L}\right] _{\alpha \alpha }</math>
353
|-
354
| 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>
355
|}
356
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
357
|}
358
359
The second Piola-Kirchhoff stress tensor can be computed by
360
361
{| class="formulaSCP" style="width: 100%; text-align: left;" 
362
|-
363
| 
364
{| style="text-align: left; margin:auto;width: 100%;" 
365
|-
366
| style="text-align: center;" | <math>\mathbf{S=}\sum _{\alpha=1}^{2}\sum _{\beta=1}^{2}\left[ S_{L}\right] _{\alpha \beta } \mathbf{r}_{\alpha }\otimes \mathbf{r}_{\beta }</math>
367
|}
368
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
369
|}
370
371
The second Piola-Kirchhoff stress vector <math display="inline">{\boldsymbol \sigma }</math> used in Eqs.([[#eq-18|18]]&#8211;[[#eq-19|19]]) can be readily extracted from the <math display="inline">\mathbf{S}</math> tensor.
372
373
==3 Enhanced Basic Shell Triangle==
374
375
The main features of the element formulation (termed EBST for Enhanced Basic Shell Triangle) are the following:
376
377
<ol>
378
379
<li>The geometry of the patch formed by an element and the three adjacent elements is ''quadratically interpolated'' from the position of the six nodes in the patch (Fig. [[#img-1|1]]). </li>
380
381
<li>The membrane strains are assumed to vary ''linearly'' within the central triangle and are expressed in terms of the (continuous) values of the deformation gradient at the mid side points of the triangle. </li>
382
383
<li>An assumed ''constant curvature'' field within the central triangle is chosen. This is computed in terms of the values of the (continuous) deformation gradient at the mid side points. </li>
384
385
</ol>
386
387
Details of the derivation of the EBST element are given below.
388
389
===3.1 Definition of the Element Geometry and Computation of Membrane Strains===
390
391
A  quadratic approximation of the geometry of the four elements patch is chosen using the position of the six nodes in the patch. It is useful to define the patch in the isoparametric space using the nodal positions given in the Table [[#table-1|1]] (see also Fig.&nbsp;[[#img-1|1]]).
392
393
394
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
395
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Isoparametric coordinates of the six nodes in the patch of Fig. [[#img-1|1]]
396
|- style="border-top: 2px solid;"
397
| style="border-left: 2px solid;border-right: 2px solid;" |  
398
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
399
| style="border-left: 2px solid;border-right: 2px solid;" | 2 
400
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
401
| style="border-left: 2px solid;border-right: 2px solid;" | 4 
402
| style="border-left: 2px solid;border-right: 2px solid;" | 5 
403
| style="border-left: 2px solid;border-right: 2px solid;" | 6
404
|- style="border-top: 2px solid;"
405
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\xi </math> 
406
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
407
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
408
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
409
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
410
| style="border-left: 2px solid;border-right: 2px solid;" | -1 
411
| style="border-left: 2px solid;border-right: 2px solid;" | 1
412
|- style="border-top: 2px solid;border-bottom: 2px solid;"
413
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\eta </math> 
414
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
415
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
416
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
417
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
418
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
419
| style="border-left: 2px solid;border-right: 2px solid;" | -1
420
421
|}
422
423
The quadratic interpolation is defined by
424
425
<span id="eq-30"></span>
426
{| class="formulaSCP" style="width: 100%; text-align: left;" 
427
|-
428
| 
429
{| style="text-align: left; margin:auto;width: 100%;" 
430
|-
431
| style="text-align: center;" | <math>{\boldsymbol \varphi }=\sum _{i=1}^{6}N_{i}{\boldsymbol \varphi }_{i}</math>
432
|}
433
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
434
|}
435
436
with (<math display="inline">\zeta=1-\xi-\eta</math>)
437
438
{| class="formulaSCP" style="width: 100%; text-align: left;" 
439
|-
440
| 
441
{| style="text-align: left; margin:auto;width: 100%;" 
442
|-
443
| style="text-align: center;" | <math>\begin{array}{ccc}N_{1}=\zeta{+\xi}\eta &  & N_{4}=\frac{\zeta }{2}\left( \zeta{-1}\right) \\ N_{2}=\xi{+\eta}\zeta &  & N_{5}=\frac{\xi }{2}\left( \xi{-1}\right) \\ N_{3}=\eta{+\zeta}\xi &  & N_{6}=\frac{\eta }{2}\left( \eta{-1}\right) \end{array} </math>
444
|}
445
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
446
|}
447
448
This interpolation allows to computing the displacement gradients at selected points in order to use an assumed strain approach. The computation of the gradients is performed at the mid side points of the central element of the patch denoted by <math display="inline">G_{1}</math>, <math display="inline">G_{2}</math> and <math display="inline">G_{3}</math> in Fig. [[#img-1|1]]. This choice has the following advantages.
449
450
<div id='img-1a'></div>
451
<div id='img-1b'></div>
452
<div id='img-1'></div>
453
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
454
|-
455
|[[Image:Draft_Samper_105745998-fig1.png|300px|]]
456
|[[Image:Draft_Samper_105745998-fig2.png|300px|]]
457
|- style="text-align: center; font-size: 75%;"
458
| (a) 
459
| (b) 
460
|- style="text-align: center; font-size: 75%;"
461
| colspan="2" | '''Figure 1:''' (a) Patch of three node triangular elements including the central   triangle (M) and three adjacent triangles (1, 2 and 3); (b) Patch of elements   in the isoparametric space
462
|}
463
464
* Gradients at the three mid side points depend only on the nodes belonging to the two elements adjacent to each side. This can be easily verified by sampling the derivatives of the shape functions at each mid-side point.
465
466
* When gradients are computed at the common mid-side point of two adjacent elements, the same values are obtained, as the coordinates of the same four points are used. This in practice means that the gradients at the mid-side points are independent of the element where they are computed. A side-oriented implementation of the finite element will therefore lead to a unique evaluation of the gradients per side.
467
468
The Cartesian derivatives of the shape functions are computed at the original configuration by the standard expression
469
470
{| class="formulaSCP" style="width: 100%; text-align: left;" 
471
|-
472
| 
473
{| style="text-align: left; margin:auto;width: 100%;" 
474
|-
475
| style="text-align: center;" | <math>\left[ \begin{array}{c}N_{i,1}\\ N_{i,2}\end{array} \right] =\mathbf{J}^{-1}\left[ \begin{array}{c}N_{i,\xi } \\ N_{i,\eta }\end{array} \right] </math>
476
|}
477
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
478
|}
479
480
where the Jacobian matrix at the original configuration is
481
482
{| class="formulaSCP" style="width: 100%; text-align: left;" 
483
|-
484
| 
485
{| style="text-align: left; margin:auto;width: 100%;" 
486
|-
487
| style="text-align: center;" | <math>\mathbf{J=}\left[ \begin{array}{cc}\mathbf{\boldsymbol \varphi }_{^{\prime }\xi }^{0}\cdot \mathbf{t}_{1} & \mathbf{\boldsymbol \varphi  }_{^{\prime }\eta }^{0}\cdot \mathbf{t}_{1}\\ \mathbf{\boldsymbol \varphi }_{^{\prime }\xi }^{0}\cdot \mathbf{t}_{2} & \mathbf{\boldsymbol \varphi  }_{^{\prime }\eta }^{0}\cdot \mathbf{t}_{2}\end{array} \right] </math>
488
|}
489
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
490
|}
491
492
The deformation gradients on the middle surface, associated to an arbitrary spatial Cartesian system and to the material cartesian system defined on the middle surface are related by
493
494
{| class="formulaSCP" style="width: 100%; text-align: left;" 
495
|-
496
| 
497
{| style="text-align: left; margin:auto;width: 100%;" 
498
|-
499
| style="text-align: center;" | <math>\left[ {\boldsymbol \varphi }_{^{\prime }1},\mathbf{\boldsymbol \varphi }_{^{\prime }2}\right] =\left[ \mathbf{\boldsymbol \varphi }_{^{\prime }\xi },\mathbf{\boldsymbol \varphi }_{^{\prime }\eta }\right]  \mathbf{J}^{-1}</math>
500
|}
501
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
502
|}
503
504
The membrane strains within the central triangle are obtained using a linear assumed strain field <math display="inline">\hat{\boldsymbol \varepsilon }_{m}</math>, i.e.
505
506
{| class="formulaSCP" style="width: 100%; text-align: left;" 
507
|-
508
| 
509
{| style="text-align: left; margin:auto;width: 100%;" 
510
|-
511
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_{m}=\hat{\boldsymbol \varepsilon }_{m}</math>
512
|}
513
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
514
|}
515
516
with
517
518
<span id="eq-36"></span>
519
{| class="formulaSCP" style="width: 100%; text-align: left;" 
520
|-
521
| 
522
{| style="text-align: left; margin:auto;width: 100%;" 
523
|-
524
| style="text-align: center;" | <math>\hat{\boldsymbol \varepsilon }_{m}=(1-2\zeta ){\boldsymbol \varepsilon }_{m}^{1}+(1-2\xi ){\boldsymbol \varepsilon  }_{m}^{2}+(1-2\eta ){\boldsymbol \varepsilon }_{m}^{3}=\sum _{i=1}^{3}\bar{N}_{i}{\boldsymbol \varepsilon }_{m}^{i}</math>
525
|}
526
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
527
|}
528
529
where <math display="inline">{\boldsymbol \varepsilon }_{m}^{i}</math> are the membrane strains computed at the three mid side points <math display="inline">G_{i}</math> (<math display="inline">i=1,2,3</math>  see Fig. [[#img-1|1]]). In Eq.([[#eq-36|36]]) <math display="inline">\bar{N}_{1}=(1-2\zeta )</math>, etc.
530
531
The gradient at each mid side point is computed from the quadratic interpolation ([[#eq-30|30]]):
532
533
<span id="eq-37"></span>
534
{| class="formulaSCP" style="width: 100%; text-align: left;" 
535
|-
536
| 
537
{| style="text-align: left; margin:auto;width: 100%;" 
538
|-
539
| style="text-align: center;" | <math>\left( {\boldsymbol \varphi }_{^{\prime }\alpha }\right) _{G_{i}}={\boldsymbol \varphi }_{^{\prime  }\alpha }^{i}=\left[ \sum _{j=1}^{3}N_{j,\alpha }^{i}{\boldsymbol \varphi }_{j}\right] +N_{i+3,\alpha }^{i}{\boldsymbol \varphi }_{i+3}\quad ,\quad \alpha=1,2\quad ,\quad  i=1,2,3</math>
540
|}
541
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
542
|}
543
544
Substituting Eq.([[#eq-11|11]]) into ([[#eq-36|36]]) and using Eq.([[#eq-9|9]]) gives the membrane strain vector as
545
546
{| class="formulaSCP" style="width: 100%; text-align: left;" 
547
|-
548
| 
549
{| style="text-align: left; margin:auto;width: 100%;" 
550
|-
551
| style="text-align: center;" | <math>{\boldsymbol \varepsilon }_{m}=\sum _{i=1}^{3}\frac{1}{2}\bar{N}_{i}\left\{ \begin{array}{c}{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}-1\\ {\boldsymbol \varphi }_{^{\prime }2}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}-1\\ 2{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}\end{array} \right\} </math>
552
|}
553
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
554
|}
555
556
and the virtual membrane strains as
557
558
<span id="eq-39"></span>
559
{| class="formulaSCP" style="width: 100%; text-align: left;" 
560
|-
561
| 
562
{| style="text-align: left; margin:auto;width: 100%;" 
563
|-
564
| style="text-align: center;" | <math>\delta{\boldsymbol \varepsilon }_{m}=\sum _{i=1}^{3}\bar{N}_{i}\left\{ \begin{array}{c}{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \delta \mathbf{\boldsymbol \varphi }_{^{\prime }1}^{i}\\ {\boldsymbol \varphi }_{2}^{i}\cdot \delta \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}\\ \delta{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \mathbf{\boldsymbol \varphi }_{^{\prime }2}^{i}+{\boldsymbol \varphi }_{^{\prime }1}^{i}\cdot \delta \mathbf{\boldsymbol \varphi }_{2}^{i}\end{array} \right\} </math>
565
|}
566
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
567
|}
568
569
We note that the gradient at each mid side point <math display="inline">G_{i}</math> depends only on the coordinates of the three nodes of the central triangle and on those of an additional node in the patch, associated to the side <math display="inline">i</math> where the gradient is computed.
570
571
Combining Eqs.([[#eq-39|39]]), ([[#eq-37|37]]) and ([[#eq-30|30]]) gives
572
573
{| class="formulaSCP" style="width: 100%; text-align: left;" 
574
|-
575
| 
576
{| style="text-align: left; margin:auto;width: 100%;" 
577
|-
578
| style="text-align: center;" | <math>\delta{\boldsymbol \varepsilon }_{m}=\mathbf{B}_{m}\delta \mathbf{a}^{p}</math>
579
|}
580
| style="width: 5px;text-align: right;white-space: nowrap;" | (40.a)
581
|}
582
583
with
584
585
<span id="eq-40.b"></span>
586
{| class="formulaSCP" style="width: 100%; text-align: left;" 
587
|-
588
| 
589
{| style="text-align: left; margin:auto;width: 100%;" 
590
|-
591
| style="text-align: center;" | <math>\underset{18\times 1}{\delta \mathbf{a}^p} =[\delta \mathbf{u}_{1}^{T},\delta \mathbf{u}_{2}^{T},\delta \mathbf{u}_{3}^{T},\delta \mathbf{u}_{4}^{T},\delta \mathbf{u}_{5}^{T},\delta \mathbf{u}_{6}^{T}]^{T}</math>
592
|}
593
| style="width: 5px;text-align: right;white-space: nowrap;" | (40.b)
594
|}
595
596
where <math display="inline">\delta \mathbf{a}^{p}</math> is the patch displacement vector and <math display="inline">\mathbf{B}_{m}</math> is the membrane strain matrix. An explicit form of this matrix is given in <span id='citeF-26'></span>[[#cite-26|[26]]].
597
598
Note that the membrane strains within the EBST element are  a function of the displacements of the six patch nodes.
599
600
===3.2 Computation of Curvatures===
601
602
We will assume the following constant curvature field within each element
603
604
<span id="eq-41"></span>
605
{| class="formulaSCP" style="width: 100%; text-align: left;" 
606
|-
607
| 
608
{| style="text-align: left; margin:auto;width: 100%;" 
609
|-
610
| style="text-align: center;" | <math>\kappa _{\alpha \beta }=\hat{\kappa }_{\alpha \beta } </math>
611
|}
612
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
613
|}
614
615
where <math display="inline">\hat{\kappa }_{\alpha \beta }</math> is the assumed constant curvature field defined by
616
617
<span id="eq-42"></span>
618
{| class="formulaSCP" style="width: 100%; text-align: left;" 
619
|-
620
| 
621
{| style="text-align: left; margin:auto;width: 100%;" 
622
|-
623
| style="text-align: center;" | <math>\hat{\kappa }_{\alpha \beta }=-\frac{1}{A_{M}^{0}}\int _{A_{M}^{0}}\mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }\beta \alpha }\;dA^{0} </math>
624
|}
625
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
626
|}
627
628
where <math display="inline">A_{M}^{0}</math> is the area (in the original configuration) of the central element in the patch.
629
630
Substituting Eq.([[#eq-42|42]]) into ([[#eq-41|41]]) and integrating by parts the area integral gives the curvature vector within the element in terms of the following line integral
631
632
<span id="eq-43"></span>
633
{| class="formulaSCP" style="width: 100%; text-align: left;" 
634
|-
635
| 
636
{| style="text-align: left; margin:auto;width: 100%;" 
637
|-
638
| style="text-align: center;" | <math>{\boldsymbol \kappa }=\left\{ \begin{array}{c}\kappa _{11}\\ \kappa _{22}\\ 2\kappa _{12}\end{array} \right\} =\frac{1}{A_{M}^{0}}{\displaystyle \oint _{\Gamma _{M}^{0}}} \left[ \begin{array}{cc}-n_{1} & 0\\ 0 & -n_{2}\\ -n_{2} & -n_{1}\end{array} \right] \left[ \begin{array}{c}\mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }1}\\ \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }2}\end{array} \right] d\Gamma </math>
639
|}
640
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
641
|}
642
643
where <math display="inline">n_{i}</math> are the components (in the local system) of the normals to the element sides in the initial configuration <math display="inline">\Gamma _{M}^{0}</math>. The integration by parts of Eq.([[#eq-42|42]]) is typical in finite volume methods for computing second derivatives over volumes by line integrals of gradient terms <span id='citeF-28'></span>[[#cite-28|[28]],<span id='citeF-29'></span>[[#cite-29|29]]].
644
645
For the definition of the normal vector <math display="inline">\mathbf{t}_{3}</math>, the linear interpolation over the central element is used. In this case the tangent plane components are
646
647
<span id="eq-44.a"></span>
648
{| class="formulaSCP" style="width: 100%; text-align: left;" 
649
|-
650
| 
651
{| style="text-align: left; margin:auto;width: 100%;" 
652
|-
653
| style="text-align: center;" | <math>{\boldsymbol \varphi }_{^{\prime }\alpha } = \sum _{i=1}^{3} L_{i,\alpha }^M {\boldsymbol \varphi }_{i}\quad ,\quad \alpha=1,2 </math>
654
|}
655
| style="width: 5px;text-align: right;white-space: nowrap;" | (44.a)
656
|}
657
658
<span id="eq-44.b"></span>
659
{| class="formulaSCP" style="width: 100%; text-align: left;" 
660
|-
661
| 
662
{| style="text-align: left; margin:auto;width: 100%;" 
663
|-
664
| style="text-align: center;" | <math>\mathbf{t}_{3}=\frac{{\boldsymbol \varphi }_{\prime{1}}\times{\boldsymbol \varphi }_{\prime{2}}}{\left\vert {\boldsymbol \varphi }_{\prime{1}}\times{\boldsymbol \varphi }_{\prime{2}}\right\vert }=\lambda \;{\boldsymbol \varphi  }_{1}\times{\boldsymbol \varphi }_{2} </math>
665
|}
666
| style="width: 5px;text-align: right;white-space: nowrap;" | (44.b)
667
|}
668
669
From these expressions it is also possible to compute in the original configuration the element area <math display="inline">A^{0}_{M}</math>, the outer normals <math display="inline">\left( n_{1},n_{2}\right) ^{i}</math> at each side and the side lengths <math display="inline">l_{i}^{M}</math>. Equation ([[#eq-44.b|44.b]]) 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>.
670
671
The numerical evaluation of the line  integral in Eq.([[#eq-43|43]]) results in a sum over the integration points at the element boundary which are, in fact, the same points used for evaluating the gradients when computing the membrane strains. As one integration point is used over each side, it is not necessary to distinguish between sides (<math display="inline">i</math>) and integration points (<math display="inline">G_{i}</math>). In this way the curvatures can be computed by
672
673
<span id="eq-45"></span>
674
{| class="formulaSCP" style="width: 100%; text-align: left;" 
675
|-
676
| 
677
{| style="text-align: left; margin:auto;width: 100%;" 
678
|-
679
| style="text-align: center;" | <math>{\boldsymbol \kappa }=\frac{1}{A_{M}^{0}} \sum ^3_{i=1} l_i^M \left[ \begin{array}{cc}-n_{1} & 0\\ 0 & -n_{2}\\ -n_{2} & -n_{1}\end{array} \right] \left[ \begin{array}{c}\mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }1}\\ \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }2}\end{array} \right] d\Gamma </math>
680
|}
681
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
682
|}
683
684
Eq.([[#eq-45|45]]) is now expressed in terms  of the shape functions of the 3-noded triangle <math display="inline">L_i^M</math> (which coincide with the area coordinates <span id='citeF-4'></span>[[#cite-4|[4]]]). Noting the property of the area coordinates
685
686
<span id="eq-46"></span>
687
{| class="formulaSCP" style="width: 100%; text-align: left;" 
688
|-
689
| 
690
{| style="text-align: left; margin:auto;width: 100%;" 
691
|-
692
| style="text-align: center;" | <math>\nabla L_{i}^{M}=\left[ \begin{array}{c}L_{i,x}^{M}\\ L_{i,y}^{M}\end{array} \right] =-\frac{l_{i}^{M}}{2A_{M}}\left[ \begin{array}{c}n_{x}^{i}\\ n_{y}^{i}\end{array} \right]  </math>
693
|}
694
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
695
|}
696
697
the expression for the curvature can be expressed as
698
699
<span id="eq-47"></span>
700
{| class="formulaSCP" style="width: 100%; text-align: left;" 
701
|-
702
| 
703
{| style="text-align: left; margin:auto;width: 100%;" 
704
|-
705
| style="text-align: center;" | <math>{\boldsymbol \kappa }=2\sum _{i=1}^{3}\left[ \begin{array}{cc}L_{i,1}^M & 0\\ 0         & L_{i,2}^M \\ L_{i,2}^M & L_{i,1}^M \end{array} \right] \left[ \begin{array}{c}\mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }1}^{i}\\ \mathbf{t}_{3}\cdot{\boldsymbol \varphi }_{^{\prime }2}^{i}\end{array} \right]  </math>
706
|}
707
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
708
|}
709
710
The gradient <math display="inline">\mathbf{\boldsymbol \varphi  }_{\prime \alpha }^{i}</math>  is evaluated at each side <math display="inline">G_{i}</math> from the quadratic interpolation
711
712
<span id="eq-48"></span>
713
{| class="formulaSCP" style="width: 100%; text-align: left;" 
714
|-
715
| 
716
{| style="text-align: left; margin:auto;width: 100%;" 
717
|-
718
| style="text-align: center;" | <math>\left[ \begin{array}{c}{\boldsymbol \varphi }_{\prime{1}}^{i}\\ {\boldsymbol \varphi }_{\prime{2}}^{i}\end{array} \right] =\left[ \begin{array}{cccc}N_{1,1}^{i} & N_{2,1}^{i} & N_{3,1}^{i} & N_{i+3,1}^{i}\\ N_{1,2}^{i} & N_{2,2}^{i} & N_{3,2}^{i} & N_{i+3,2}^{i}\end{array} \right] \left[ \begin{array}{c}{\boldsymbol \varphi }_{1}\\ {\boldsymbol \varphi }_{2}\\ {\boldsymbol \varphi }_{3}\\ {\boldsymbol \varphi }_{i+3}\end{array} \right]  </math>
719
|}
720
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
721
|}
722
723
This is a basic difference with respect of the computation of the curvature field in the original Basic Shell Triangle (BST) where the gradient at the side mid-point is computed as the average value between the values at two adjacent elements <span id='citeF-20'></span><span id='citeF-23'></span><span id='citeF-26'></span><span id='citeF-27'></span>[[#cite-20|[20]],[[#cite-23|23]],[[#cite-26|26]],[[#cite-27|27]]].
724
725
Note again than at each side the gradients depend only on the positions of the three nodes of the central triangle and of an extra node (<math display="inline">i+3</math>), associated precisely to the side (<math display="inline">G_{i}</math>) where the gradient is computed.
726
727
Direction '''t'''<math display="inline">_{3}</math> in Eq.([[#eq-47|47]]) can be seen as a reference direction. If a different direction than that given by Eq.([[#eq-44.b|44.b]]) 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-44.b|44.b]]) for the definition of '''t'''<math display="inline">_{3}</math> as a function exclusively of the three nodes of the central triangle, instead of using the 6-node isoparametric interpolation.
728
729
The variation of the curvatures can be obtained as
730
731
<span id="eq-49"></span>
732
{| class="formulaSCP" style="width: 100%; text-align: left;" 
733
|-
734
| 
735
{| style="text-align: left; margin:auto;width: 100%;" 
736
|-
737
| style="text-align: center;" | <math>\delta{\boldsymbol \kappa }   =2\sum _{i=1}^{3}\left[ \begin{array}{cc}L_{i,1}^M & 0\\ 0         & L_{i,2}^M\\ L_{i,2}^M & L_{i,1}^M\end{array} \right] \left\{ \sum _{i=1}^{3}\left[ \begin{array}{c}N_{j,1}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}_{j})\\ N_{j,2}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}_{j}) \end{array} \right] +\left[ \begin{array}{c}N_{i+3,1}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}^{i+3})\\ N_{i+3,2}^{i}(\mathbf{t}_{3}\cdot \delta \mathbf{u}^{i+3}) \end{array} \right] \right\} -</math>
738
|-
739
| style="text-align: center;" | <math>   -\sum _{i=1}^{3}\left[ \begin{array}{c}(L_{i,1}^M\rho _{11}^{1}+L_{i,2}^M\rho _{11}^{2})\\ (L_{i,1}^M\rho _{22}^{1}+L_{i,2}^M\rho _{22}^{2})\\ (L_{i,1}^M\rho _{12}^{1}+L_{i,2}^M\rho _{12}^{2}) \end{array} \right] (\mathbf{t}_{3}\cdot \delta \mathbf{u}_{i})=\mathbf{B}_{b}\delta \mathbf{a}^{p}</math>
740
|}
741
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
742
|}
743
744
In Eq.([[#eq-49|49]])
745
746
<span id="eq-50"></span>
747
{| class="formulaSCP" style="width: 100%; text-align: left;" 
748
|-
749
| 
750
{| style="text-align: left; margin:auto;width: 100%;" 
751
|-
752
| style="text-align: center;" | <math>\mathbf{B}_{b}=[\mathbf{B}_{b_{1}},\mathbf{B}_{b_{2}},\cdots ,\mathbf{B}_{b_{6}}]</math>
753
|}
754
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
755
|}
756
757
Details of the derivation of the curvature matrix <math display="inline">\mathbf{B}_b</math> are given in <span id='citeF-26'></span><span id='citeF-27'></span>[[#cite-26|[26]],[[#cite-27|27]]].
758
759
===3.3 The EBST1 Element===
760
761
A simplified and yet very effective version of the EBST element can be obtained by using ''one point quadrature'' for the computation of all the element integrals. This element is termed EBST1. Note that this only affects the membrane stiffness matrices and it is equivalent to using a assumed constant membrane strain field defined by an average of the metric tensors computed at each side.
762
763
Numerical experiments have shown that both the EBST and the EBST1 elements are free of spurious energy modes.
764
765
==4 Boundary Conditions==
766
767
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 the bending  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">\boldsymbol \varphi _{^{\prime }n}^{0}</math> that does not change during the process.
768
769
The tangent plane at the boundary (mid-side point) is expressed in terms of two orthogonal unit vectors referred to a local-to-the-boundary Cartesian system (see Fig. [[#img-2|2]]) defined as
770
771
{| class="formulaSCP" style="width: 100%; text-align: left;" 
772
|-
773
| 
774
{| style="text-align: left; margin:auto;width: 100%;" 
775
|-
776
| style="text-align: center;" | <math>\left[\boldsymbol \varphi _{^{\prime }n}^{0},\;\bar{\boldsymbol \varphi }_{^{\prime }s}\right] </math>
777
|}
778
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
779
|}
780
781
where vector <math display="inline">\boldsymbol \varphi _{^{\prime }n}^{0}</math> is fixed during the process while direction <math display="inline">\bar{\boldsymbol \varphi }_{^{\prime }s}</math> emerges from the intersection of the symmetry plane with the plane defined by the central element (<math display="inline">M</math>). The plane (gradient) defined by the central element in the selected original convective Cartesian system (<math display="inline">\mathbf{t}_{1},\mathbf{t}_{2} </math>) is
782
783
{| class="formulaSCP" style="width: 100%; text-align: left;" 
784
|-
785
| 
786
{| style="text-align: left; margin:auto;width: 100%;" 
787
|-
788
| style="text-align: center;" | <math>\left[\boldsymbol \varphi _{^{\prime }1}^{M},\;\boldsymbol \varphi _{^{\prime  }2}^{M}\right] </math>
789
|}
790
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
791
|}
792
793
the intersection line (side <math display="inline">i</math>) of this plane with the plane of symmetry can be written in terms of the position of the nodes that define the side (<math display="inline">j </math> and <math display="inline">k</math>) and the original length of the side <math display="inline">l_{i}^{M}</math>, i.e.
794
795
{| class="formulaSCP" style="width: 100%; text-align: left;" 
796
|-
797
| 
798
{| style="text-align: left; margin:auto;width: 100%;" 
799
|-
800
| style="text-align: center;" | <math>\boldsymbol \varphi _{^{\prime }s}^{i}=\frac{1}{l_{i}^{M}}\left(\boldsymbol \varphi _{k}-\boldsymbol \varphi _{j}\right) </math>
801
|}
802
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
803
|}
804
805
That together with the outer normal to the side <math display="inline">\mathbf{n}^{i} =\left[n_{1},n_{2}\right]^{T}=\left[\mathbf{n\cdot t}_{1},\mathbf{n\cdot t}_{2}\right]^{T}</math> (resolved in the selected original convective Cartesian system) leads to
806
807
{| class="formulaSCP" style="width: 100%; text-align: left;" 
808
|-
809
| 
810
{| style="text-align: left; margin:auto;width: 100%;" 
811
|-
812
| style="text-align: center;" | <math>\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }1}^{iT} \\ \boldsymbol \varphi _{^{\prime }2}^{iT}\end{array}\right]=\left[ \begin{array}{cc}n_{1} & -n_{2} \\ n_{2} & n_{1}\end{array}\right]\left[ \begin{array}{c}\boldsymbol \varphi _{^{\prime }n}^{iT} \\ \boldsymbol \varphi _{^{\prime }s}^{iT}\end{array}\right] </math>
813
|}
814
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
815
|}
816
817
where, noting  that <math display="inline">\lambda </math> is the determinant of the gradient, the normal component of the gradient <math display="inline">\boldsymbol \varphi _{^{\prime }n}^{i}</math> can be approximated by
818
819
{| class="formulaSCP" style="width: 100%; text-align: left;" 
820
|-
821
| 
822
{| style="text-align: left; margin:auto;width: 100%;" 
823
|-
824
| style="text-align: center;" | <math>\boldsymbol \varphi _{^{\prime }n}^{i}=\frac{\boldsymbol \varphi _{^{\prime }n}^{0}}{\lambda |\boldsymbol \varphi _{^{\prime }s}^{i}|} </math>
825
|}
826
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
827
|}
828
829
<div id='img-2'></div>
830
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
831
|-
832
|[[Image:Draft_Samper_105745998-fig3.png|500px|Local Cartesian system for the treatment of symmetry boundary conditions]]
833
|- style="text-align: center; font-size: 75%;"
834
| colspan="1" | '''Figure 2:''' Local Cartesian system for the treatment of symmetry boundary conditions
835
|}
836
837
For a simple supported (hinged) side, the problem is not completely defined. The simplest choice is to neglect the contribution to the side rotations from the adjacent element missing in the patch in the evaluation of the curvatures via Eq.([[#eq-43|43]]) <span id='citeF-20'></span><span id='citeF-23'></span><span id='citeF-26'></span>[[#cite-20|[20]],[[#cite-23|23]],[[#cite-26|26]]]. This is equivalent to assume that the gradient at the side is equal to the gradient in the central element, i.e.
838
839
{| class="formulaSCP" style="width: 100%; text-align: left;" 
840
|-
841
| 
842
{| style="text-align: left; margin:auto;width: 100%;" 
843
|-
844
| style="text-align: center;" | <math>\left[\boldsymbol \varphi _{^{\prime }1}^{i},\;\boldsymbol \varphi _{^{\prime }2}^{i}\right]=\left[\boldsymbol \varphi _{^{\prime }1}^{M},\;\boldsymbol \varphi _{^{\prime }2}^{M}\right] </math>
845
|}
846
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
847
|}
848
849
More precise changes can be however introduced to account for the different natural boundary conditions. One may assume that the curvature normal to the side is zero, and consider a contribution of the missing side to introduce this constraint. As the change of curvature parallel to the side is also zero along the hinged side, this obviously leads to zero curvatures in both directions.
850
851
We note finally that for the membrane formulation of element EBST, the gradient at the mid-side point of the boundary is assumed equal to the gradient of the main triangle.
852
853
More details on the specification of the boundary conditions on the EBST element can be found in <span id='citeF-26'></span><span id='citeF-27'></span>[[#cite-26|[26]],[[#cite-27|27]]].
854
855
==5 Explicit Solution Scheme==
856
857
For simulations including large non-linearities, such as those occuring in sheet metal forming processes involving frictional contact conditions on complex geometries or large instabilities, convergence is difficult to achieve with implicit schemes. In those cases an explicit solution algorithm is typically most advantageous. This scheme provides the solution for dynamic problems and also for quasi-static problems if an adequate damping is chosen.
858
859
The dynamic equations of motion to solve are of the form
860
861
{| class="formulaSCP" style="width: 100%; text-align: left;" 
862
|-
863
| 
864
{| style="text-align: left; margin:auto;width: 100%;" 
865
|-
866
| style="text-align: center;" | <math>\mathbf{r}(\mathbf{u}) + \mathbf{D} \dot{\mathbf{u}} + \mathbf{M} \ddot{\mathbf{u}} = 0 </math>
867
|}
868
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
869
|}
870
871
where <math display="inline">\mathbf{M}</math> is the mass matrix, <math display="inline">\mathbf{D}</math> is the damping matrix and the dot means the time derivative. The solution is performed using the ''central difference method''. To make the method competitive a diagonal (lumped) <math display="inline">\mathbf{M}</math> matrix is typically used and <math display="inline">\mathbf{D}</math> is taken proportional to <math display="inline">\mathbf{M}</math>. As usual, mass lumping is performed by assigning one third of the triangular element mass to each node in the central element.
872
873
The explicit solution scheme can be summarized as follows. At each time step <math display="inline">n</math> where displacements have been computed:
874
875
<ol>
876
877
<li>Compute the internal forces <math display="inline">\mathbf{r}^{n}</math>. This follows the  steps described in Box [[#Box-1|1]]. </li>
878
879
<li>Compute the accelerations at time <math display="inline">t_{n}</math>
880
881
{| class="formulaSCP" style="width: 100%; text-align: left;" 
882
|-
883
| 
884
{| style="text-align: left; margin:auto;width: 100%;" 
885
|-
886
| style="text-align: center;" | <math>
887
888
\ddot{\mathbf{u}}^{n} = {\boldsymbol M}_d^{-1} [ \mathbf{r}^{n} - \mathbf{D} \dot{\mathbf{u}}^{n-1/2} ] </math>
889
|}
890
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
891
|}</li>
892
893
where <math display="inline">{\boldsymbol M}_d</math> is the diagonal (lumped) mass matrix.
894
895
<li>Compute the velocities at time <math display="inline">t_{n+1/2}</math>
896
897
{| class="formulaSCP" style="width: 100%; text-align: left;" 
898
|-
899
| 
900
{| style="text-align: left; margin:auto;width: 100%;" 
901
|-
902
| style="text-align: center;" | <math>
903
904
\dot{\mathbf{u}}^{n+1/2} = \dot{\mathbf{u}}^{n-1/2}+ \ddot{\mathbf{u}}^{n} \delta t </math>
905
|}
906
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
907
|}</li>
908
909
<li>Compute the displacements at  time <math display="inline">t_{n+1}</math>
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>
917
918
\mathbf{u}^{n+1} = \mathbf{u}^{n} +\dot{\mathbf{u}}^{n+1/2} \delta t </math>
919
|}
920
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
921
|}</li>
922
<li>Update the shell geometry </li>
923
<li>Check frictional contact conditions. </li>
924
925
</ol>
926
927
==6 Example 1. Cylindrical Panel under Impulse Loading==
928
929
The geometry of the cylinder and the material properties are shown in Fig. [[#img-3|3]]. A prescribed initial normal velocity of <math display="inline">v_{o}=-5650</math> in/sec is applied to the points in the region shown modelling the effect of the detonation of an explosive layer. The panel is assumed to be clamped along all the boundary. One half of the cylinder is discretized only due to symmetry conditions. Three different meshes of <math display="inline">6\times{12}</math>, <math display="inline">12\times{32}</math> and <math display="inline">18\times{48}</math>  triangles were used for the analysis. The deformed configurations for <math display="inline">time =1 msec</math> are shown for the three meshes in Fig. [[#img-3|3]].
930
931
The analysis was performed assuming an elastic-perfect plastic material behaviour (<math display="inline">\sigma _y = k_y</math> <math display="inline">k_y'=0</math>). A study of the convergence of the solution with the number of thickness layers showed again that four layers suffice to capture accurately the non linear material response <span id='citeF-25'></span>[[#cite-25|[25]]].
932
933
A comparison of the results obtained with the BST and EBST1 elements using the coarse mesh and the finer mesh is shown in Fig. [[#img-3|3]] where experimental results reported in <span id='citeF-32'></span>[[#cite-32|[32]]] have also been plotted for comparison purposes. Good agreement between the numerical and experimental results is obtained. Figs. [[#img-4|4]] show the time evolution of the vertical displacement of two reference points along the center line located at <math display="inline">y=6.28</math>in and <math display="inline">y=9.42</math>in, respectively. For the finer mesh results between both elements are almost identical. For the coarse mesh it can been seen  that the  BST element is more flexible than the  EBST1.
934
935
<span id="Box-1"></span>
936
937
{|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
938
939
|- style="border-top: 2px solid;"
940
941
<ol>
942
943
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |<li>Generate the actual configuration <math display="inline">\mathbf{\boldsymbol \varphi }^{n+1}=\mathbf{\boldsymbol \varphi }^{n}+\Delta \mathbf{u}^{n}</math>  </li>
944
<li>Compute the metric tensor <math display="inline">a_{\alpha \beta }^{n+1}\mathbf{ }</math>and the curvatures <math display="inline">\kappa _{\alpha \beta }^{n+1}</math>. Then at each layer <math display="inline">k</math> compute the (approximate) right Cauchy-Green tensor. From Eq.([[#eq-14|14]])
945
946
{| class="formulaSCP" style="width: 100%; text-align: left;" 
947
|-
948
| 
949
{| style="text-align: left; margin:auto;width: 100%;" 
950
|-
951
| style="text-align: center;" | <math> \mathbf{C}_{k}^{n+1}=\mathbf{a}^{n+1}+z_{k}{\boldsymbol \chi }^{n+1} </math>
952
|}
953
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
954
|}</li>
955
<li>Compute the total ([[#eq-21|21]]) and elastic ([[#eq-22|22]]) deformations at each layer <math display="inline">k</math>
956
957
{| class="formulaSCP" style="width: 100%; text-align: left;" 
958
|-
959
| 
960
{| style="text-align: left; margin:auto;width: 100%;" 
961
|-
962
| style="text-align: center;" | <math> {\boldsymbol \varepsilon }_{k}^{n+1}   = \frac{1}{2}\ln{\mathbf{C}_{k}^{n+1}} </math>
963
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
964
|-
965
| style="text-align: center;" | <math> \left[ {\boldsymbol \varepsilon }_{e}\right] _{k}^{n+1}   ={\boldsymbol \varepsilon  }_{k}^{n+1}-\left[ {\boldsymbol \varepsilon }_{p}\right] _{k}^{n} </math>
966
|}
967
|}</li>
968
<li>Compute the trial Hencky elastic stresses ([[#eq-23|23]]) at each layer <math display="inline">k</math>
969
970
{| class="formulaSCP" style="width: 100%; text-align: left;" 
971
|-
972
| 
973
{| style="text-align: left; margin:auto;width: 100%;" 
974
|-
975
| style="text-align: center;" | <math> \mathbf{T} _{k}^{n+1}=\mathbf{H}\left[ {\boldsymbol \varepsilon }_{e}\right] _{k}^{n+1} </math>
976
|}
977
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
978
|}</li>
979
<li>Check the plasticity condition and return to the plasticity surface. If necessary correct the plastic strains <math display="inline">\left[{\boldsymbol \varepsilon }_{p}\right] _{k}^{n+1}</math> at each layer  </li>
980
<li>Compute the second Piola-Kirchhoff stress vector <math display="inline">\boldsymbol \sigma _k^{n+1}</math> and the generalized stresses
981
982
{| class="formulaSCP" style="width: 100%; text-align: left;" 
983
|-
984
| 
985
{| style="text-align: left; margin:auto;width: 100%;" 
986
|-
987
| style="text-align: center;" | <math>\begin{array}{l} {\boldsymbol \sigma }^{n+1}_{m}  &  =\frac{h^{0}}{N_{L}}\sum _{k=1}^{N_{L}}\boldsymbol \sigma _{k}^{n+1} w_{k}\\[.25cm] {\boldsymbol \sigma }^{n+1}_{b}  &  =\frac{h^{0}}{N_{L}}\sum _{k=1}^{N_{L}}\boldsymbol \sigma _{k}^{n+1}z_{k} w_{k}\end{array}</math>
988
|}
989
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
990
|}</li>
991
992
Where <math display="inline"> w_{k}</math> is the weight of the through-the-thickness integration point. Recall that <math display="inline">z_{k}</math> is the current distance of the layer to the mid-surface and not the original distance. However, for small strain plasticity this distinction is not important.  This computation of stresses is  exact for an elastic problem.
993
<li>Compute the residual force vector from
994
995
<span id="eq-65"></span>
996
{| class="formulaSCP" style="width: 100%; text-align: left;" 
997
|-
998
| 
999
{| style="text-align: left; margin:auto;width: 100%;" 
1000
|-
1001
| style="text-align: center;" | <math> \mathbf{r}^e_i =\iint _A L_i {\boldsymbol t}\, dA - \iint _{A^\circ } ({\boldsymbol    B}_{m_i}^T {\boldsymbol \sigma }_m + {\boldsymbol B}_{b_i}^T {\boldsymbol \sigma }_b)dA  </math>
1002
|}
1003
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
1004
|}</li>
1005
1006
</ol>
1007
1008
|- style="border-bottom: 2px solid;"
1009
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;"|
1010
|}
1011
1012
<br/><div class="center" style="width: auto; margin-left: auto; margin-right: auto;">'''Box 1.''' Computation of the internal forces vector</div>
1013
1014
1015
<div id='img-3'></div>
1016
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1017
|-
1018
|[[Image:Draft_Samper_105745998-fig14.png|400px|]]
1019
|-
1020
|[[Image:Draft_Samper_105745998-fig15.png|600px|]]
1021
|- style="text-align: center; font-size: 75%;"
1022
| colspan="2" |'''Figure 3.''' Cylindrical panel under impulse loading. Geometry and material   properties. Deformed meshes for <math>time =1 msec</math>
1023
|}
1024
1025
<div id='img-4'></div>
1026
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1027
|-
1028
|[[Image:Draft_Samper_105745998-fig16.png|400px|Cylindrical panel under impulse loading. Time evolution of the   displacement of two points along the crown line. Upper lines y=6.28in. Lower lines y=9.42 in. Comparison of results obtained with BST and EBST1 elements (mesh 1: 6×12 elements and mesh 3: 18×48 elements) and experimental values]]
1029
|- style="text-align: center; font-size: 75%;"
1030
|'''Figure 4.''' Cylindrical panel under impulse loading. Time evolution of the   displacement of two points along the crown line. Upper lines <math>y=6.28</math>in. Lower lines <math>y=9.42</math> in. Comparison of results obtained with BST and EBST1 elements (mesh 1: <math>6\times{12}</math> elements and mesh 3: <math>18\times{48}</math> elements) and experimental values
1031
|-
1032
|}
1033
1034
1035
The numerical values of the vertical displacement at the two reference points obtained with the BST and EBST1  elements after a time of 0.4 ms using the <math display="inline">16\times{32}</math> mesh are compared in Table [[#table-2|2]]  with a numerical solution obtained by Stolarski ''et al.'' <span id='citeF-31'></span>[[#cite-31|[31]]] using a curved triangular shell element and the <math display="inline">16\times{32}</math> mesh. Experimental results reported in <span id='citeF-32'></span>[[#cite-32|[32]]] are also given for comparison. It is interesting to note the reasonable agreement of the results for <math display="inline">y=6.28</math>in. and the discrepancy of present and other published numerical solutions with the experimental value for <math display="inline">y=9.42</math>in.
1036
1037
1038
{|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
1039
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Cylindrical panel under impulse load. Comparison of vertical displacement values of two central points for <math>t=0.4</math> ms
1040
|- style="border-top: 2px solid;"
1041
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  
1042
| colspan='2' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Vertical displacement (in.)
1043
|- style="border-top: 2px solid;"
1044
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  element/mesh                
1045
| style="border-left: 2px solid;border-right: 2px solid;" | <math>y=6.28</math>in 
1046
| style="border-left: 2px solid;border-right: 2px solid;" | <math>y=9.42</math>in 
1047
|- style="border-top: 2px solid;"
1048
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  BST  (<math display="inline"> 6\times 12</math> el.)    
1049
| style="border-left: 2px solid;border-right: 2px solid;" | -1.310     
1050
| style="border-left: 2px solid;border-right: 2px solid;" | -0.679      
1051
|-
1052
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | BST  (<math display="inline">18\times 48</math> el.)    
1053
| style="border-left: 2px solid;border-right: 2px solid;" | -1.181     
1054
| style="border-left: 2px solid;border-right: 2px solid;" | -0.587      
1055
|-
1056
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | EBST1 (<math display="inline"> 6\times 12</math> el.)    
1057
| style="border-left: 2px solid;border-right: 2px solid;" | -1.147     
1058
| style="border-left: 2px solid;border-right: 2px solid;" | -0.575      
1059
|-
1060
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | EBST1 (<math display="inline">18\times 48</math> el.)    
1061
| style="border-left: 2px solid;border-right: 2px solid;" | -1.171     
1062
| style="border-left: 2px solid;border-right: 2px solid;" | -0.584      
1063
|-
1064
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Stolarski ''et al.'' <span id='citeF-31'></span>[[#cite-31|[31]]] 
1065
| style="border-left: 2px solid;border-right: 2px solid;" | -1.183     
1066
| style="border-left: 2px solid;border-right: 2px solid;" | -0.530      
1067
|- style="border-bottom: 2px solid;"
1068
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Experimental <span id='citeF-32'></span>[[#cite-32|[32]]] 
1069
| style="border-left: 2px solid;border-right: 2px solid;" | -1.280     
1070
| style="border-left: 2px solid;border-right: 2px solid;" | -0.700      
1071
1072
|}
1073
1074
<div id='img-5'></div>
1075
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1076
|-
1077
|[[Image:Draft_Samper_105745998-fig17.png|300px|Cylindrical panel under impulse loading. Final deformation (t=1 msec) of the panel at the cross section y=6.28 in. Comparison with experimental values]]
1078
|- style="text-align: center; font-size: 75%;"
1079
| colspan="1" | '''Figure 5:''' Cylindrical panel under impulse loading. Final deformation (<math>t=1 msec</math>) of the panel at the cross section <math>y=6.28 in</math>. Comparison with experimental values
1080
|}
1081
1082
<div id='img-6'></div>
1083
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1084
|-
1085
|[[Image:Draft_Samper_105745998-fig18.png|500px|Cylindrical panel under impulse loading. Final deformation (t=1 msec) of the panel at the crown line (x=0.00 in). Comparison with experimental values]]
1086
|- style="text-align: center; font-size: 75%;"
1087
| colspan="1" | '''Figure 6:''' Cylindrical panel under impulse loading. Final deformation (<math>t=1 msec</math>) of the panel at the crown line (<math>x=0.00 in</math>). Comparison with experimental values
1088
|}
1089
1090
The deformed shapes of the transverse section for <math display="inline">y=6.28</math>in. and the longitudinal section for <math display="inline">x=0</math> obtained with the both elements for the coarse and the fine meshes after 1ms. are compared with the experimental results in Figs. [[#img-5|5]] and [[#img-6|6]].  Excellent agreement is observed for the fine mesh for both elements.
1091
1092
==7 Application to Sheet Metal Forming Problems==
1093
1094
The features of tghe EBST1 element make it ideal for analysis of sheet metal stamping processes. A number of examples of simulations of practical problems of this kind are presented. Numerical results have been obtained with the sheet stamping simulation code STAMPACK where the EBST1 element has been implemented <span id='citeF-35'></span>[[#cite-35|[35]]].
1095
1096
===7.1 S-rail Sheet Stamping===
1097
1098
The next problem corresponds to one of the sheet stamping benchmark tests proposed in NUMISHEET'96 <span id='citeF-33'></span>[[#cite-33|[33]]].  The analysis comprises two parts, namely, simulation of the stamping of a S-rail sheet component and springback computations once the stamping tools are removed.  Figure [[#img-7|7]] shows the deformed sheet after springback.
1099
1100
<div id='img-7'></div>
1101
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1102
|-
1103
|[[Image:Draft_Samper_105745998-fig24.png|400px|Stamping of a S-rail. Final deformation of the sheet after springback obtained in the simulation. The triangular mesh of the deformed sheet is also shown]]
1104
|- style="text-align: center; font-size: 75%;"
1105
| colspan="1" | '''Figure 7:''' Stamping of a S-rail. Final deformation of the sheet after springback obtained in the simulation. The triangular mesh of the deformed sheet is also shown
1106
|}
1107
1108
The detailed geometry and material data can be found in the proceedings of the conference <span id='citeF-33'></span>[[#cite-33|[33]]] or in the web <span id='citeF-34'></span>[[#cite-34|[34]]]. The mesh used for the sheet has 6000 triangles and 3111 points (Fig. [[#img-7|7]]). The tools are treated as rigid bodies. The meshes used for the sheet and the tools are those provided by the  benchmark organizers. The material considered here is a mild steel (IF) with Young Modulus <math display="inline">E=2.06 GPa</math> and Poisson ratio <math display="inline">\nu=0.3</math>. Mises yield criterion was used for plasticity behaviour with non-linear isotropic hardening defined by <math display="inline">\sigma _y(e^p) = 545(0.13+e^p)^{0.267} [MPa]</math>. A uniform friction of 0.15 was used for all the tools. A low (10kN) blank holder force was considered in this simulation.
1109
1110
Figure [[#img-8|8]] compares the punch force during the stamping stage obtained with both BST and EBST1 elements for the simulation and experimental values. Also for reference the average values of the simulations presented in the conference are included. Explicit and implicit simulations are considered as different curves. There is a remarkable coincidence between the experimental values and the results obtained with both the BST and EBST1 elements.
1111
1112
<div id='img-8'></div>
1113
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1114
|-
1115
|[[Image:Draft_Samper_105745998-fig25.png|400px|Stamping of a S-rail. Punch force versus punch travel. Average of explicit and implicit results reported at the benchmark conference are also shown]]
1116
|- style="text-align: center; font-size: 75%;"
1117
| colspan="1" | '''Figure 8:''' Stamping of a S-rail. Punch force versus punch travel. Average of explicit and implicit results reported at the benchmark conference are also shown
1118
|}
1119
1120
<div id='img-9'></div>
1121
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1122
|-
1123
|[[Image:Draft_Samper_105745998-fig26.png|400px|Stamping of a S-rail. Z-coordinate along line B''&#8211;-G'' after springback. Average of explicit and implicit results reported at the benchmark conference are also shown]]
1124
|- style="text-align: center; font-size: 75%;"
1125
| colspan="1" | '''Figure 9:''' Stamping of a S-rail. Z-coordinate along line B''&#8211;-G'' after springback. Average of explicit and implicit results reported at the benchmark conference are also shown
1126
|}
1127
1128
Figure [[#img-9|9]] plots the <math display="inline">Z</math> coordinate along line B"&#8211;G" after springback. The top surface of the sheet does not remain plane due to some instabilities due to the low blank holder force used. Results obtained with the simulations compare very well with the experimental values.
1129
1130
===7.2 Stamping of Industrial Automotive Part===
1131
1132
Figure [[#img-10|10]]  shows the geometry of the lateral panel of a car and the mesh of 457760 EBST1 elements used for the computation. Results of the stamping simulation are shown in Fig. [[#img-11|11]]. Note that the outpus of the simulation have been translated into graphical plots indicating the quality of the stamping process and the risk of failure in the different zones of the panel. This helps designers to taking decissions on the adequacy of the stamping process and for introducing changes in the design of the stamping tools (dies, punch, blankholders, etc.) and the process parameters if needed.
1133
1134
<div id='img-10'></div>
1135
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1136
|-
1137
|[[Image:Draft_Samper_105745998-lateral-panel.png|351px|Lateral panel of an automotive. Finite element mesh of 457760   triangles used   for the simulation]]
1138
|- style="text-align: center; font-size: 75%;"
1139
| colspan="1" | '''Figure 10:''' Lateral panel of an automotive. Finite element mesh of 457760   triangles used   for the simulation
1140
|}
1141
1142
Figure [[#img-12|12]] shows the geometry mesh and results of the stamping of a front fender part of an automotive. The initial mesh had 121960 EBST1 elements. Adaptive mesh refinement was used along the simulation process leading to a final mesh of 389870 elements. Finally, Figs. [[#img-13|13]] and [[#img-14|14]] show the same type of information for the stamping of a car tail gate. The initial and final meshes (after adaptive mesh refinement) had 186528 and 489560 EBST1 elements, respectively. The simulation results are displayed in both problems with an “engineering insight” in order to help the design and manufacturing of the stamping tools and the definition of the stamping process as previously mentioned.
1143
1144
<div id='img-11'></div>
1145
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1146
|-
1147
|[[Image:Draft_Samper_105745998-f_zone_237.png|550px|]]
1148
|[[Image:Draft_Samper_105745998-s_zone_237.png|550px|]]
1149
|-
1150
|[[Image:Draft_Samper_105745998-mesh_det_2_237.png|300px|]]
1151
|[[Image:Draft_Samper_105745998-relth_det_237.png|300px|]]
1152
|- style="text-align: center; font-size: 75%;"
1153
| colspan="2" | '''Figure 11:''' Lateral panel of a car. Results of the stamping analysis 
1154
|}
1155
1156
<div id='img-12'></div>
1157
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1158
|-
1159
|[[Image:Draft_Samper_105745998-guardabarros_2.png|500px|Front fender. Results of the stamping analysis using an   initial mesh of   121960 EBST1 elements. The final adapted mesh had 389870 elements]]
1160
|- style="text-align: center; font-size: 75%;"
1161
| colspan="1" | '''Figure 12:''' Front fender. Results of the stamping analysis using an   initial mesh of   121960 EBST1 elements. The final adapted mesh had 389870 elements
1162
|}
1163
1164
<div id='img-13'></div>
1165
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1166
|-
1167
|[[Image:Draft_Samper_105745998-chapas-juntas.png|351px|Car tail gate. Geometry and final adapted mesh of 489560 EBST1 elements used for the  stamping simulation]]
1168
|- style="text-align: center; font-size: 75%;"
1169
| colspan="1" | '''Figure 13:''' Car tail gate. Geometry and final adapted mesh of 489560 EBST1 elements used for the  stamping simulation
1170
|}
1171
1172
<div id='img-14'></div>
1173
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1174
|-
1175
|[[Image:Draft_Samper_105745998-chapa_map_thick.png|351px|]]
1176
|-
1177
|[[Image:Draft_Samper_105745998-chapa_map_form.png|351px|Car tail gate. Map of relative thickness distribution and forming zones on the   stamped part]]
1178
|- style="text-align: center; font-size: 75%;"
1179
| colspan="2" | '''Figure 14:''' Car tail gate. Map of relative thickness distribution and forming zones on the   stamped part
1180
|}
1181
1182
==8 Concluding Remarks==
1183
1184
An enhanced rotation-free shell triangle (termed EBST) is obtained by using a quadratic interpolation of the geometry in terms of the six nodes belonging to the  four elements patch associated to each triangle.  This allows to computing an assumed constant curvature field and an assumed linear membrane strain field which improves the in-plane behaviour of the original element.  A simple and economic version of the   element using a single integration point has been presented.  The efficiency of the  rotation-free shell triangle has been demonstrated in examples of application including the analysis of a cylinder under impulse loading and practical sheet stamping problems.
1185
1186
The enhanced rotation-free basic shell triangle element with a single integration point (the EBST1 element) is an excellent candidate for solving practical  sheet metal stamping problems and other non linear shell problems in engineering involving complex geometry, dynamics, material and geometrical non linearities and frictional contact conditions.
1187
1188
==ACKNOWLEDGEMENTS==
1189
1190
The support of the company QUANTECH (www.quantech.es) providing the code STAMPACK <span id='citeF-35'></span>[[#cite-35|[35]]] is gratefully acknowledged.
1191
1192
This research was partially supported by project SEDUREC of the Consolider Programme of the Ministerio de Educación y Ciencia of Spain.
1193
1194
==References==
1195
1196
<div id="cite-1"></div>
1197
[[#citeF-1|[1]]]  Oñate E (1994) A review of some finite element families for thick and thin plate and shell analysis. Publication CIMNE N. 53, May
1198
1199
<div id="cite-2"></div>
1200
[[#citeF-2|[2]]]  Flores FG, Oñate E, Zárate F (1995) New assumed strain triangles for non-linear shell analysis. Computational Mechanics 17:107&#8211;114
1201
1202
<div id="cite-3"></div>
1203
[[#citeF-3|[3]]]  Argyris JH, Papadrakakis M, Apostolopoulou C, Koutsourelakis S (2000) The TRIC element. Theoretical and numerical investigation. Comput Meth Appl Mech Engrg 182:217&#8211;245
1204
1205
<div id="cite-4"></div>
1206
[[#citeF-4|[4]]]  Zienkiewicz OC, Taylor RL (2005) The finite element method. Solid Mechanics. Vol II, Elsevier
1207
1208
<div id="cite-5"></div>
1209
[[#citeF-5|[5]]]  Stolarski H, Belytschko T, Lee S-H  (1995) A review of shell finite elements and corotational theories. Computational Mechanics Advances vol. 2 (2), North-Holland
1210
1211
<div id="cite-6"></div>
1212
[[#citeF-6|[6]]]  Ramm E, Wall WA (2002) Shells in advanced computational environment. In V World Congress on Computational Mechanics, Eberhardsteiner J, Mang H, Rammerstorfer F (eds), Vienna, Austria, July 7&#8211;12. http://wccm.tuwien.ac.at.
1213
1214
<div id="cite-7"></div>
1215
[[#citeF-7|[7]]]  Bushnell D, Almroth BO (1971) Finite difference energy method for non linear shell analysis. J Computers and Structures 1:361
1216
1217
<div id="cite-8"></div>
1218
[[#citeF-8|[8]]] Timoshenko SP (1971)  Theory of Plates and Shells, McGraw Hill, New York
1219
1220
<div id="cite-9"></div>
1221
[[#citeF-9|[9]]] Ugural AC (1981) Stresses in  Plates and Shells, McGraw Hill, New York
1222
1223
<div id="cite-10"></div>
1224
[[#citeF-10|[10]]]  Nay RA, Utku S (1972) An alternative to the finite element method. Variational Methods Eng vol. 1
1225
1226
<div id="cite-11"></div>
1227
[[#citeF-11|[11]]]  Hampshire JK, Topping BHV, Chan HC (1992) Three node triangular elements with one degree of freedom per node. Engng Comput  9:49&#8211;62
1228
1229
<div id="cite-12"></div>
1230
[[#citeF-12|[12]]] Phaal R, Calladine CR (1992) A simple class of finite elements for plate and shell problems. I: Elements for beams and thin plates. Int J Num Meth Engng 35:955&#8211;977
1231
1232
<div id="cite-13"></div>
1233
[[#citeF-13|[13]]] Phaal R, Calladine CR (1992)  A simple class of finite elements for plate and shell problems. II: An element for thin shells with only translational degrees of freedom. Int J Num Meth Engng 35:979&#8211;996
1234
1235
<div id="cite-14"></div>
1236
[[#citeF-14|[14]]] Oñate E, Cervera M (1993) Derivation of thin plate bending elements with one degree of freedom per node. Engineering Computations 10:553&#8211;561
1237
1238
<div id="cite-15"></div>
1239
[[#citeF-15|[15]]] Brunet M, Sabourin F (1994) Prediction of necking and wrinkles with a simplified shell element in sheet forming. Int Conf of Metal Forming Simulation in Industry II:27&#8211;48, Kröplin B (ed)
1240
1241
<div id="cite-16"></div>
1242
[[#citeF-16|[16]]] Rio G, Tathi B, Laurent H (1994) A new efficient finite element model of shell with only three degrees of freedom per node. Applications to industrial deep drawing test. In Recent Developments in Sheet Metal Forming Technology, Barata Marques MJM (ed), 18th IDDRG Biennial Congress, Lisbon  
1243
1244
<div id="cite-17"></div>
1245
[[#citeF-17|[17]]]  Rojek J, Oñate E (1998)  Sheet springback analysis using a simple shell triangle with translational degrees of freedom only. Int J of Forming Processes 1(3):275&#8211;296  
1246
1247
<div id="cite-18"></div>
1248
[[#citeF-18|[18]]]  Rojek J, Oñate E, Postek E  (1998) Application of explicit FE codes to simulation of sheet and bulk forming processes. J of Materials Processing Technology 80-81:620&#8211;627  
1249
1250
<div id="cite-19"></div>
1251
[[#citeF-19|[19]]]  Jovicevic J, Oñate E (1999) Analysis of beams and shells using a rotation-free finite element-finite volume formulation, Monograph 43, CIMNE, Barcelona
1252
1253
<div id="cite-20"></div>
1254
[[#citeF-20|[20]]] Oñate E, Zárate F (2000) Rotation-free plate and shell triangles. Int J Num Meth Engng 47:557&#8211;603
1255
1256
<div id="cite-21"></div>
1257
[[#citeF-21|[21]]] Cirak F, Ortiz M  (2000) Subdivision surfaces: A new paradigm for thin-shell finite element analysis. Int J Num Meths in Engng 47:2039-2072
1258
1259
<div id="cite-22"></div>
1260
[[#citeF-22|[22]]] Cirak F, Ortiz M  (2001) Fully <math display="inline">C^{1}</math>-conforming subdivision elements for finite deformations thin-shell analysis. Int J Num Meths in Engng 51:813-833
1261
1262
<div id="cite-23"></div>
1263
[[#citeF-23|[23]]] Flores FG, Oñate E (2001) A basic thin shell triangle with only translational DOFs for large strain plasticity. Int J Num Meths in Engng 51:57&#8211;83
1264
1265
<div id="cite-24"></div>
1266
[[#citeF-24|[24]]]  Engel G, Garikipati K, Hughes TJR,  Larson MG, Mazzei L, Taylor RL (2002). Continuous/discontinuous finite element approximation of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity. Comput Methods Appl Mech Engrg 191:3669&#8211;3750
1267
1268
<div id="cite-25"></div>
1269
[[#citeF-25|[25]]]  Oñate E, Cendoya P, Miquel J (2002) Non linear explicit dynamic analysis of shells using the BST rotation-free triangle. Engineering Computations 19(6):662&#8211;706
1270
1271
<div id="cite-26"></div>
1272
[[#citeF-26|[26]]] Flores FG, Oñate E (2005)  Improvements in the membrane behaviour of the three node rotation-free BST shell triangle  using an assumed strain approach. Comput Meth Appl Mech Engng 194(6-8):907&#8211;932
1273
1274
<div id="cite-27"></div>
1275
[[#citeF-27|[27]]]  Oñate E, Flores FG (2005)  Advances in the formulation     of the rotation-free basic shell triangle. Comput Meth Appl Mech Engng      194(21&#8211;24):2406&#8211;2443
1276
1277
<div id="cite-28"></div>
1278
[[#citeF-28|[28]]]   Zienkiewicz OC, Oñate E (1991)  Finite Elements vs. Finite Volumes. Is there really a choice?. Nonlinear Computational Mechanics. State of the Art. Wriggers P, Wagner R (eds), Springer Verlag, Heidelberg
1279
1280
<div id="cite-29"></div>
1281
[[#citeF-29|[29]]] Oñate E, Cervera M,  Zienkiewicz OC (1994) A finite volume format for structural mechanics. Int J Num Meth Engng 37:181&#8211;201
1282
1283
<div id="cite-30"></div>
1284
[[#citeF-30|[30]]] Hill R (1948) A Theory of the Yielding and Plastic Flow of Anisotropic Metals. Proc Royal Society London A193:281
1285
1286
<div id="cite-31"></div>
1287
[[#citeF-31|[31]]]  Stolarski H, Belytschko T, Carpenter N (1984) A simple triangular curved shell element. Eng Comput 1:210&#8211;218
1288
1289
<div id="cite-32"></div>
1290
[[#citeF-32|[32]]]  Balmer HA, Witmer EA (1964) Theoretical experimental correlation of large dynamic and permanent deformation of impulsively loaded simple structures. Air force flight Dynamic Lab Rep FDQ-TDR-64-108, Wright-Patterson AFB, Ohio, USA
1291
1292
<div id="cite-33"></div>
1293
[[#citeF-33|[33]]] NUMISHEET'96 Third International Conference and Workshop on Numerical Simulation of 3D Sheet Forming Processes, Lee EH, Kinzel GL, Wagoner RH (eds), Dearbon-Michigan, USA, 1996
1294
1295
<div id="cite-34"></div>
1296
[[#citeF-34|[34]]]  <code>http://rclsgi.eng.ohio-state.edu/%Elee-j-k/numisheet96/</code>
1297
1298
<div id="cite-35"></div>
1299
[[#citeF-35|[35]]] STAMPACK. A General Finite Element System for Sheet Stamping and Forming Problems, Quantech ATZ, Barcelona, Spain, 2007. (www.quantech.es)
1300

Return to Onate et al 2007c.

Back to Top