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 13 // ============================================================================= 14 // Material for 6-field Reissner shells 15 // ============================================================================= 16 17 #ifndef CHMATERIALSHELLREISSNER_H 18 #define CHMATERIALSHELLREISSNER_H 19 20 #include <array> 21 #include <vector> 22 23 #include "chrono/fea/ChElementShell.h" 24 25 namespace chrono { 26 namespace fea { 27 28 /// @addtogroup fea_elements 29 /// @{ 30 31 32 // forward 33 class ChMaterialShellReissner; 34 35 36 37 /// Base interface for elasticity of 6-field Reissner-Mindlin shells (kinematically-exact shell theory 38 /// as in Witkowski et al.) to be used in a ChMaterialShellReissner. 39 /// Children classes must implement the ComputeStress function to get 40 /// {n_u,n_v,m_u,m_v}=f({e_u,e_v,k_u,k_v}) 41 /// Inherited materials do not define any thickness, which should be 42 /// a property of the element or its layer(s) using this material. 43 44 class ChApi ChElasticityReissner { 45 public: ChElasticityReissner()46 ChElasticityReissner() : section(nullptr) {} 47 ~ChElasticityReissner()48 virtual ~ChElasticityReissner() {} 49 50 /// Compute the generalized force and torque, given actual deformation and curvature. 51 /// This MUST be implemented by subclasses. 52 virtual void ComputeStress(ChVector<>& n_u, ///< forces along \e u direction (per unit length) 53 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 54 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 55 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 56 const ChVector<>& eps_u, ///< strains along \e u direction 57 const ChVector<>& eps_v, ///< strains along \e v direction 58 const ChVector<>& kur_u, ///< curvature along \e u direction 59 const ChVector<>& kur_v, ///< curvature along \e v direction 60 const double z_inf, ///< layer lower z value (along thickness coord) 61 const double z_sup, ///< layer upper z value (along thickness coord) 62 const double angle ///< layer angle respect to x (if needed) 63 ) = 0; 64 65 /// Compute the 12x12 stiffness matrix [Km] , that is [ds/de], the tangent of the constitutive relation 66 /// stresses/strains. 67 /// By default, it is computed by backward differentiation from the ComputeStress() function, 68 /// but inherited classes should better provide an analytical form, if possible. 69 virtual void ComputeStiffnessMatrix(ChMatrixRef mC, ///< tangent matrix 70 const ChVector<>& eps_u, ///< strains along \e u direction 71 const ChVector<>& eps_v, ///< strains along \e v direction 72 const ChVector<>& kur_u, ///< curvature along \e u direction 73 const ChVector<>& kur_v, ///< curvature along \e v direction 74 const double z_inf, ///< layer lower z value (along thickness coord) 75 const double z_sup, ///< layer upper z value (along thickness coord) 76 const double angle ///< layer angle respect to x (if needed) 77 ); 78 79 ChMaterialShellReissner* section; 80 }; 81 82 83 84 /// Elasticity of 6-field Reissner-Mindlin shells (kinematically-exact shell theory 85 /// as in Witkowski et al.) to be used in a ChMaterialShellReissner. 86 /// This class implements material properties for a layer from the Reissner theory, 87 /// for the case of isotropic linear linear elastic material. 88 /// This is probably the material that you need most often when using 6-field shells. 89 /// Previously: ChMaterialShellReissnerIsothropic 90 91 class ChApi ChElasticityReissnerIsothropic : public ChElasticityReissner { 92 public: 93 /// Construct an isotropic material. 94 ChElasticityReissnerIsothropic (double E, ///< Young's modulus 95 double nu, ///< Poisson ratio 96 double alpha = 1.0, ///< shear factor 97 double beta = 0.1 ///< torque factor 98 ); 99 100 /// Return the elasticity moduli Get_E()101 double Get_E() const { return m_E; } 102 /// Return the Poisson ratio Get_nu()103 double Get_nu() const { return m_nu; } 104 /// Return the shear factor Get_alpha()105 double Get_alpha() const { return m_alpha; } 106 /// Return the torque factor Get_beta()107 double Get_beta() const { return m_beta; } 108 109 /// The FE code will evaluate this function to compute 110 /// per-unit-length forces/torques given the u,v strains/curvatures. 111 virtual void ComputeStress( 112 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 113 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 114 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 115 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 116 const ChVector<>& eps_u, ///< strains along \e u direction 117 const ChVector<>& eps_v, ///< strains along \e v direction 118 const ChVector<>& kur_u, ///< curvature along \e u direction 119 const ChVector<>& kur_v, ///< curvature along \e v direction 120 const double z_inf, ///< layer lower z value (along thickness coord) 121 const double z_sup, ///< layer upper z value (along thickness coord) 122 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 123 ); 124 125 /// Compute 12x12 stiffness matrix [Km] , that is [ds/de], the tangent of the constitutive relation 126 /// per-unit-length forces/torques vs generalized strains. 127 virtual void ComputeStiffnessMatrix( 128 ChMatrixRef mC, ///< tangent matrix 129 const ChVector<>& eps_u, ///< strains along \e u direction 130 const ChVector<>& eps_v, ///< strains along \e v direction 131 const ChVector<>& kur_u, ///< curvature along \e u direction 132 const ChVector<>& kur_v, ///< curvature along \e v direction 133 const double z_inf, ///< layer lower z value (along thickness coord) 134 const double z_sup, ///< layer upper z value (along thickness coord) 135 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 136 ); 137 138 private: 139 double m_E; ///< elasticity moduli 140 double m_nu; ///< Poisson ratio 141 double m_alpha; ///< shear factor 142 double m_beta; ///< torque factor 143 }; 144 145 146 /// Elasticity of 6-field Reissner-Mindlin shells (kinematically-exact shell theory 147 /// as in Witkowski et al.) to be used in a ChMaterialShellReissner. 148 /// This class implements material properties for a layer from the Reissner theory, 149 /// for the case of orthotropic linear elastic material. 150 /// This is useful for laminated shells. One direction can be made softer than the other. 151 /// Note that the angle and the thickness are defined when adding a material with this elasticity to 152 /// a shell finite element (ex. ChElementShellReissner4) as a layer. 153 /// Previously: ChMaterialShellReissnerOrthotropic 154 155 class ChApi ChElasticityReissnerOrthotropic : public ChElasticityReissner { 156 public: 157 /// Construct an orthotropic material 158 ChElasticityReissnerOrthotropic (double m_E_x, ///< Young's modulus on x 159 double m_E_y, ///< Young's modulus on y 160 double m_nu_xy, ///< Poisson ratio xy (for yx it holds: nu_yx*E_x = nu_xy*E_y) 161 double m_G_xy, ///< Shear modulus, in plane 162 double m_G_xz, ///< Shear modulus, transverse 163 double m_G_yz, ///< Shear modulus, transverse 164 double m_alpha = 1.0, ///< shear factor 165 double m_beta = 0.1 ///< torque factor 166 ); 167 /// Construct an orthotropic material as sub case isotropic 168 ChElasticityReissnerOrthotropic (double m_E, ///< Young's modulus on x 169 double m_nu, ///< Poisson ratio 170 double m_alpha = 1.0, ///< shear factor 171 double m_beta = 0.1 ///< torque factor 172 ); 173 174 /// Return the elasticity moduli, on x Get_E_x()175 double Get_E_x() const { return E_x; } 176 /// Return the elasticity moduli, on y Get_E_y()177 double Get_E_y() const { return E_y; } 178 /// Return the Poisson ratio, for xy Get_nu_xy()179 double Get_nu_xy() const { return nu_xy; } 180 /// Return the Poisson ratio, for yx (follows xy as it must be nu_yx*E_x = nu_xy*E_y) Get_nu_yx()181 double Get_nu_yx() const { return nu_xy * (E_y / E_x); } 182 /// Return the shear mod, in plane Get_G_xy()183 double Get_G_xy() const { return G_xy; } 184 /// Return the shear mod, transverse Get_G_xz()185 double Get_G_xz() const { return G_xz; } 186 /// Return the shear mod, transverse Get_G_yz()187 double Get_G_yz() const { return G_yz; } 188 /// Return the shear factor Get_alpha()189 double Get_alpha() const { return alpha; } 190 /// Return the torque factor Get_beta()191 double Get_beta() const { return beta; } 192 193 /// The FE code will evaluate this function to compute 194 /// per-unit-length forces/torques given the u,v strains/curvatures. 195 virtual void ComputeStress( 196 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 197 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 198 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 199 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 200 const ChVector<>& eps_u, ///< strains along \e u direction 201 const ChVector<>& eps_v, ///< strains along \e v direction 202 const ChVector<>& kur_u, ///< curvature along \e u direction 203 const ChVector<>& kur_v, ///< curvature along \e v direction 204 const double z_inf, ///< layer lower z value (along thickness coord) 205 const double z_sup, ///< layer upper z value (along thickness coord) 206 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 207 ); 208 209 /// /// Compute the 12x12 stiffness matrix [Km] , that is [ds/de], the tangent of the constitutive relation 210 /// stresses/strains. 211 virtual void ComputeStiffnessMatrix( 212 ChMatrixRef mC, ///< tangent matrix 213 const ChVector<>& eps_u, ///< strains along \e u direction 214 const ChVector<>& eps_v, ///< strains along \e v direction 215 const ChVector<>& kur_u, ///< curvature along \e u direction 216 const ChVector<>& kur_v, ///< curvature along \e v direction 217 const double z_inf, ///< layer lower z value (along thickness coord) 218 const double z_sup, ///< layer upper z value (along thickness coord) 219 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 220 ); 221 222 private: 223 double E_x; ///< elasticity moduli 224 double E_y; ///< elasticity moduli 225 double nu_xy; ///< Poisson ratio 226 double G_xy; ///< Shear factor, in plane 227 double G_xz; ///< Shear factor, out of plane 228 double G_yz; ///< Shear factor, out of plane 229 double alpha; ///< shear factor 230 double beta; ///< torque factor 231 }; 232 233 234 // ---------------------------------------------------------------------------- 235 236 /// Generic linear elasticity for 6-field Reissner-Mindlin shells (kinematically-exact shell theory 237 /// as in Witkowski et al.) to be used in a ChMaterialShellReissner. 238 /// This uses a 12x12 matrix [E] from user-input data. The [E] matrix can be 239 /// computed from a preprocessing stage using a FEA analysis over a detailed 3D model 240 /// of a slab of shell, hence recovering the 6x6 matrix in the linear mapping: 241 /// {n,m}=[E]{e,k}. 242 243 class ChApi ChElasticityReissnerGeneric : public ChElasticityReissner { 244 public: 245 ChElasticityReissnerGeneric(); 246 ~ChElasticityReissnerGeneric()247 virtual ~ChElasticityReissnerGeneric() {} 248 249 /// Access the 12x12 [E] matrix, for getting/setting its values. 250 /// This is the matrix that defines the linear elastic constitutive model 251 /// as it maps yxz displacements "e" and xyz rotations "k" 252 /// to the "n" force and "m" torque as in 253 /// {n_u,n_v,m_u,m_v}=[E]{e_u,e_v,k_u,k_v}. Ematrix()254 ChMatrixNM<double, 12, 12>& Ematrix() { return this->mE; } 255 256 /// The FE code will evaluate this function to compute 257 /// per-unit-length forces/torques given the u,v strains/curvatures. 258 virtual void ComputeStress( 259 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 260 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 261 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 262 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 263 const ChVector<>& eps_u, ///< strains along \e u direction 264 const ChVector<>& eps_v, ///< strains along \e v direction 265 const ChVector<>& kur_u, ///< curvature along \e u direction 266 const ChVector<>& kur_v, ///< curvature along \e v direction 267 const double z_inf, ///< layer lower z value (along thickness coord) 268 const double z_sup, ///< layer upper z value (along thickness coord) 269 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 270 ) override; 271 272 /// /// Compute the 12x12 stiffness matrix [Km] , that is [ds/de], the tangent of the constitutive relation 273 /// stresses/strains. 274 virtual void ComputeStiffnessMatrix( 275 ChMatrixRef mC, ///< tangent matrix 276 const ChVector<>& eps_u, ///< strains along \e u direction 277 const ChVector<>& eps_v, ///< strains along \e v direction 278 const ChVector<>& kur_u, ///< curvature along \e u direction 279 const ChVector<>& kur_v, ///< curvature along \e v direction 280 const double z_inf, ///< layer lower z value (along thickness coord) 281 const double z_sup, ///< layer upper z value (along thickness coord) 282 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 283 ) override; 284 285 private: 286 ChMatrixNM<double, 12, 12> mE; 287 288 public: 289 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 290 }; 291 292 293 294 295 296 297 // ---------------------------------------------------------------------------- 298 299 /// Base class for internal variables of Reissner shells materials. 300 /// Especially useful for plasticity, where internal variables are used 301 /// to carry information on plastic flow, accumulated flow, etc. 302 class ChApi ChShellReissnerInternalData { 303 public: ChShellReissnerInternalData()304 ChShellReissnerInternalData() : p_strain_acc(0) {} 305 ~ChShellReissnerInternalData()306 virtual ~ChShellReissnerInternalData(){}; 307 Copy(const ChShellReissnerInternalData & other)308 virtual void Copy(const ChShellReissnerInternalData& other) { p_strain_acc = other.p_strain_acc; } 309 310 double p_strain_acc; // accumulated flow, \overbar\eps^p in Neto-Owen book 311 }; 312 313 /// Base interface for plasticity of 6-field Reissner-Mindlin shells (kinematically-exact shell theory 314 /// as in Witkowski et al.) to be used in a ChMaterialShellReissner. 315 /// Children classes must implement the ComputeStressWithReturnMapping to compute 316 /// effective stress and strain given a tentative strain that might violate the yeld function. 317 /// Inherited materials do not define any thickness, which should be 318 /// a property of the element or its layer(s) using this material. 319 class ChApi ChPlasticityReissner { 320 public: 321 ChPlasticityReissner(); 322 ~ChPlasticityReissner()323 virtual ~ChPlasticityReissner() {} 324 325 /// Given a trial strain, it computes the effective stress and strain by 326 /// clamping against the yeld surface. An implicit return mapping integration 327 /// step is computed automatically per each call of this function. 328 /// Note: for the elastic part, it must use the elasticity model in this->section->elasticity. 329 /// If not beyond yeld, simply: 330 /// elastic strain = tot strain - plastic strain 331 /// If it is beyond yeld: 332 /// elastic strain is computed by fully implicit strain integration with return mapping, 333 /// and plastic strains in "data_new" are updated. 334 /// Returns true if it had to do return mapping, false if it was in elastic regime 335 /// This MUST be implemented by subclasses. 336 virtual bool ComputeStressWithReturnMapping( 337 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 338 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 339 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 340 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 341 ChShellReissnerInternalData& data_new, ///< updated material internal variables, at this point, including {p_strain_e, p_strain_k, p_strain_acc} 342 const ChVector<>& eps_u_trial, ///< trial strains along \e u direction 343 const ChVector<>& eps_v_trial, ///< trial strains along \e v direction 344 const ChVector<>& kur_u_trial, ///< trial curvature along \e u direction 345 const ChVector<>& kur_v_trial, ///< trial curvature along \e v direction 346 const ChShellReissnerInternalData& data, ///< trial material internal variables, at this point, including {p_strain_e, p_strain_k, p_strain_acc} 347 const double z_inf, ///< layer lower z value (along thickness coord) 348 const double z_sup, ///< layer upper z value (along thickness coord) 349 const double angle ///< layer angle respect to x (if needed) 350 ) = 0; 351 352 /// Compute the 12x12 tangent material stiffness matrix [Km] = dσ/dε, 353 /// given actual internal data and deformation and curvature (if needed). If in 354 /// plastic regime, uses elastoplastic matrix, otherwise uses elastic. 355 /// This must be overridden by subclasses if an analytical solution is 356 /// known (preferred for high performance), otherwise the base behaviour here is to compute 357 /// [Km] by numerical differentiation calling ComputeStressWithReturnMapping() multiple times. 358 virtual void ComputeStiffnessMatrixElastoplastic( 359 ChMatrixRef K, ///< 12x12 material elastoplastic stiffness matrix values here 360 const ChVector<>& eps_u, ///< strains along \e u direction 361 const ChVector<>& eps_v, ///< strains along \e v direction 362 const ChVector<>& kur_u, ///< curvature along \e u direction 363 const ChVector<>& kur_v, ///< curvature along \e v direction 364 const ChShellReissnerInternalData& data, ///< updated material internal variables, at this point including {p_strain_e, p_strain_k, p_strain_acc} 365 const double z_inf, ///< layer lower z value (along thickness coord) 366 const double z_sup, ///< layer upper z value (along thickness coord) 367 const double angle ///< layer angle respect to x (if needed) 368 ); 369 370 // Populate a vector with the appropriate ChBeamSectionPlasticity data structures. 371 // Children classes may override this. By default uses ChBeamMaterialInternalData for basic plasticity. 372 // Thanks to unique_ptr there is no need to call delete for the pointed objects. 373 virtual void CreatePlasticityData(int numpoints, 374 std::vector<std::unique_ptr<ChShellReissnerInternalData>>& plastic_data); 375 376 377 ChMaterialShellReissner* section; 378 379 double nr_yeld_tolerance; 380 int nr_yeld_maxiters; 381 }; 382 383 384 // ---------------------------------------------------------------------------- 385 386 387 /// Base interface for damping of 6-field Reissner-Mindlin shells (kinematically-exact shell theory 388 /// as in Witkowski et al.) to be used in a ChMaterialShellReissner. 389 /// Children classes should implement a ComputeStress function that returns generalized stresses 390 /// given time derivatives of strains as: 391 /// {n_u,n_v,m_u.m_v}=f({e_u',e_v',k_u',k_v'}) 392 393 class ChApi ChDampingReissner { 394 public: ChDampingReissner()395 ChDampingReissner() : section(nullptr) {} 396 ~ChDampingReissner()397 virtual ~ChDampingReissner() {} 398 399 /// Compute the generalized cut force and cut torque, caused by structural damping, 400 /// given actual deformation speed and curvature speed. 401 /// This MUST be implemented by subclasses. 402 virtual void ComputeStress( 403 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 404 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 405 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 406 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 407 const ChVector<>& deps_u, ///< time derivative of strains along \e u direction 408 const ChVector<>& deps_v, ///< time derivative of strains along \e v direction 409 const ChVector<>& dkur_u, ///< time derivative of curvature along \e u direction 410 const ChVector<>& dkur_v, ///< time derivative of curvature along \e v direction 411 const double z_inf, ///< layer lower z value (along thickness coord) 412 const double z_sup, ///< layer upper z value (along thickness coord) 413 const double angle ///< layer angle respect to x (if needed) 414 ) = 0; 415 416 /// Compute the 12x12 tangent material damping matrix, ie the jacobian [Rm]=dstress/dstrainspeed. 417 /// This must be overridden by subclasses if an analytical solution is 418 /// known (preferred for high performance), otherwise the base behaviour here is to compute 419 /// [Rm] by numerical differentiation calling ComputeStress() multiple times. 420 virtual void ComputeDampingMatrix( ChMatrixRef R, ///< 12x12 material damping matrix values here 421 const ChVector<>& deps_u, ///< time derivative of strains along \e u direction 422 const ChVector<>& deps_v, ///< time derivative of strains along \e v direction 423 const ChVector<>& dkur_u, ///< time derivative of curvature along \e u direction 424 const ChVector<>& dkur_v, ///< time derivative of curvature along \e v direction 425 const double z_inf, ///< layer lower z value (along thickness coord) 426 const double z_sup, ///< layer upper z value (along thickness coord) 427 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 428 ); 429 430 ChMaterialShellReissner* section; 431 }; 432 433 434 435 436 /// Simple Rayleight damping of a Reissner-mindlin shell, 437 /// where damping is proportional to stiffness via a beta coefficient. 438 /// In order to generalize it also in case of nonlinearity, the full 439 /// element tangent stiffness matrix cannot be used (it may contain negative eigenvalues) 440 /// and it can't be used to recover instant nodal caused by damping as F=beta*K*q_dt 441 /// so it is generalized to the following implementation at the material stress level 442 /// <pre> 443 /// {n,m}=beta*[E]*{e',k'} 444 /// </pre> 445 /// where 446 /// - beta is the 2nd Rayleigh damping parameter 447 /// - [E] is the 6x6 shell stiffness matrix at the undeformed unstressed case (hence assumed constant) 448 /// - {e',k'} is the speed of deformation/curvature 449 /// Note that the alpha mass-proportional parameter (the first of the alpha,beta parameters of the original 450 /// Rayleigh model) is not supported. 451 452 class ChApi ChDampingReissnerRayleigh : public ChDampingReissner { 453 public: 454 /// Construct the Rayleigh damping model from the stiffness model used by the shell layer. 455 /// This is important because the Rayleigh damping is proportional to the stiffness, 456 /// so the model must know which is the stiffness matrix of the material. 457 /// Note: melasticity must be alreay set with proper values: its [E] stiffness matrix will be 458 /// fetched just once for all. 459 ChDampingReissnerRayleigh(std::shared_ptr<ChElasticityReissner> melasticity, const double& mbeta = 0); 460 ~ChDampingReissnerRayleigh()461 virtual ~ChDampingReissnerRayleigh() {} 462 463 /// Compute the generalized cut force and cut torque, caused by structural damping, 464 /// given actual deformation speed and curvature speed. 465 virtual void ComputeStress( 466 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 467 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 468 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 469 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 470 const ChVector<>& deps_u, ///< time derivative of strains along \e u direction 471 const ChVector<>& deps_v, ///< time derivative of strains along \e v direction 472 const ChVector<>& dkur_u, ///< time derivative of curvature along \e u direction 473 const ChVector<>& dkur_v, ///< time derivative of curvature along \e v direction 474 const double z_inf, ///< layer lower z value (along thickness coord) 475 const double z_sup, ///< layer upper z value (along thickness coord) 476 const double angle ///< layer angle respect to x (if needed) 477 ); 478 479 /// Compute the 6x6 tangent material damping matrix, ie the jacobian [Rm]=dstress/dstrainspeed. 480 /// In this model, it is beta*[E] where [E] is the 12x12 stiffness matrix at material level, assumed constant 481 virtual void ComputeDampingMatrix( ChMatrixRef R, ///< 12x12 material damping matrix values here 482 const ChVector<>& deps_u, ///< time derivative of strains along \e u direction 483 const ChVector<>& deps_v, ///< time derivative of strains along \e v direction 484 const ChVector<>& dkur_u, ///< time derivative of curvature along \e u direction 485 const ChVector<>& dkur_v, ///< time derivative of curvature along \e v direction 486 const double z_inf, ///< layer lower z value (along thickness coord) 487 const double z_sup, ///< layer upper z value (along thickness coord) 488 const double angle ///< layer angle respect to x (if needed) -not used in this, isotropic 489 ); 490 491 /// Get the beta Rayleigh parameter (stiffness proportional damping) GetBeta()492 double GetBeta() { return beta; } 493 /// Set the beta Rayleigh parameter (stiffness proportional damping) SetBeta(const double mbeta)494 void SetBeta(const double mbeta) { beta = mbeta; } 495 496 497 private: 498 std::shared_ptr<ChElasticityReissner> section_elasticity; 499 ChMatrixNM<double, 12, 12> E_const; // to store the precomputed stiffness matrix at undeformed unstressed initial state 500 double beta; 501 bool updated; 502 503 public: 504 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 505 }; 506 507 508 509 510 // ---------------------------------------------------------------------------- 511 512 /// Material for a single layer of a 6-field Reissner-Mindlin shells 513 /// (kinematically-exact shell theory as in Witkowski et al). 514 /// This base implementation assumes that one creates a ChMaterialShellReissner 515 /// by providing three components: 516 /// 517 /// - an elasticity model (from ChElasticityReissner classes) 518 /// - a plasticity model (optional, from ChPlasticityReissner classes) 519 /// - a damping model (optional, from ChDampingReissner classes) 520 /// 521 /// Thickness is defined when adding a ChMaterialShellReissner material as a layer 522 /// in a shell finite element (ex. ChElementShellReissner4). 523 /// A material can be shared between multiple layers. 524 525 class ChApi ChMaterialShellReissner { 526 public: 527 ChMaterialShellReissner(std::shared_ptr<ChElasticityReissner> melasticity ///< elasticity model 528 ); 529 530 ChMaterialShellReissner( 531 std::shared_ptr<ChElasticityReissner> melasticity, ///< elasticity model 532 std::shared_ptr<ChPlasticityReissner> mplasticity ///< plasticity model, if any 533 ); 534 535 ChMaterialShellReissner( 536 std::shared_ptr<ChElasticityReissner> melasticity, ///< elasticity model 537 std::shared_ptr<ChPlasticityReissner> mplasticity, ///< plasticity model, if any 538 std::shared_ptr<ChDampingReissner> mdamping ///< damping model, if any 539 ); 540 ~ChMaterialShellReissner()541 virtual ~ChMaterialShellReissner() {} 542 543 /// Compute the generalized cut force and cut torque, given the actual generalized section strain 544 /// expressed as deformation vector e and curvature k, that is: {n_u,n_v,m_u,m_v}=f({e_u,e_v,k_u,k_v}), and 545 /// given the actual material state required for plasticity if any (but if mdata=nullptr, 546 /// computes only the elastic force). 547 /// If there is plasticity, the stress is clamped by automatically performing an implicit return mapping. 548 /// In sake of generality, if possible this is the function that should be used by beam finite elements 549 /// to compute internal forces, ex.by some Gauss quadrature. 550 virtual void ComputeStress( 551 ChVector<>& n_u, ///< forces along \e u direction (per unit length) 552 ChVector<>& n_v, ///< forces along \e v direction (per unit length) 553 ChVector<>& m_u, ///< torques along \e u direction (per unit length) 554 ChVector<>& m_v, ///< torques along \e v direction (per unit length) 555 const ChVector<>& eps_u, ///< strains along \e u direction 556 const ChVector<>& eps_v, ///< strains along \e v direction 557 const ChVector<>& kur_u, ///< curvature along \e u direction 558 const ChVector<>& kur_v, ///< curvature along \e v direction 559 const double z_inf, ///< layer lower z value (along thickness coord) 560 const double z_sup, ///< layer upper z value (along thickness coord) 561 const double angle, ///< layer angle respect to x (if needed) 562 ChShellReissnerInternalData* mdata_new = nullptr, ///< updated material internal variables, at this 563 ///< point, including {p_strain_e, p_strain_k, p_strain_acc} 564 const ChShellReissnerInternalData* mdata = nullptr ///< current material internal variables, at this point, 565 ///< including {p_strain_e, p_strain_k, p_strain_acc} 566 ); 567 568 /// Compute the 6x6 tangent material stiffness matrix [Km] = dσ/dε 569 /// at a given strain state, and at given internal data state (if mdata=nullptr, 570 /// computes only the elastic tangent stiffenss, regardless of plasticity). 571 virtual void ComputeStiffnessMatrix( 572 ChMatrixRef K, ///< 12x12 stiffness matrix 573 const ChVector<>& eps_u, ///< strains along \e u direction 574 const ChVector<>& eps_v, ///< strains along \e v direction 575 const ChVector<>& kur_u, ///< curvature along \e u direction 576 const ChVector<>& kur_v, ///< curvature along \e v direction 577 const double z_inf, ///< layer lower z value (along thickness coord) 578 const double z_sup, ///< layer upper z value (along thickness coord) 579 const double angle, ///< layer angle respect to x (if needed) 580 const ChShellReissnerInternalData* mdata = nullptr ///< material internal variables, at this point, if any, 581 ///< including {p_strain_e, p_strain_k, p_strain_acc} 582 ); 583 584 /// Set the elasticity model for this section. 585 /// By default it uses a simple centered linear elastic model, but you can set more complex models. 586 void SetElasticity(std::shared_ptr<ChElasticityReissner> melasticity); 587 588 /// Get the elasticity model for this section. 589 /// Use this function to access parameters such as stiffness, Young modulus, etc. 590 /// By default it uses a simple centered linear elastic model. GetElasticity()591 std::shared_ptr<ChElasticityReissner> GetElasticity() { return this->elasticity; } 592 593 /// Set the plasticity model for this section. 594 /// This is independent from the elasticity model. 595 /// Note that by default there is no plasticity model, 596 /// so by default plasticity never happens. 597 void SetPlasticity(std::shared_ptr<ChPlasticityReissner> mplasticity); 598 599 /// Get the elasticity model for this section, if any. 600 /// Use this function to access parameters such as yeld limit, etc. GetPlasticity()601 std::shared_ptr<ChPlasticityReissner> GetPlasticity() { return this->plasticity; } 602 603 /// Set the damping model for this section. 604 /// By default no damping. 605 void SetDamping(std::shared_ptr<ChDampingReissner> mdamping); 606 607 /// Get the damping model for this section. 608 /// By default no damping. GetDamping()609 std::shared_ptr<ChDampingReissner> GetDamping() { return this->damping; } 610 611 /// Set the density of the shell (kg/m^3) SetDensity(double md)612 void SetDensity(double md) { this->density = md; } GetDensity()613 double GetDensity() const { return this->density; } 614 615 ///for backward compatibility - use GetDensity instead Get_rho()616 double Get_rho() const { return this->density; } 617 618 private: 619 620 std::shared_ptr<ChElasticityReissner> elasticity; 621 std::shared_ptr<ChPlasticityReissner> plasticity; 622 std::shared_ptr<ChDampingReissner> damping; 623 624 double density; 625 }; 626 627 628 //------------------------------------------------------------------ 629 630 // 631 // ONLY FOR BACKWARD COMPATIBILITY 632 // 633 634 /// For backward compatibility only! 635 /// New approach: create a ChElasticityReissnerOrthotropic and create a ChMaterialShellReissner by passing the elasticity as a parameter. 636 637 class ChApi ChMaterialShellReissnerIsothropic : public ChMaterialShellReissner { 638 public: 639 /// Construct an isotropic material. 640 ChMaterialShellReissnerIsothropic(double mdensity, ///< material density 641 double E, ///< Young's modulus 642 double nu, ///< Poisson ratio 643 double alpha = 1.0, ///< shear factor 644 double beta = 0.1 ///< torque factor 645 ) ChMaterialShellReissner(chrono_types::make_shared<ChElasticityReissnerIsothropic> (E,nu,alpha,beta))646 : ChMaterialShellReissner(chrono_types::make_shared<ChElasticityReissnerIsothropic>(E, nu, alpha, beta)) { 647 this->SetDensity(mdensity); 648 } 649 650 }; 651 652 /// For backward compatibility only! 653 /// New approach: create a ChElasticityReissnerOrthotropic and create a ChMaterialShellReissner by passing the elasticity as a parameter. 654 655 class ChApi ChMaterialShellReissnerOrthotropic : public ChMaterialShellReissner { 656 public: 657 /// Construct an orthotropic material 658 ChMaterialShellReissnerOrthotropic(double mdensity, ///< material density 659 double m_E_x, ///< Young's modulus on x 660 double m_E_y, ///< Young's modulus on y 661 double m_nu_xy, ///< Poisson ratio xy (for yx it holds: nu_yx*E_x = nu_xy*E_y) 662 double m_G_xy, ///< Shear modulus, in plane 663 double m_G_xz, ///< Shear modulus, transverse 664 double m_G_yz, ///< Shear modulus, transverse 665 double m_alpha = 1.0, ///< shear factor 666 double m_beta = 0.1 ///< torque factor 667 ) ChMaterialShellReissner(chrono_types::make_shared<ChElasticityReissnerOrthotropic> (m_E_x,m_E_y,m_nu_xy,m_G_xy,m_G_xz,m_G_yz,m_alpha,m_beta))668 : ChMaterialShellReissner(chrono_types::make_shared<ChElasticityReissnerOrthotropic>(m_E_x, m_E_y, m_nu_xy, m_G_xy, m_G_xz, m_G_yz, m_alpha, m_beta)) 669 { 670 this->SetDensity(mdensity); 671 } 672 673 /// Construct an orthotropic material as sub case isotropic 674 ChMaterialShellReissnerOrthotropic(double mdensity, ///< material density 675 double m_E, ///< Young's modulus on x 676 double m_nu, ///< Poisson ratio 677 double m_alpha = 1.0, ///< shear factor 678 double m_beta = 0.1 ///< torque factor 679 ) ChMaterialShellReissner(chrono_types::make_shared<ChElasticityReissnerOrthotropic> (m_E,m_nu,m_alpha,m_beta))680 : ChMaterialShellReissner( 681 chrono_types::make_shared<ChElasticityReissnerOrthotropic>(m_E, m_nu, m_alpha, m_beta)) 682 { 683 this->SetDensity(mdensity); 684 } 685 686 }; 687 688 689 690 /// @} fea_elements 691 692 } // end of namespace fea 693 } // end of namespace chrono 694 695 #endif 696