1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Alessandro Tasora, Radu Serban
13 // =============================================================================
14
15 #include "chrono/fea/ChMaterialShellReissner.h"
16
17 namespace chrono {
18 namespace fea {
19
20
ComputeStiffnessMatrix(ChMatrixRef mC,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)21 void ChElasticityReissner::ComputeStiffnessMatrix(ChMatrixRef mC,
22 const ChVector<>& eps_u,
23 const ChVector<>& eps_v,
24 const ChVector<>& kur_u,
25 const ChVector<>& kur_v,
26 const double z_inf,
27 const double z_sup,
28 const double angle) {
29 assert(mC.rows() == 12);
30 assert(mC.cols() == 12);
31
32 mC.setZero();
33
34 ChVectorN<double, 12> strain_0;
35 strain_0.segment(0, 3) = eps_u.eigen();
36 strain_0.segment(3, 3) = eps_v.eigen();
37 strain_0.segment(6, 3) = kur_u.eigen();
38 strain_0.segment(9, 3) = kur_v.eigen();
39
40 ChVector<> nu, nv, mu, mv;
41 this->ComputeStress(nu, nv, mu, mv, eps_u, eps_v, kur_u, kur_v, z_inf, z_sup, angle);
42
43 ChVectorN<double, 12> stress_0;
44 stress_0.segment(0,3) = nu.eigen();
45 stress_0.segment(3,3) = nv.eigen();
46 stress_0.segment(6,3) = mu.eigen();
47 stress_0.segment(9,3) = mv.eigen();
48
49 double delta = 1e-9;
50 for (int i = 0; i < 12; ++i) {
51 strain_0(i, 0) += delta;
52 ChVector<> deps_u(strain_0.segment(0, 3));
53 ChVector<> deps_v(strain_0.segment(3, 3));
54 ChVector<> dkur_u(strain_0.segment(6, 3));
55 ChVector<> dkur_v(strain_0.segment(9, 3));
56 this->ComputeStress(nu, nv, mu, mv, deps_u, deps_v, dkur_u, dkur_v, z_inf, z_sup, angle);
57 ChVectorN<double, 12> stress_1;
58 stress_1.segment(0,3) = nu.eigen();
59 stress_1.segment(3,3) = nv.eigen();
60 stress_1.segment(6,3) = mu.eigen();
61 stress_1.segment(9,3) = mv.eigen();
62 ChVectorN<double, 12> stress_d = (1. / delta) * (stress_1 - stress_0);
63 mC.block(0, i, 12, 1) = stress_d;
64 strain_0(i, 0) -= delta;
65 }
66 }
67
68
69
70 //--------------------------------------------------------------
71
ChElasticityReissnerIsothropic(double E,double nu,double alpha,double beta)72 ChElasticityReissnerIsothropic::ChElasticityReissnerIsothropic( double E,
73 double nu,
74 double alpha,
75 double beta) {
76 m_E = E;
77 m_nu = nu;
78 m_alpha = alpha;
79 m_beta = beta;
80 }
81
ComputeStress(ChVector<> & n_u,ChVector<> & n_v,ChVector<> & m_u,ChVector<> & m_v,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)82 void ChElasticityReissnerIsothropic::ComputeStress( ChVector<>& n_u,
83 ChVector<>& n_v,
84 ChVector<>& m_u,
85 ChVector<>& m_v,
86 const ChVector<>& eps_u,
87 const ChVector<>& eps_v,
88 const ChVector<>& kur_u,
89 const ChVector<>& kur_v,
90 const double z_inf,
91 const double z_sup,
92 const double angle) {
93 if (z_inf == -z_sup) {
94 // simplified computation for centered layer
95 double h = z_sup - z_inf;
96 double G = m_E / (2. * (1. + m_nu));
97 double C = m_E * h / (1. - m_nu * m_nu);
98 double D = C * h * h / 12.;
99 double F = G * h * h * h / 12.;
100
101 n_u.x() = eps_u.x() * C + eps_v.y() * m_nu * C;
102 n_u.y() = eps_u.y() * 2 * G * h;
103 n_u.z() = eps_u.z() * m_alpha * G * h;
104 n_v.x() = eps_v.x() * 2 * G * h;
105 n_v.y() = eps_v.y() * C + eps_u.x() * m_nu * C;
106 n_v.z() = eps_v.z() * m_alpha * G * h;
107
108 m_u.x() = kur_u.x() * 2 * F;
109 m_u.y() = kur_u.y() * D + kur_v.x() * (-m_nu * D);
110 m_u.z() = kur_u.z() * m_beta * F;
111 m_v.x() = kur_v.x() * D + kur_u.y() * (-m_nu * D);
112 m_v.y() = kur_v.y() * 2 * F;
113 m_v.z() = kur_v.z() * m_beta * F;
114 } else {
115 // throw ChException("ComputeStiffnessMatrix not yet implemented for non-centered layers");
116 double G = m_E / (2. * (1. + m_nu));
117 double Q11 = m_E / (1. - m_nu * m_nu);
118 double Q22 = Q11;
119 double Q12 = m_nu * Q11;
120 double Q33 = 2 * G;
121 double Q44 = 2 * G;
122 double Qss = m_alpha * 2 * G;
123 double Qdd = m_beta * 2 * G;
124 double h1 = z_sup - z_inf;
125 double h2 = 0.5 * (pow(z_sup, 2) - pow(z_inf, 2));
126 double h3 = (1. / 3.) * (pow(z_sup, 3) - pow(z_inf, 3));
127
128 n_u.x() = h1 * (eps_u.x() * Q11 + eps_v.y() * Q12) + h2 * (kur_u.y() * Q11 + kur_v.x() * Q12);
129 n_u.y() = h1 * (eps_u.y() * Q33) + h2 * (kur_u.x() * Q33);
130 n_u.z() = h1 * (eps_u.z() * Qss);
131 n_v.x() = h1 * (eps_v.x() * Q44) + h2 * (kur_v.y() * Q44);
132 n_v.y() = h1 * (eps_u.x() * Q12 + eps_v.y() * Q22) + h2 * (kur_u.y() * Q12 + kur_v.x() * Q22);
133 n_v.z() = h1 * (eps_v.z() * Qss);
134
135 m_u.x() = h2 * (eps_u.y() * Q33) + h3 * (kur_u.x() * Q33);
136 m_u.y() = h2 * (eps_u.x() * Q11 + eps_v.y() * Q12) + h3 * (kur_u.y() * Q11 + kur_v.x() * Q12);
137 m_u.z() = h3 * (eps_u.z() * Qdd);
138 m_v.x() = h2 * (eps_u.x() * Q12 + eps_v.y() * Q22) + h3 * (kur_u.y() * Q12 + kur_v.x() * Q22);
139 m_v.y() = h2 * (eps_v.x() * Q44) + h3 * (kur_v.y() * Q44);
140 m_v.z() = h3 * (eps_v.z() * Qdd);
141 }
142 }
143
ComputeStiffnessMatrix(ChMatrixRef mC,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)144 void ChElasticityReissnerIsothropic::ComputeStiffnessMatrix(ChMatrixRef mC,
145 const ChVector<>& eps_u,
146 const ChVector<>& eps_v,
147 const ChVector<>& kur_u,
148 const ChVector<>& kur_v,
149 const double z_inf,
150 const double z_sup,
151 const double angle) {
152 assert(mC.rows() == 12);
153 assert(mC.cols() == 12);
154
155 mC.setZero();
156
157 if (z_inf == -z_sup) {
158 // simplified computation for centered layer
159 double h = z_sup - z_inf;
160 double G = m_E / (2. * (1. + m_nu));
161 double C = m_E * h / (1. - m_nu * m_nu);
162 double D = C * h * h / 12.;
163 double F = G * h * h * h / 12.;
164 mC(0, 0) = C;
165 mC(0, 4) = m_nu * C;
166 mC(4, 0) = m_nu * C;
167 mC(1, 1) = 2. * G * h;
168 mC(2, 2) = m_alpha * G * h;
169 mC(3, 3) = 2. * G * h;
170 mC(4, 4) = C;
171 mC(5, 5) = m_alpha * G * h;
172 mC(6, 6) = 2. * F;
173 mC(7, 7) = D;
174 mC(7, 9) = -m_nu * D;
175 mC(9, 7) = -m_nu * D;
176 mC(8, 8) = m_beta * F;
177 mC(9, 9) = D;
178 mC(10, 10) = 2. * F;
179 mC(11, 11) = m_beta * F;
180 } else {
181 // throw ChException("ComputeStiffnessMatrix not yet implemented for non-centered layers");
182 double G = m_E / (2. * (1. + m_nu));
183 double Q11 = m_E / (1. - m_nu * m_nu);
184 double Q22 = Q11;
185 double Q12 = m_nu * Q11;
186 double Q33 = 2 * G;
187 double Q44 = 2 * G;
188 double Qss = m_alpha * 2 * G;
189 double Qdd = m_beta * 2 * G;
190 double h1 = z_sup - z_inf;
191 double h2 = 0.5 * (pow(z_sup, 2) - pow(z_inf, 2));
192 double h3 = (1. / 3.) * (pow(z_sup, 3) - pow(z_inf, 3));
193
194 mC(0, 0) = h1 * Q11;
195 mC(0, 4) = h1 * Q12;
196 mC(0, 7) = h2 * Q11;
197 mC(0, 9) = h2 * Q12;
198 mC(1, 1) = h1 * Q33;
199 mC(1, 6) = h2 * Q33;
200 mC(2, 2) = h1 * Qss;
201 mC(3, 3) = h1 * Q44;
202 mC(3, 10) = h2 * Q44;
203 mC(4, 0) = h1 * Q12;
204 mC(4, 4) = h1 * Q22;
205 mC(4, 7) = h2 * Q12;
206 mC(4, 9) = h2 * Q22;
207 mC(5, 5) = h1 * Qss;
208 mC(6, 1) = h2 * Q33;
209 mC(6, 6) = h3 * Q33;
210 mC(7, 0) = h2 * Q11;
211 mC(7, 4) = h2 * Q12;
212 mC(7, 7) = h3 * Q11;
213 mC(7, 9) = h3 * Q12;
214 mC(8, 8) = h3 * Qdd;
215 mC(9, 0) = h2 * Q12;
216 mC(9, 4) = h2 * Q22;
217 mC(9, 7) = h3 * Q12;
218 mC(9, 9) = h3 * Q22;
219 mC(10, 3) = h2 * Q44;
220 mC(10, 10) = h3 * Q44;
221 mC(11, 11) = h3 * Qdd;
222 }
223 }
224
225 //--------------------------------------------------------------
226
227 /// Construct an orthotropic material
ChElasticityReissnerOrthotropic(double m_E_x,double m_E_y,double m_nu_xy,double m_G_xy,double m_G_xz,double m_G_yz,double m_alpha,double m_beta)228 ChElasticityReissnerOrthotropic::ChElasticityReissnerOrthotropic( double m_E_x,
229 double m_E_y,
230 double m_nu_xy,
231 double m_G_xy,
232 double m_G_xz,
233 double m_G_yz,
234 double m_alpha,
235 double m_beta) {
236 E_x = m_E_x;
237 E_y = m_E_y;
238 nu_xy = m_nu_xy;
239 G_xy = m_G_xy;
240 G_xz = m_G_xz;
241 G_yz = m_G_yz;
242 alpha = m_alpha;
243 beta = m_beta;
244 }
245
ChElasticityReissnerOrthotropic(double m_E,double m_nu,double m_alpha,double m_beta)246 ChElasticityReissnerOrthotropic::ChElasticityReissnerOrthotropic( double m_E,
247 double m_nu,
248 double m_alpha,
249 double m_beta) {
250
251 double m_G = m_E / (2. * (1. + m_nu)); // default value of G for special subcase of isotropic constructor
252 this->E_x = m_E;
253 this->E_y = m_E;
254 this->nu_xy = m_nu;
255 this->G_xy = m_G;
256 this->G_xz = m_G;
257 this->G_yz = m_G;
258 this->alpha = m_alpha;
259 this->beta = m_beta;
260 }
261
ComputeStress(ChVector<> & n_u,ChVector<> & n_v,ChVector<> & m_u,ChVector<> & m_v,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)262 void ChElasticityReissnerOrthotropic::ComputeStress(ChVector<>& n_u,
263 ChVector<>& n_v,
264 ChVector<>& m_u,
265 ChVector<>& m_v,
266 const ChVector<>& eps_u,
267 const ChVector<>& eps_v,
268 const ChVector<>& kur_u,
269 const ChVector<>& kur_v,
270 const double z_inf,
271 const double z_sup,
272 const double angle) {
273 // Since it is a linear material, just compute S by using the
274 // constitutive matrix, as S = C*eps, where S={n_u, n_v, m_u, m_v}
275 ChMatrixNM<double, 12, 12> mC;
276 this->ComputeStiffnessMatrix(mC, eps_u, eps_v, kur_u, kur_v, z_inf, z_sup, angle);
277
278 ChVectorN<double, 12> eps;
279 eps.segment(0,3) = eps_u.eigen();
280 eps.segment(3,3) = eps_v.eigen();
281 eps.segment(6,3) = kur_u.eigen();
282 eps.segment(9,3) = kur_v.eigen();
283
284 ChVectorN<double, 12> Sigma = mC * eps;
285 n_u = Sigma.segment(0, 3);
286 n_v = Sigma.segment(3, 3);
287 m_u = Sigma.segment(6, 3);
288 m_v = Sigma.segment(9, 3);
289 }
290
ComputeStiffnessMatrix(ChMatrixRef mC,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)291 void ChElasticityReissnerOrthotropic::ComputeStiffnessMatrix(ChMatrixRef mC,
292 const ChVector<>& eps_u,
293 const ChVector<>& eps_v,
294 const ChVector<>& kur_u,
295 const ChVector<>& kur_v,
296 const double z_inf,
297 const double z_sup,
298 const double angle) {
299 assert(mC.rows() == 12);
300 assert(mC.cols() == 12);
301
302 mC.setZero();
303
304 // Compute Qm_local for inplane stresses as in sigma_local = Qm_local * eps_local
305 double nu_yx = this->Get_nu_yx(); // follows xy as it must be nu_yx*E_x = nu_xy*E_y
306 ChMatrixNM<double, 4, 4> Qm_local;
307 Qm_local.setZero();
308 Qm_local(0, 0) = E_x / (1. - nu_xy * nu_yx);
309 Qm_local(0, 1) = (nu_xy * E_y) / (1. - nu_xy * nu_yx);
310 Qm_local(1, 0) = (nu_yx * E_x) / (1. - nu_xy * nu_yx);
311 Qm_local(1, 1) = E_y / (1. - nu_xy * nu_yx);
312 Qm_local(2, 2) = 2. * G_xy;
313 Qm_local(3, 3) = 2. * G_xy;
314 // Compute Qs_local for transverse shear stresses as in sigma_local = Qs_local * eps_local
315 ChMatrixNM<double, 2, 2> Qs_local;
316 Qs_local(0, 0) = 2. * G_xz;
317 Qs_local(1, 1) = 2. * G_yz;
318 Qs_local(0, 1) = 0;
319 Qs_local(1, 0) = 0;
320
321 // Rotate Qm_local into Qm, as Qm = Tm'*Qm_local*Tm
322 double Co = cos(angle);
323 double Si = sin(angle);
324 double CC = Co * Co;
325 double SS = Si * Si;
326 double SC = Si * Co;
327 ChMatrixNM<double, 4, 4> Tm;
328 Tm(0, 0) = CC;
329 Tm(0, 1) = SS;
330 Tm(0, 2) = SC;
331 Tm(0, 3) = SC;
332 Tm(1, 0) = SS;
333 Tm(1, 1) = CC;
334 Tm(1, 2) = -SC;
335 Tm(1, 3) = -SC;
336 Tm(2, 0) = -SC;
337 Tm(2, 1) = SC;
338 Tm(2, 2) = CC;
339 Tm(2, 3) = -SS;
340 Tm(3, 0) = -SC;
341 Tm(3, 1) = SC;
342 Tm(3, 2) = -SS;
343 Tm(3, 3) = CC;
344
345 ChMatrixNM<double, 4, 4> Qm = Tm.transpose() * Qm_local * Tm;
346
347 // Rotate Qs_local into Qs, as Qs = Ts'*Qs_local*Ts
348 ChMatrixNM<double, 2, 2> Ts;
349 Ts(0, 0) = Co;
350 Ts(0, 1) = -Si;
351 Ts(1, 0) = Si;
352 Ts(1, 1) = Co;
353
354 ChMatrixNM<double, 2, 2> Qs = Ts.transpose() * Qs_local * Ts;
355
356 // Fill the 12x12 constitutive matrix
357 // upper left part
358 double h = z_sup - z_inf;
359 mC(0, 0) = h * Qm(0, 0);
360 mC(0, 1) = h * Qm(0, 2);
361 mC(0, 3) = h * Qm(0, 3);
362 mC(0, 4) = h * Qm(0, 1);
363 mC(1, 0) = h * Qm(2, 0);
364 mC(1, 1) = h * Qm(2, 2);
365 mC(1, 3) = h * Qm(2, 3);
366 mC(1, 4) = h * Qm(2, 1);
367 mC(3, 0) = h * Qm(3, 0);
368 mC(3, 1) = h * Qm(3, 2);
369 mC(3, 3) = h * Qm(3, 3);
370 mC(3, 4) = h * Qm(3, 1);
371 mC(4, 0) = h * Qm(1, 0);
372 mC(4, 1) = h * Qm(1, 2);
373 mC(4, 3) = h * Qm(1, 3);
374 mC(4, 4) = h * Qm(1, 1);
375 mC(2, 2) = h * alpha * Qs(0, 0);
376 mC(2, 5) = h * alpha * Qs(0, 1);
377 mC(5, 2) = h * alpha * Qs(1, 0);
378 mC(5, 5) = h * alpha * Qs(1, 1);
379 // lower right part
380 double H = (1. / 3.) * (pow(z_sup, 3) - pow(z_inf, 3));
381 mC(6, 6) = H * Qm(2, 2);
382 mC(6, 7) = H * Qm(2, 0);
383 mC(6, 9) = H * Qm(2, 1);
384 mC(6, 10) = H * Qm(2, 3);
385 mC(7, 6) = H * Qm(0, 2);
386 mC(7, 7) = H * Qm(0, 0);
387 mC(7, 9) = H * Qm(0, 1);
388 mC(7, 10) = H * Qm(0, 3);
389 mC(9, 6) = H * Qm(1, 2);
390 mC(9, 7) = H * Qm(1, 0);
391 mC(9, 9) = H * Qm(1, 1);
392 mC(9, 10) = H * Qm(1, 3);
393 mC(10, 6) = H * Qm(3, 2);
394 mC(10, 7) = H * Qm(3, 0);
395 mC(10, 9) = H * Qm(3, 1);
396 mC(10, 10) = H * Qm(3, 3);
397 mC(8, 8) = H * beta * Qs(0, 0);
398 mC(8, 11) = H * beta * Qs(0, 1);
399 mC(11, 8) = H * beta * Qs(1, 0);
400 mC(11, 11) = H * beta * Qs(1, 1);
401 // upper right part and lower right part
402 double hh = (0.5) * (pow(z_sup, 2) - pow(z_inf, 2));
403 mC(0, 6) = hh * Qm(0, 2);
404 mC(0, 7) = hh * Qm(0, 0);
405 mC(0, 9) = hh * Qm(0, 1);
406 mC(0, 10) = hh * Qm(0, 3);
407 mC(1, 6) = hh * Qm(2, 2);
408 mC(1, 7) = hh * Qm(2, 0);
409 mC(1, 9) = hh * Qm(2, 1);
410 mC(1, 10) = hh * Qm(2, 3);
411 mC(3, 6) = hh * Qm(3, 2);
412 mC(3, 7) = hh * Qm(3, 0);
413 mC(3, 9) = hh * Qm(3, 1);
414 mC(3, 10) = hh * Qm(3, 3);
415 mC(4, 6) = hh * Qm(1, 2);
416 mC(4, 7) = hh * Qm(1, 0);
417 mC(4, 9) = hh * Qm(1, 1);
418 mC(4, 10) = hh * Qm(1, 3);
419 mC(6, 0) = hh * Qm(2, 0);
420 mC(6, 1) = hh * Qm(2, 2);
421 mC(6, 3) = hh * Qm(2, 3);
422 mC(6, 4) = hh * Qm(2, 1);
423 mC(7, 0) = hh * Qm(0, 0);
424 mC(7, 1) = hh * Qm(0, 2);
425 mC(7, 3) = hh * Qm(0, 3);
426 mC(7, 4) = hh * Qm(0, 1);
427 mC(9, 0) = hh * Qm(1, 0);
428 mC(9, 1) = hh * Qm(1, 2);
429 mC(9, 3) = hh * Qm(1, 3);
430 mC(9, 4) = hh * Qm(1, 1);
431 mC(10, 0) = hh * Qm(3, 0);
432 mC(10, 1) = hh * Qm(3, 2);
433 mC(10, 3) = hh * Qm(3, 3);
434 mC(10, 4) = hh * Qm(3, 1);
435 }
436
437
438 //----------------------------------------------------------------
439
ChElasticityReissnerGeneric()440 ChElasticityReissnerGeneric::ChElasticityReissnerGeneric() {
441 mE.setIdentity();
442 }
443
ComputeStress(ChVector<> & n_u,ChVector<> & n_v,ChVector<> & m_u,ChVector<> & m_v,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)444 void ChElasticityReissnerGeneric::ComputeStress(
445 ChVector<>& n_u, ///< forces along \e u direction (per unit length)
446 ChVector<>& n_v, ///< forces along \e v direction (per unit length)
447 ChVector<>& m_u, ///< torques along \e u direction (per unit length)
448 ChVector<>& m_v, ///< torques along \e v direction (per unit length)
449 const ChVector<>& eps_u, ///< strains along \e u direction
450 const ChVector<>& eps_v, ///< strains along \e v direction
451 const ChVector<>& kur_u, ///< curvature along \e u direction
452 const ChVector<>& kur_v, ///< curvature along \e v direction
453 const double z_inf, ///< layer lower z value (along thickness coord)
454 const double z_sup, ///< layer upper z value (along thickness coord)
455 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic
456 ) {
457 ChVectorN<double, 12> mstrain;
458 ChVectorN<double, 12> mstress;
459 mstrain.segment(0 , 3) = eps_u.eigen();
460 mstrain.segment(3 , 3) = eps_v.eigen();
461 mstrain.segment(6 , 3) = kur_u.eigen();
462 mstrain.segment(9 , 3) = kur_v.eigen();
463 mstress = this->mE * mstrain;
464 n_u = mstress.segment(0, 3);
465 n_v = mstress.segment(3, 3);
466 m_u = mstress.segment(6, 3);
467 m_v = mstress.segment(9, 3);
468 }
469
ComputeStiffnessMatrix(ChMatrixRef mC,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle)470 void ChElasticityReissnerGeneric::ComputeStiffnessMatrix(
471 ChMatrixRef mC, ///< tangent matrix
472 const ChVector<>& eps_u, ///< strains along \e u direction
473 const ChVector<>& eps_v, ///< strains along \e v direction
474 const ChVector<>& kur_u, ///< curvature along \e u direction
475 const ChVector<>& kur_v, ///< curvature along \e v direction
476 const double z_inf, ///< layer lower z value (along thickness coord)
477 const double z_sup, ///< layer upper z value (along thickness coord)
478 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic
479 ) {
480 mC = this->mE;
481 }
482
483
484
485 //----------------------------------------------------------------
486
487
ChPlasticityReissner()488 ChPlasticityReissner::ChPlasticityReissner() : section(nullptr), nr_yeld_tolerance(1e-7), nr_yeld_maxiters(5) {}
489
ComputeStiffnessMatrixElastoplastic(ChMatrixRef K,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const ChShellReissnerInternalData & data,const double z_inf,const double z_sup,const double angle)490 void ChPlasticityReissner::ComputeStiffnessMatrixElastoplastic(
491 ChMatrixRef K, ///< 12x12 material elastoplastic stiffness matrix values here
492 const ChVector<>& eps_u, ///< strains along \e u direction
493 const ChVector<>& eps_v, ///< strains along \e v direction
494 const ChVector<>& kur_u, ///< curvature along \e u direction
495 const ChVector<>& kur_v, ///< curvature along \e v direction
496 const ChShellReissnerInternalData& data, ///< updated material internal variables, at this point including {p_strain_e, p_strain_k, p_strain_acc}
497 const double z_inf, ///< layer lower z value (along thickness coord)
498 const double z_sup, ///< layer upper z value (along thickness coord)
499 const double angle ///< layer angle respect to x (if needed)
500 ) {
501 ChVector<> n_u;
502 ChVector<> n_v;
503 ChVector<> m_u;
504 ChVector<> m_v;
505
506 std::vector<std::unique_ptr<ChShellReissnerInternalData>> a_plastic_data;
507 this->CreatePlasticityData(1, a_plastic_data);
508 std::vector<std::unique_ptr<ChShellReissnerInternalData>> b_plastic_data;
509 this->CreatePlasticityData(1, b_plastic_data);
510
511 bool in_plastic = ComputeStressWithReturnMapping(n_u, n_v, m_u, m_v,
512 *a_plastic_data[0], eps_u, eps_v, kur_u, kur_v, data, z_inf, z_sup, angle
513 );
514
515
516 if (!in_plastic) {
517
518 // if no return mapping is needed at this strain state, just use elastic matrix:
519
520 return this->section->GetElasticity()->ComputeStiffnessMatrix(K,eps_u, eps_v, kur_u, kur_v, z_inf, z_sup, angle);
521
522 } else {
523
524 // if return mapping is needed at this strain state, compute the elastoplastic stiffness by brute force BDF
525
526 ChVectorN<double, 12> strain_0;
527 strain_0.segment(0, 3) = eps_u.eigen();
528 strain_0.segment(3, 3) = eps_v.eigen();
529 strain_0.segment(6, 3) = kur_u.eigen();
530 strain_0.segment(9, 3) = kur_v.eigen();
531
532 ChVectorN<double, 12> stress_0;
533 stress_0.segment(0,3) = n_u.eigen();
534 stress_0.segment(3,3) = n_v.eigen();
535 stress_0.segment(6,3) = m_u.eigen();
536 stress_0.segment(9,3) = m_v.eigen();
537
538 double delta = 1e-6;
539 double invdelta = 1.0 / delta;
540 for (int i = 0; i < 12; ++i) {
541 strain_0(i, 0) += delta;
542 ChVector<> deps_u(strain_0.segment(0, 3));
543 ChVector<> deps_v(strain_0.segment(3, 3));
544 ChVector<> dkur_u(strain_0.segment(6, 3));
545 ChVector<> dkur_v(strain_0.segment(9, 3));
546 this->ComputeStressWithReturnMapping(n_u, n_v, m_u, m_v, *b_plastic_data[0], deps_u, deps_v, dkur_u, dkur_v, data, z_inf, z_sup, angle);
547 ChVectorN<double, 12> stress_1;
548 stress_1.segment(0,3) = n_u.eigen();
549 stress_1.segment(3,3) = n_v.eigen();
550 stress_1.segment(6,3) = m_u.eigen();
551 stress_1.segment(9,3) = m_v.eigen();
552 ChVectorN<double, 12> stress_d = invdelta * (stress_1 - stress_0);
553 K.block(0, i, 12, 1) = stress_d;
554 strain_0(i, 0) -= delta;
555 }
556
557
558 }
559 }
560
CreatePlasticityData(int numpoints,std::vector<std::unique_ptr<ChShellReissnerInternalData>> & plastic_data)561 void ChPlasticityReissner::CreatePlasticityData(
562 int numpoints,
563 std::vector<std::unique_ptr<ChShellReissnerInternalData>>& plastic_data) {
564 plastic_data.resize(numpoints);
565 for (int i = 0; i < numpoints; ++i) {
566 plastic_data[i] = std::unique_ptr<ChShellReissnerInternalData>(new ChShellReissnerInternalData());
567 }
568 }
569
570
571 //-----------------------------------------------------------------------
572
573
ComputeDampingMatrix(ChMatrixRef R,const ChVector<> & deps_u,const ChVector<> & deps_v,const ChVector<> & dkur_u,const ChVector<> & dkur_v,const double z_inf,const double z_sup,const double angle)574 void ChDampingReissner::ComputeDampingMatrix(
575 ChMatrixRef R, // 12x12 material damping matrix values here
576 const ChVector<>& deps_u, // time derivative of strains along \e u direction
577 const ChVector<>& deps_v, // time derivative of strains along \e v direction
578 const ChVector<>& dkur_u, // time derivative of curvature along \e u direction
579 const ChVector<>& dkur_v, // time derivative of curvature along \e v direction
580 const double z_inf, // layer lower z value (along thickness coord)
581 const double z_sup, // layer upper z value (along thickness coord)
582 const double angle // layer angle respect to x (if needed) -not used in this, isotropic+
583 ) {
584 assert(R.rows() == 12);
585 assert(R.cols() == 12);
586
587 R.setZero();
588
589 ChVectorN<double, 12> dstrain_0;
590 dstrain_0.segment(0, 3) = deps_u.eigen();
591 dstrain_0.segment(3, 3) = deps_v.eigen();
592 dstrain_0.segment(6, 3) = dkur_u.eigen();
593 dstrain_0.segment(9, 3) = dkur_v.eigen();
594
595 ChVector<> nu, nv, mu, mv;
596 this->ComputeStress(nu, nv, mu, mv, deps_u, deps_v, dkur_u, dkur_v, z_inf, z_sup, angle);
597
598 ChVectorN<double, 12> stress_0;
599 stress_0.segment(0,3) = nu.eigen();
600 stress_0.segment(3,3) = nv.eigen();
601 stress_0.segment(6,3) = mu.eigen();
602 stress_0.segment(9,3) = mv.eigen();
603
604 double delta = 1e-9;
605 for (int i = 0; i < 12; ++i) {
606 dstrain_0(i, 0) += delta;
607 ChVector<> ddeps_u(dstrain_0.segment(0, 3));
608 ChVector<> ddeps_v(dstrain_0.segment(3, 3));
609 ChVector<> ddkur_u(dstrain_0.segment(6, 3));
610 ChVector<> ddkur_v(dstrain_0.segment(9, 3));
611 this->ComputeStress(nu, nv, mu, mv, ddeps_u, ddeps_v, ddkur_u, ddkur_v, z_inf, z_sup, angle);
612 ChVectorN<double, 12> stress_1;
613 stress_1.segment(0,3) = nu.eigen();
614 stress_1.segment(3,3) = nv.eigen();
615 stress_1.segment(6,3) = mu.eigen();
616 stress_1.segment(9,3) = mv.eigen();
617 ChVectorN<double, 12> stress_d = (1. / delta) * (stress_1 - stress_0);
618 R.block(0, i, 12, 1) = stress_d;
619 dstrain_0(i, 0) -= delta;
620 }
621 }
622
623
624
625 // -----------------------------------------------------------------------------
626
627
ChDampingReissnerRayleigh(std::shared_ptr<ChElasticityReissner> melasticity,const double & mbeta)628 ChDampingReissnerRayleigh::ChDampingReissnerRayleigh(
629 std::shared_ptr<ChElasticityReissner> melasticity,
630 const double& mbeta) {
631 this->beta = mbeta;
632 this->section_elasticity = melasticity;
633 this->updated = false;
634 }
635
636
637
ComputeStress(ChVector<> & n_u,ChVector<> & n_v,ChVector<> & m_u,ChVector<> & m_v,const ChVector<> & deps_u,const ChVector<> & deps_v,const ChVector<> & dkur_u,const ChVector<> & dkur_v,const double z_inf,const double z_sup,const double angle)638 void ChDampingReissnerRayleigh::ComputeStress(
639 ChVector<>& n_u, ///< forces along \e u direction (per unit length)
640 ChVector<>& n_v, ///< forces along \e v direction (per unit length)
641 ChVector<>& m_u, ///< torques along \e u direction (per unit length)
642 ChVector<>& m_v, ///< torques along \e v direction (per unit length)
643 const ChVector<>& deps_u, ///< time derivative of strains along \e u direction
644 const ChVector<>& deps_v, ///< time derivative of strains along \e v direction
645 const ChVector<>& dkur_u, ///< time derivative of curvature along \e u direction
646 const ChVector<>& dkur_v, ///< time derivative of curvature along \e v direction
647 const double z_inf, ///< layer lower z value (along thickness coord)
648 const double z_sup, ///< layer upper z value (along thickness coord)
649 const double angle ///< layer angle respect to x (if needed)
650 ) {
651 if (!this->updated && this->section_elasticity->section) {
652 this->section_elasticity->ComputeStiffnessMatrix(this->E_const, VNULL, VNULL, VNULL, VNULL, z_inf, z_sup, angle);
653 this->updated = true;
654 }
655 ChVectorN<double, 12> mdstrain;
656 ChVectorN<double, 12> mstress;
657 mdstrain.segment(0 , 3) = deps_u.eigen();
658 mdstrain.segment(3 , 3) = deps_v.eigen();
659 mdstrain.segment(6 , 3) = dkur_u.eigen();
660 mdstrain.segment(9 , 3) = dkur_v.eigen();
661 mstress = this->beta * this->E_const * mdstrain;
662 n_u = mstress.segment(0, 3);
663 n_v = mstress.segment(3, 3);
664 m_u = mstress.segment(6, 3);
665 m_v = mstress.segment(9, 3);
666 }
667
ComputeDampingMatrix(ChMatrixRef R,const ChVector<> & deps_u,const ChVector<> & deps_v,const ChVector<> & dkur_u,const ChVector<> & dkur_v,const double z_inf,const double z_sup,const double angle)668 void ChDampingReissnerRayleigh::ComputeDampingMatrix( ChMatrixRef R, ///< 12x12 material damping matrix values here
669 const ChVector<>& deps_u, ///< time derivative of strains along \e u direction
670 const ChVector<>& deps_v, ///< time derivative of strains along \e v direction
671 const ChVector<>& dkur_u, ///< time derivative of curvature along \e u direction
672 const ChVector<>& dkur_v, ///< time derivative of curvature along \e v direction
673 const double z_inf, ///< layer lower z value (along thickness coord)
674 const double z_sup, ///< layer upper z value (along thickness coord)
675 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic
676 ) {
677 R = this->beta * this->E_const;
678 }
679
680
681
682 // -----------------------------------------------------------------------------
683
ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity)684 ChMaterialShellReissner::ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity) {
685 this->SetElasticity(melasticity);
686 }
687
ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity,std::shared_ptr<ChPlasticityReissner> mplasticity)688 ChMaterialShellReissner::ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity,
689 std::shared_ptr<ChPlasticityReissner> mplasticity) {
690 this->SetElasticity(melasticity);
691 this->SetPlasticity(mplasticity);
692 }
693
ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity,std::shared_ptr<ChPlasticityReissner> mplasticity,std::shared_ptr<ChDampingReissner> mdamping)694 ChMaterialShellReissner::ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity,
695 std::shared_ptr<ChPlasticityReissner> mplasticity,
696 std::shared_ptr<ChDampingReissner> mdamping) {
697 this->SetElasticity(melasticity);
698
699 if (mplasticity)
700 this->SetPlasticity(mplasticity);
701
702 if (mdamping)
703 this->SetDamping(mdamping);
704 }
705
ComputeStress(ChVector<> & n_u,ChVector<> & n_v,ChVector<> & m_u,ChVector<> & m_v,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle,ChShellReissnerInternalData * mdata_new,const ChShellReissnerInternalData * mdata)706 void ChMaterialShellReissner::ComputeStress(ChVector<>& n_u,
707 ChVector<>& n_v,
708 ChVector<>& m_u,
709 ChVector<>& m_v,
710 const ChVector<>& eps_u,
711 const ChVector<>& eps_v,
712 const ChVector<>& kur_u,
713 const ChVector<>& kur_v,
714 const double z_inf,
715 const double z_sup,
716 const double angle,
717 ChShellReissnerInternalData* mdata_new,
718 const ChShellReissnerInternalData* mdata) {
719 if (!plasticity || !mdata || !mdata_new)
720 this->elasticity->ComputeStress(n_u, n_v, m_u, m_v, eps_u, eps_v, kur_u, kur_v, z_inf, z_sup, angle);
721 else {
722 this->plasticity->ComputeStressWithReturnMapping(n_u, n_v, m_u, m_v, *mdata_new, eps_u, eps_v, kur_u, kur_v, *mdata, z_inf, z_sup, angle);
723 }
724 }
725
ComputeStiffnessMatrix(ChMatrixRef K,const ChVector<> & eps_u,const ChVector<> & eps_v,const ChVector<> & kur_u,const ChVector<> & kur_v,const double z_inf,const double z_sup,const double angle,const ChShellReissnerInternalData * mdata)726 void ChMaterialShellReissner::ComputeStiffnessMatrix(ChMatrixRef K,
727 const ChVector<>& eps_u,
728 const ChVector<>& eps_v,
729 const ChVector<>& kur_u,
730 const ChVector<>& kur_v,
731 const double z_inf,
732 const double z_sup,
733 const double angle,
734 const ChShellReissnerInternalData* mdata) {
735 if (!plasticity || !mdata)
736 this->elasticity->ComputeStiffnessMatrix(K, eps_u, eps_v, kur_u, kur_v, z_inf, z_sup, angle);
737 else {
738 this->plasticity->ComputeStiffnessMatrixElastoplastic(K, eps_u, eps_v, kur_u, kur_v, *mdata, z_inf, z_sup, angle);
739 }
740 }
741
SetElasticity(std::shared_ptr<ChElasticityReissner> melasticity)742 void ChMaterialShellReissner::SetElasticity(std::shared_ptr<ChElasticityReissner> melasticity) {
743 elasticity = melasticity;
744 elasticity->section = this;
745 }
746
SetPlasticity(std::shared_ptr<ChPlasticityReissner> mplasticity)747 void ChMaterialShellReissner::SetPlasticity(std::shared_ptr<ChPlasticityReissner> mplasticity) {
748 plasticity = mplasticity;
749 mplasticity->section = this;
750 }
751
SetDamping(std::shared_ptr<ChDampingReissner> mdamping)752 void ChMaterialShellReissner::SetDamping(std::shared_ptr<ChDampingReissner> mdamping) {
753 damping = mdamping;
754 damping->section = this;
755 }
756
757
758
759
760
761
762
763
764 } // end of namespace fea
765 } // end of namespace chrono
766