MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
shell.cc
Go to the documentation of this file.
1 /*
2  * MBDyn (C) is a multibody analysis code.
3  * http://www.mbdyn.org
4  *
5  * Copyright (C) 2010-2017
6  *
7  * Marco Morandini <morandini@aero.polimi.it>
8  * Riccardo Vescovini <vescovini@aero.polimi.it>
9  *
10  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
11  * via La Masa, 34 - 20156 Milano, Italy
12  * http://www.aero.polimi.it
13  *
14  * Changing this copyright notice is forbidden.
15  *
16  * This program is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation (version 2 of the License).
19  *
20  *
21  * This program is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with this program; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  */
30 
31 /*
32  * Inspired by
33  * Wojciech Witkowski
34  * "4-Node combined shell element with semi-EAS-ANS strain interpolations
35  * in 6-parameter shell theories with drilling degrees of freedom"
36  * Comput Mech (2009) 43:307­319 DOI 10.1007/s00466-008-0307-x
37  */
38 
39 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
40 
41 #include "constltp_impl.h"
42 #include "tpldrive_impl.h"
43 #include "shell.h"
44 #include "mynewmem.h"
45 
46 
47 Shell::Shell(unsigned uLabel, const DofOwner* pDO, flag fOut)
48 : Elem(uLabel, fOut),
49 ElemGravityOwner(uLabel, fOut),
50 ElemWithDofs(uLabel, pDO, fOut),
51 InitialAssemblyElem(uLabel, fOut)
52 {
53  NO_OP;
54 }
55 
56 Shell::~Shell(void)
57 {
58  NO_OP;
59 }
60 
61 int
62 ReadShellConstLaw(MBDynParser& HP, Shell::fmh& pD, Shell::vh& PreStress)
63 {
64  ASSERT(pD.iGetNumRows() == 12);
65  ASSERT(pD.iGetNumCols() == 12);
66 
67  if (HP.IsKeyWord("diag")) {
68  pD.Reset();
69  for (unsigned ir = 1; ir <= 12; ir++) {
70  pD(ir, ir) = HP.GetReal();
71  }
72 
73  } else if (HP.IsKeyWord("sym")) {
74  for (unsigned ir = 1; ir <= 12; ir++) {
75  pD(ir, ir) = HP.GetReal();
76  for (unsigned ic = ir + 1; ic <= 12; ic++) {
77  doublereal d = HP.GetReal();
78  pD(ir, ic) = d;
79  pD(ic, ir) = d;
80  }
81  }
82 
83  } else if (HP.IsKeyWord("isotropic")) {
84  doublereal dE;
85  doublereal dnu;
86  doublereal dG;
87  doublereal dh;
88  doublereal das = 1.;
89  doublereal dat = .01;
90  bool bGot_E(false);
91  bool bGot_nu(false);
92  bool bGot_G(false);
93  bool bGot_h(false);
94  bool bGot_as(false);
95  bool bGot_at(false);
96 
97  while (HP.IsArg()) {
98  if (HP.IsKeyWord("E") || HP.IsKeyWord("Young" "modulus")) {
99  if (bGot_E) {
100  silent_cerr("Shell isotropic constitutive law: Young's modulus already provided at line " << HP.GetLineData() << std::endl);
101  return -1;
102  }
103  bGot_E = true;
104  dE = HP.GetReal();
105  if (dE <= 0.) {
106  silent_cerr("Shell isotropic constitutive law: invalid Young's modulus " << dE << " at line " << HP.GetLineData() << std::endl);
107  return -1;
108  }
109 
110  } else if (HP.IsKeyWord("nu") || HP.IsKeyWord("Poisson" "modulus")) {
111  if (bGot_nu) {
112  silent_cerr("Shell isotropic constitutive law: Poisson's modulus already provided at line " << HP.GetLineData() << std::endl);
113  return -1;
114  }
115  bGot_nu = true;
116  dnu = HP.GetReal();
117  if (dnu < 0. || dnu >= .5) {
118  silent_cerr("Shell isotropic constitutive law: invalid Poisson's modulus " << dnu << " at line " << HP.GetLineData() << std::endl);
119  return -1;
120  }
121 
122  } else if (HP.IsKeyWord("G") || HP.IsKeyWord("shear" "modulus")) {
123  if (bGot_G) {
124  silent_cerr("Shell isotropic constitutive law: shear modulus already provided at line " << HP.GetLineData() << std::endl);
125  return -1;
126  }
127  bGot_G = true;
128  dG = HP.GetReal();
129  if (dG <= 0.) {
130  silent_cerr("Shell isotropic constitutive law: invalid shear modulus " << dG << " at line " << HP.GetLineData() << std::endl);
131  return -1;
132  }
133 
134  } else if (HP.IsKeyWord("thickness")) {
135  if (bGot_h) {
136  silent_cerr("Shell isotropic constitutive law: thickness already provided at line " << HP.GetLineData() << std::endl);
137  return -1;
138  }
139  bGot_h = true;
140  dh = HP.GetReal();
141  if (dh <= 0.) {
142  silent_cerr("Shell isotropic constitutive law: invalid thickness " << dh << " at line " << HP.GetLineData() << std::endl);
143  return -1;
144  }
145 
146  } else if (HP.IsKeyWord("as") /* better name? */ ) {
147  if (bGot_as) {
148  silent_cerr("Shell isotropic constitutive law: as (?) already provided at line " << HP.GetLineData() << std::endl);
149  return -1;
150  }
151  bGot_as = true;
152  das = HP.GetReal();
153  if (das <= 0.) {
154  silent_cerr("Shell isotropic constitutive law: invalid as " << das << " at line " << HP.GetLineData() << std::endl);
155  return -1;
156  }
157 
158  } else if (HP.IsKeyWord("at") /* better name? */ ) {
159  if (bGot_at) {
160  silent_cerr("Shell isotropic constitutive law: at (?) already provided at line " << HP.GetLineData() << std::endl);
161  return -1;
162  }
163  bGot_at = true;
164  dat = HP.GetReal();
165  if (dat <= 0.) {
166  silent_cerr("Shell isotropic constitutive law: invalid at " << dat << " at line " << HP.GetLineData() << std::endl);
167  return -1;
168  }
169 
170  } else {
171  break;
172  }
173  }
174 
175  int got = 0;
176  if (bGot_E) got++;
177  if (bGot_nu) got++;
178  if (bGot_G) got++;
179 
180  if (got < 2) {
181  silent_cerr("Shell isotropic constitutive law: incomplete material data (need at least two among Young's modulus, Poisson's modulus, shear modulus)" << std::endl);
182  return -1;
183  }
184 
185  if (got > 2) {
186  if (std::abs(dE/(2.*(1 + dnu)) - dG) >
187  std::numeric_limits<doublereal>::epsilon())
188  {
189  silent_cerr("Shell isotropic constitutive law: inconsistent material data" << std::endl);
190  return -1;
191  }
192  } else {
193  if (!bGot_G) {
194  dG = dE/(2.*(1 + dnu));
195 
196  } else if (!bGot_E) {
197  dE = 2.*(1 + dnu)*dG;
198 
199  } else if (!bGot_nu) {
200  dnu = dE/(2.*dG) - 1.;
201  if (dnu <= 0. || dnu >= .5) {
202  silent_cerr("Shell isotropic constitutive law: inconsistent material data (computed Poisson's modulus: " << dnu << ")" << std::endl);
203  return -1;
204  }
205  }
206  }
207 
208  if (!bGot_h) {
209  silent_cerr("Shell isotropic constitutive law: shell thickness missing" << std::endl);
210  return -1;
211  }
212 
213  doublereal C = dE/(1. - dnu*dnu)*dh;
214  doublereal D = C*dh*dh/12;
215  doublereal G = dG*dh;
216  doublereal F = G*dh*dh/12;
217 
218  pD.Reset();
219 
220  pD(1, 1) = C;
221  pD(1, 5) = C*dnu;
222  pD(2, 2) = 2*G;
223  pD(3, 3) = G*das;
224  pD(4, 4) = 2*G;
225  pD(5, 1) = C*dnu;
226  pD(5, 5) = C;
227  pD(6, 6) = G*das;
228 
229  pD(7, 7) = 2*F;
230  pD(8, 8) = D;
231  pD(8, 10) = -D*dnu;
232  pD(9, 9) = F*dat;
233  pD(10, 10) = D;
234  pD(10, 8) = -D*dnu;
235  pD(11, 11) = 2*F;
236  pD(12, 12) = F*dat;
237 
238  } else if (HP.IsKeyWord("plane" "stress" "orthotropic")) {
239 /*
240 Eshell =
241 [ h/(1-nu_lt^2/E_l*E_t)*E_l, 0, 0, 0, h/(1-nu_lt^2/E_l*E_t)*nu_lt*E_t, 0, 0, 0, 0, 0, 0, 0]
242 [ 0, 2*h*G, 0, 2*h*G, 0, 0, 0, 0, 0, 0, 0, 0]
243 [ 0, 0, h*G*as, 0, 0, 0, 0, 0, 0, 0, 0, 0]
244 [ 0, 2*h*G, 0, 2*h*G, 0, 0, 0, 0, 0, 0, 0, 0]
245 [ h/(1-nu_lt^2/E_l*E_t)*nu_lt*E_t, 0, 0, 0, h/(1-nu_lt^2/E_l*E_t)*E_t, 0, 0, 0, 0, 0, 0, 0]
246 [ 0, 0, 0, 0, 0, h*G*as, 0, 0, 0, 0, 0, 0]
247 [ 0, 0, 0, 0, 0, 0, 1/12*h^3/(1-nu_lt^2/E_l*E_t)*E_l, 0, 0, 0, 1/12*h^3/(1-nu_lt^2/E_l*E_t)*nu_lt*E_t, 0]
248 [ 0, 0, 0, 0, 0, 0, 0, 1/6*h^3*G, 0, 1/6*h^3*G, 0, 0]
249 [ 0, 0, 0, 0, 0, 0, 0, 0, 1/12*h^3*G*dat, 0, 0, 0]
250 [ 0, 0, 0, 0, 0, 0, 0, 1/6*h^3*G, 0, 1/6*h^3*G, 0, 0]
251 [ 0, 0, 0, 0, 0, 0, 1/12*h^3/(1-nu_lt^2/E_l*E_t)*nu_lt*E_t, 0, 0, 0, 1/12*h^3/(1-nu_lt^2/E_l*E_t)*E_t, 0]
252 [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/12*h^3*G*dat]
253 
254 */
255 
256  /* TODO:
257  - check singularity in shear coefficients
258  - add (optional) re-orientation?
259  */
260 
261  doublereal dE_l;
262  doublereal dE_t;
263  doublereal dnu_lt;
264  doublereal dnu_tl;
265  doublereal dG;
266  doublereal dBeta;
267  doublereal dh;
268  doublereal das = 1.;
269  doublereal dat = .01;
270  bool bGot_E_l(false);
271  bool bGot_E_t(false);
272  bool bGot_nu_lt(false);
273  bool bGot_nu_tl(false);
274  bool bGot_G(false);
275  bool bGot_Beta(false);
276  bool bGot_h(false);
277  bool bGot_as(false);
278  bool bGot_at(false);
279 
280  while (HP.IsArg()) {
281  if (HP.IsKeyWord("E_l") || HP.IsKeyWord("longitudinal" "Young" "modulus")) {
282  if (bGot_E_l) {
283  silent_cerr("Shell plane stress orthotropic constitutive law: longitudinal Young's modulus already provided at line " << HP.GetLineData() << std::endl);
284  return -1;
285  }
286  bGot_E_l = true;
287  dE_l = HP.GetReal();
288  if (dE_l <= 0.) {
289  silent_cerr("Shell plane stress orthotropic constitutive law: invalid longitudinal Young's modulus " << dE_l << " at line " << HP.GetLineData() << std::endl);
290  return -1;
291  }
292 
293  } else if (HP.IsKeyWord("E_t") || HP.IsKeyWord("transverse" "Young" "modulus")) {
294  if (bGot_E_t) {
295  silent_cerr("Shell plane stress orthotropic constitutive law: transverse Young's modulus already provided at line " << HP.GetLineData() << std::endl);
296  return -1;
297  }
298  bGot_E_t = true;
299  dE_t = HP.GetReal();
300  if (dE_t <= 0.) {
301  silent_cerr("Shell plane stress orthotropic constitutive law: invalid transverse Young's modulus " << dE_t << " at line " << HP.GetLineData() << std::endl);
302  return -1;
303  }
304 
305  } else if (HP.IsKeyWord("nu_lt") || HP.IsKeyWord("longitudinal" "transverse" "Poisson" "modulus")) {
306  if (bGot_nu_lt) {
307  silent_cerr("Shell plane stress orthotropic constitutive law: longitudinal transverse Poisson's modulus already provided at line " << HP.GetLineData() << std::endl);
308  return -1;
309  }
310  bGot_nu_lt = true;
311  dnu_lt = HP.GetReal();
312 #if 0
313  if (dnu_lt <= 0. || dnu_lt >= .5) {
314  silent_cerr("Shell plane stress orthotropic constitutive law: invalid longitudinal transverse Poisson's modulus " << dnu_lt << " at line " << HP.GetLineData() << std::endl);
315  return -1;
316  }
317 #endif
318 
319  } else if (HP.IsKeyWord("nu_tl") || HP.IsKeyWord("transverse" "longitudinal" "Poisson" "modulus")) {
320  if (bGot_nu_lt) {
321  silent_cerr("Shell plane stress orthotropic constitutive law: transverse longitudinal Poisson's modulus already provided at line " << HP.GetLineData() << std::endl);
322  return -1;
323  }
324  bGot_nu_tl = true;
325  dnu_tl = HP.GetReal();
326 #if 0
327  if (dnu_tl <= 0. || dnu_tl >= .5) {
328  silent_cerr("Shell plane stress orthotropic constitutive law: invalid transverse longitudinal Poisson's modulus " << dnu_lt << " at line " << HP.GetLineData() << std::endl);
329  return -1;
330  }
331 #endif
332 
333  } else if (HP.IsKeyWord("G") || HP.IsKeyWord("shear" "modulus")) {
334  if (bGot_G) {
335  silent_cerr("Shell plane stress orthotropic constitutive law: shear modulus already provided at line " << HP.GetLineData() << std::endl);
336  return -1;
337  }
338  bGot_G = true;
339  dG = HP.GetReal();
340  if (dG <= 0.) {
341  silent_cerr("Shell plane stress orthotropic constitutive law: invalid shear modulus " << dG << " at line " << HP.GetLineData() << std::endl);
342  return -1;
343  }
344 
345  } else if (HP.IsKeyWord("beta")) {
346  if (bGot_Beta) {
347  silent_cerr("Shell plane stress orthotropic constitutive law: fiber angle \"beta\" already provided at line " << HP.GetLineData() << std::endl);
348  return -1;
349  }
350  bGot_Beta = true;
351  dBeta = HP.GetReal();
352 
353  // temporary
354  silent_cerr("Shell plane stress orthotropic constitutive law: fiber angle \"beta\" not supported yet" << std::endl);
355  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
356 
357  } else if (HP.IsKeyWord("thickness")) {
358  if (bGot_h) {
359  silent_cerr("Shell plane stress orthotropic constitutive law: thickness already provided at line " << HP.GetLineData() << std::endl);
360  return -1;
361  }
362  bGot_h = true;
363  dh = HP.GetReal();
364  if (dh <= 0.) {
365  silent_cerr("Shell plane stress orthotropic constitutive law: invalid thickness " << dh << " at line " << HP.GetLineData() << std::endl);
366  return -1;
367  }
368 
369  } else if (HP.IsKeyWord("as") /* better name? */ ) {
370  if (bGot_as) {
371  silent_cerr("Shell plane stress orthotropic constitutive law: as (?) already provided at line " << HP.GetLineData() << std::endl);
372  return -1;
373  }
374  bGot_as = true;
375  das = HP.GetReal();
376  if (das <= 0.) {
377  silent_cerr("Shell plane stress orthotropic constitutive law: invalid as " << das << " at line " << HP.GetLineData() << std::endl);
378  return -1;
379  }
380 
381  } else if (HP.IsKeyWord("at") /* better name? */ ) {
382  if (bGot_at) {
383  silent_cerr("Shell plane stress orthotropic constitutive law: at (?) already provided at line " << HP.GetLineData() << std::endl);
384  return -1;
385  }
386  bGot_at = true;
387  dat = HP.GetReal();
388  if (dat <= 0.) {
389  silent_cerr("Shell plane stress orthotropic constitutive law: invalid at " << dat << " at line " << HP.GetLineData() << std::endl);
390  return -1;
391  }
392 
393  } else {
394  break;
395  }
396  }
397 
398  if (!bGot_E_l) {
399  silent_cerr("Shell plane stress orthotropic constitutive law: longitudinal Young's modulus missing at line " << HP.GetLineData() << std::endl);
400  return -1;
401  }
402 
403  if (!bGot_E_t) {
404  silent_cerr("Shell plane stress orthotropic constitutive law: transverse Young's modulus missing at line " << HP.GetLineData() << std::endl);
405  return -1;
406  }
407 
408  if (!bGot_nu_lt && !bGot_nu_tl) {
409  silent_cerr("Shell plane stress orthotropic constitutive law: Poisson's modulus missing at line " << HP.GetLineData() << std::endl);
410  return -1;
411 
412  } else if (bGot_nu_lt && bGot_nu_tl) {
413  if (std::abs(dnu_tl*dE_l - dnu_lt*dE_t) >
414  std::numeric_limits<doublereal>::epsilon())
415  {
416  silent_cerr("Shell plane stress orthotropic constitutive law: inconsistent material data" << std::endl);
417  return -1;
418  }
419 
420  } else if (bGot_nu_tl) {
421  dnu_lt = dnu_tl*dE_l/dE_t;
422  }
423 
424  if (!bGot_G) {
425  silent_cerr("Shell plane stress orthotropic constitutive law: shear modulus missing at line " << HP.GetLineData() << std::endl);
426  return -1;
427  }
428 
429  doublereal C = dh/(1. - dnu_lt*dnu_lt/dE_l*dE_t);
430  doublereal G = dh*dG;
431  doublereal D = 1./12.*dh*dh*dh/(1. - dnu_lt*dnu_lt/dE_l*dE_t);
432  doublereal F = 1./12.*dh*dh*dh*dG;
433 
434  pD(1, 1) = C*dE_l;
435  pD(1, 5) = C*dnu_lt*dE_t;
436  pD(2, 2) = 2.*G;
437  // pD(2, 4) = 2.*G;
438  pD(3, 3) = G*das;
439  // pD(4, 2) = 2.*G;
440  pD(4, 4) = 2.*G;
441  pD(5, 1) = C*dnu_lt*dE_t;
442  pD(5, 5) = C*dE_t;
443  pD(6, 6) = G*das;
444 
445 #if 0
446  pD(7, 7) = D*dE_l;
447  pD(7, 11) = D*dnu_lt*dE_t;
448  pD(8, 8) = 2.*F;
449  // pD(8, 10) = 2*F;
450  pD(9, 9) = F*dat;
451  // pD(10, 8) = 2*F;
452  pD(10, 10) = 2*F;
453  pD(11, 7) = D*dnu_lt*dE_t;
454  pD(11, 11) = D*dE_t;
455  pD(12, 12) = F*dat;
456 #endif
457 
458  pD(7, 7) = 2.*F;
459  pD(8, 8) = D*dE_l;
460  pD(8, 10) = -D*dnu_lt*dE_t;
461  // pD(8, 10) = 2*F;
462  pD(9, 9) = F*dat;
463  // pD(10, 8) = 2*F;
464  pD(10, 8) = -D*dnu_lt*dE_t;
465  pD(10, 10) = D*dE_t;
466  pD(11, 11) = 2*F;
467  pD(12, 12) = F*dat;
468 
469  } else {
470  if (HP.IsKeyWord("matr")) {
471  // tolerate it
472  }
473 
474  for (unsigned ir = 1; ir <= 12; ir++) {
475  for (unsigned ic = 1; ic <= 12; ic++) {
476  pD(ir, ic) = HP.GetReal();
477  }
478  }
479  }
480 
481  if (HP.IsKeyWord("prestress")) {
482  ASSERT(PreStress.iGetSize() == 12);
483  for (unsigned ir = 1; ir <= 12; ir++) {
484  PreStress(ir) = HP.GetReal();
485  }
486  }
487 
488  return 0;
489 }
490 
constltp_impl.h
flag
long int flag
Definition: mbdyn.h:43
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
tpldrive_impl.h
Shell::Shell
Shell(unsigned uLabel, const DofOwner *pDO, flag fOut)
Definition: shell.cc:47
NO_OP
#define NO_OP
Definition: myassert.h:74
MyVectorHandler::iGetSize
virtual integer iGetSize(void) const
Definition: vh.h:255
ReadShellConstLaw
int ReadShellConstLaw(MBDynParser &HP, Shell::fmh &pD, Shell::vh &PreStress)
Definition: shell.cc:62
Shell::~Shell
virtual ~Shell(void)
Definition: shell.cc:56
HighParser::IsKeyWord
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
FullMatrixHandler::iGetNumCols
virtual integer iGetNumCols(void) const
Definition: fullmh.h:229
HighParser::IsArg
virtual bool IsArg(void)
Definition: parser.cc:807
Elem
Definition: elem.h:75
doublereal
double doublereal
Definition: colamd.c:52
IncludeParser::GetLineData
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
FullMatrixHandler::iGetNumRows
virtual integer iGetNumRows(void) const
Definition: fullmh.h:225
FullMatrixHandler::Reset
void Reset(void)
Definition: fullmh.cc:44
HighParser::GetReal
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056