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: Bryan Peterson, Milad Rakhsha, Antonio Recuero, Radu Serban 13 // ============================================================================= 14 // ANCF laminated shell element with four nodes. 15 // ============================================================================= 16 17 #ifndef CHELEMENTSHELLANCF3423_H 18 #define CHELEMENTSHELLANCF3423_H 19 20 #include <vector> 21 22 #include "chrono/fea/ChElementShell.h" 23 #include "chrono/fea/ChMaterialShellANCF.h" 24 #include "chrono/fea/ChNodeFEAxyzD.h" 25 26 namespace chrono { 27 namespace fea { 28 29 /// @addtogroup fea_elements 30 /// @{ 31 32 /// ANCF laminated shell element with four nodes. 33 /// This class implements composite material elastic force formulations. 34 /// 35 /// The node numbering is in ccw fashion as in the following scheme: 36 /// <pre> 37 /// v 38 /// ^ 39 /// | 40 /// D o-----+-----o C 41 /// | | | 42 /// --+-----+-----+----> u 43 /// | | | 44 /// A o-----+-----o B 45 /// </pre> 46 class ChApi ChElementShellANCF_3423 : public ChElementShell, public ChLoadableUV, public ChLoadableUVW { 47 public: 48 static const int NSF = 8; ///< number of shape functions 49 50 using ShapeVector = ChMatrixNM<double, 1, NSF>; 51 using VectorN = ChVectorN<double, NSF>; 52 using MatrixNx3 = ChMatrixNM<double, NSF, 3>; 53 54 ChElementShellANCF_3423(); ~ChElementShellANCF_3423()55 ~ChElementShellANCF_3423() {} 56 57 /// Definition of a layer 58 class Layer { 59 public: 60 /// Return the layer thickness. Get_thickness()61 double Get_thickness() const { return m_thickness; } 62 63 /// Return the fiber angle. Get_theta()64 double Get_theta() const { return m_theta; } 65 66 /// Return the layer material. GetMaterial()67 std::shared_ptr<ChMaterialShellANCF> GetMaterial() const { return m_material; } 68 69 private: 70 /// Private constructor (a layer can be created only by adding it to an element) 71 Layer(ChElementShellANCF_3423* element, ///< containing element 72 double thickness, ///< layer thickness 73 double theta, ///< fiber angle 74 std::shared_ptr<ChMaterialShellANCF> material ///< layer material 75 ); 76 Get_detJ0C()77 double Get_detJ0C() const { return m_detJ0C; } Get_T0()78 const ChMatrixNM<double, 6, 6>& Get_T0() const { return m_T0; } 79 80 /// Initial setup for this layer: calculate T0 and detJ0 at the element center. 81 void SetupInitial(); 82 83 ChElementShellANCF_3423* m_element; ///< containing ANCF shell element 84 std::shared_ptr<ChMaterialShellANCF> m_material; ///< layer material 85 double m_thickness; ///< layer thickness 86 double m_theta; ///< fiber angle 87 88 double m_detJ0C; 89 ChMatrixNM<double, 6, 6> m_T0; 90 91 public: 92 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 93 94 friend class ChElementShellANCF_3423; 95 friend class ShellANCF_Force; 96 friend class ShellANCF_Jacobian; 97 }; 98 99 /// Get the number of nodes used by this element. GetNnodes()100 virtual int GetNnodes() override { return 4; } 101 102 /// Get the number of coordinates in the field used by the referenced nodes. GetNdofs()103 virtual int GetNdofs() override { return 4 * 6; } 104 105 /// Get the number of coordinates from the n-th node used by this element. GetNodeNdofs(int n)106 virtual int GetNodeNdofs(int n) override { return 6; } 107 108 /// Specify the nodes of this element. 109 void SetNodes(std::shared_ptr<ChNodeFEAxyzD> nodeA, 110 std::shared_ptr<ChNodeFEAxyzD> nodeB, 111 std::shared_ptr<ChNodeFEAxyzD> nodeC, 112 std::shared_ptr<ChNodeFEAxyzD> nodeD); 113 114 /// Specify the element dimensions. SetDimensions(double lenX,double lenY)115 void SetDimensions(double lenX, double lenY) { 116 m_lenX = lenX; 117 m_lenY = lenY; 118 } 119 120 /// Access the n-th node of this element. GetNodeN(int n)121 virtual std::shared_ptr<ChNodeFEAbase> GetNodeN(int n) override { return m_nodes[n]; } 122 123 /// Get a handle to the first node of this element. GetNodeA()124 std::shared_ptr<ChNodeFEAxyzD> GetNodeA() const { return m_nodes[0]; } 125 126 /// Get a handle to the second node of this element. GetNodeB()127 std::shared_ptr<ChNodeFEAxyzD> GetNodeB() const { return m_nodes[1]; } 128 129 /// Get a handle to the third node of this element. GetNodeC()130 std::shared_ptr<ChNodeFEAxyzD> GetNodeC() const { return m_nodes[2]; } 131 132 /// Get a handle to the fourth node of this element. GetNodeD()133 std::shared_ptr<ChNodeFEAxyzD> GetNodeD() const { return m_nodes[3]; } 134 135 /// Add a layer. 136 void AddLayer(double thickness, ///< layer thickness 137 double theta, ///< fiber angle (radians) 138 std::shared_ptr<ChMaterialShellANCF> material ///< layer material 139 ); 140 141 /// Get the number of layers. GetNumLayers()142 size_t GetNumLayers() const { return m_numLayers; } 143 144 /// Get a handle to the specified layer. GetLayer(size_t i)145 const Layer& GetLayer(size_t i) const { return m_layers[i]; } 146 147 /// Set the structural damping. SetAlphaDamp(double a)148 void SetAlphaDamp(double a) { m_Alpha = a; } 149 150 /// Get the element length in the X direction. GetLengthX()151 double GetLengthX() const { return m_lenX; } 152 /// Get the element length in the Y direction. GetLengthY()153 double GetLengthY() const { return m_lenY; } 154 /// Get the total thickness of the shell element. GetThickness()155 double GetThickness() { return m_thickness; } 156 157 // Shape functions 158 // --------------- 159 160 /// Fills the N shape function matrix. 161 /// NOTE! actually N should be a 3row, 24 column sparse matrix, 162 /// as N = [s1*eye(3) s2*eye(3) s3*eye(3) s4*eye(3)...]; , 163 /// but to avoid wasting zero and repeated elements, here 164 /// it stores only the s1 through s8 values in a 1 row, 8 columns matrix! 165 void ShapeFunctions(ShapeVector& N, double x, double y, double z); 166 167 /// Fills the Nx shape function derivative matrix with respect to X. 168 /// NOTE! to avoid wasting zero and repeated elements, here 169 /// it stores only the four values in a 1 row, 8 columns matrix! 170 void ShapeFunctionsDerivativeX(ShapeVector& Nx, double x, double y, double z); 171 172 /// Fills the Ny shape function derivative matrix with respect to Y. 173 /// NOTE! to avoid wasting zero and repeated elements, here 174 /// it stores only the four values in a 1 row, 8 columns matrix! 175 void ShapeFunctionsDerivativeY(ShapeVector& Ny, double x, double y, double z); 176 177 /// Fills the Nz shape function derivative matrix with respect to Z. 178 /// NOTE! to avoid wasting zero and repeated elements, here 179 /// it stores only the four values in a 1 row, 8 columns matrix! 180 void ShapeFunctionsDerivativeZ(ShapeVector& Nz, double x, double y, double z); 181 182 /// Return a struct with 6-component strain and stress vectors evaluated at a 183 /// given quadrature point and layer number. 184 ChStrainStress3D EvaluateSectionStrainStress(const ChVector<>& loc, int layer_id); 185 void EvaluateDeflection(double& defVec); 186 187 public: 188 // Interface to ChElementBase base class 189 // ------------------------------------- 190 191 /// Fill the D vector with the current field values at the nodes of the element, with proper ordering. 192 /// If the D vector has not the size of this->GetNdofs(), it will be resized. 193 /// {x_a y_a z_a Dx_a Dx_a Dx_a x_b y_b z_b Dx_b Dy_b Dz_b} 194 virtual void GetStateBlock(ChVectorDynamic<>& mD) override; 195 196 // Set H as a linear combination of M, K, and R. 197 // H = Mfactor * [M] + Kfactor * [K] + Rfactor * [R], 198 // where [M] is the mass matrix, [K] is the stiffness matrix, and [R] is the damping matrix. 199 virtual void ComputeKRMmatricesGlobal(ChMatrixRef H, 200 double Kfactor, 201 double Rfactor = 0, 202 double Mfactor = 0) override; 203 204 /// Set M as the global mass matrix. 205 virtual void ComputeMmatrixGlobal(ChMatrixRef M) override; 206 207 /// Add contribution of element inertia to total nodal masses 208 virtual void ComputeNodalMass() override; 209 210 /// Compute the generalized force vector due to gravity using the efficient ANCF specific method 211 virtual void ComputeGravityForces(ChVectorDynamic<>& Fg, const ChVector<>& G_acc) override; 212 213 /// Computes the internal forces. 214 /// (E.g. the actual position of nodes is not in relaxed reference position) and set values in the Fi vector. 215 virtual void ComputeInternalForces(ChVectorDynamic<>& Fi) override; 216 217 /// Update the state of this element. 218 virtual void Update() override; 219 220 // Interface to ChElementShell base class 221 // -------------------------------------- 222 223 virtual void EvaluateSectionDisplacement(const double u, 224 const double v, 225 ChVector<>& u_displ, 226 ChVector<>& u_rotaz) override; 227 228 virtual void EvaluateSectionFrame(const double u, const double v, ChVector<>& point, ChQuaternion<>& rot) override; 229 230 virtual void EvaluateSectionPoint(const double u, const double v, ChVector<>& point) override; 231 232 // Internal computations 233 // --------------------- 234 235 /// Compute Jacobians of the internal forces. 236 /// This function calculates a linear combination of the stiffness (K) and damping (R) matrices, 237 /// J = Kfactor * K + Rfactor * R 238 /// for given coefficients Kfactor and Rfactor. 239 /// This Jacobian will be further combined with the global mass matrix M and included in the global 240 /// stiffness matrix H in the function ComputeKRMmatricesGlobal(). 241 void ComputeInternalJacobians(double Kfactor, double Rfactor); 242 243 /// Compute the mass matrix of the element. 244 /// Note: in this 'basic' implementation, constant section and 245 /// constant material are assumed 246 void ComputeMassMatrix(); 247 248 /// Compute the matrix to scale gravity by to get the generalized gravitational force. 249 void ComputeGravityForceScale(); 250 251 // [ANS] Shape function for Assumed Naturals Strain (Interpolation of strain and strainD in a thickness direction) 252 void ShapeFunctionANSbilinearShell(ChMatrixNM<double, 1, 4>& S_ANS, double x, double y); 253 254 // [ANS] Calculate the ANS strain and strain derivatives. 255 void CalcStrainANSbilinearShell(); 256 257 // [EAS] Basis function of M for Enhanced Assumed Strain. 258 void Basis_M(ChMatrixNM<double, 6, 5>& M, double x, double y, double z); 259 260 // Calculate the determinant of the initial configuration position vector gradient matrix 261 // at the specified point. 262 double Calc_detJ0(double x, double y, double z); 263 264 // Same as above, but also return the dense shape function vector derivatives. 265 double Calc_detJ0(double x, 266 double y, 267 double z, 268 ShapeVector& Nx, 269 ShapeVector& Ny, 270 ShapeVector& Nz, 271 ChMatrixNM<double, 1, 3>& Nx_d0, 272 ChMatrixNM<double, 1, 3>& Ny_d0, 273 ChMatrixNM<double, 1, 3>& Nz_d0); 274 275 // Calculate the current 8x3 matrix of nodal coordinates. 276 void CalcCoordMatrix(ChMatrixNM<double, 8, 3>& d); 277 278 // Calculate the current 24x1 matrix of nodal coordinate derivatives. 279 void CalcCoordDerivMatrix(ChVectorN<double, 24>& dt); 280 281 // Functions for ChLoadable interface 282 // ---------------------------------- 283 284 /// Gets the number of DOFs affected by this element (position part). LoadableGet_ndof_x()285 virtual int LoadableGet_ndof_x() override { return 4 * 6; } 286 287 /// Gets the number of DOFs affected by this element (velocity part). LoadableGet_ndof_w()288 virtual int LoadableGet_ndof_w() override { return 4 * 6; } 289 290 /// Gets all the DOFs packed in a single vector (position part). 291 virtual void LoadableGetStateBlock_x(int block_offset, ChState& mD) override; 292 293 /// Gets all the DOFs packed in a single vector (velocity part). 294 virtual void LoadableGetStateBlock_w(int block_offset, ChStateDelta& mD) override; 295 296 /// Increment all DOFs using a delta. 297 virtual void LoadableStateIncrement(const unsigned int off_x, 298 ChState& x_new, 299 const ChState& x, 300 const unsigned int off_v, 301 const ChStateDelta& Dv) override; 302 303 /// Number of coordinates in the interpolated field, ex=3 for a 304 /// tetrahedron finite element or a cable, = 1 for a thermal problem, etc. Get_field_ncoords()305 virtual int Get_field_ncoords() override { return 6; } 306 307 /// Get the number of DOFs sub-blocks. GetSubBlocks()308 virtual int GetSubBlocks() override { return 4; } 309 310 /// Get the offset of the specified sub-block of DOFs in global vector. GetSubBlockOffset(int nblock)311 virtual unsigned int GetSubBlockOffset(int nblock) override { return m_nodes[nblock]->NodeGetOffset_w(); } 312 313 /// Get the size of the specified sub-block of DOFs in global vector. GetSubBlockSize(int nblock)314 virtual unsigned int GetSubBlockSize(int nblock) override { return 6; } 315 316 /// Check if the specified sub-block of DOFs is active. IsSubBlockActive(int nblock)317 virtual bool IsSubBlockActive(int nblock) const override { return !m_nodes[nblock]->GetFixed(); } 318 319 virtual void EvaluateSectionVelNorm(double U, double V, ChVector<>& Result) override; 320 321 /// Get the pointers to the contained ChVariables, appending to the mvars vector. 322 virtual void LoadableGetVariables(std::vector<ChVariables*>& mvars) override; 323 324 /// Evaluate N'*F , where N is some type of shape function 325 /// evaluated at U,V coordinates of the surface, each ranging in -1..+1 326 /// F is a load, N'*F is the resulting generalized load 327 /// Returns also det[J] with J=[dx/du,..], that might be useful in gauss quadrature. 328 /// For this ANCF element, only the first 6 entries in F are used in the calculation. The first three entries is 329 /// the applied force in global coordinates and the second 3 entries is the applied moment in global space. 330 virtual void ComputeNF(const double U, ///< parametric coordinate in surface 331 const double V, ///< parametric coordinate in surface 332 ChVectorDynamic<>& Qi, ///< Return result of Q = N'*F here 333 double& detJ, ///< Return det[J] here 334 const ChVectorDynamic<>& F, ///< Input F vector, size is =n. field coords. 335 ChVectorDynamic<>* state_x, ///< if != 0, update state (pos. part) to this, then evaluate Q 336 ChVectorDynamic<>* state_w ///< if != 0, update state (speed part) to this, then evaluate Q 337 ) override; 338 339 /// Evaluate N'*F , where N is some type of shape function 340 /// evaluated at U,V,W coordinates of the volume, each ranging in -1..+1 341 /// F is a load, N'*F is the resulting generalized load 342 /// Returns also det[J] with J=[dx/du,..], that might be useful in gauss quadrature. 343 /// For this ANCF element, only the first 6 entries in F are used in the calculation. The first three entries is 344 /// the applied force in global coordinates and the second 3 entries is the applied moment in global space. 345 virtual void ComputeNF(const double U, ///< parametric coordinate in volume 346 const double V, ///< parametric coordinate in volume 347 const double W, ///< parametric coordinate in volume 348 ChVectorDynamic<>& Qi, ///< Return result of N'*F here, maybe with offset block_offset 349 double& detJ, ///< Return det[J] here 350 const ChVectorDynamic<>& F, ///< Input F vector, size is = n.field coords. 351 ChVectorDynamic<>* state_x, ///< if != 0, update state (pos. part) to this, then evaluate Q 352 ChVectorDynamic<>* state_w ///< if != 0, update state (speed part) to this, then evaluate Q 353 ) override; 354 355 /// This is needed so that it can be accessed by ChLoaderVolumeGravity. 356 /// Density is mass per unit surface. 357 virtual double GetDensity() override; 358 359 /// Gets the normal to the surface at the parametric coordinate U,V. 360 /// Each coordinate ranging in -1..+1. 361 virtual ChVector<> ComputeNormal(const double U, const double V) override; 362 363 private: 364 /// Initial setup. This is used to precompute matrices that do not change during the simulation, such as the local 365 /// stiffness of each element (if any), the mass, etc. 366 virtual void SetupInitial(ChSystem* system) override; 367 368 //// RADU 369 //// Why is m_d_dt inconsistent with m_d? Why not keep it as an 8x3 matrix? 370 371 std::vector<std::shared_ptr<ChNodeFEAxyzD>> m_nodes; ///< element nodes 372 std::vector<Layer, Eigen::aligned_allocator<Layer>> m_layers; ///< element layers 373 size_t m_numLayers; ///< number of layers for this element 374 double m_lenX; ///< element length in X direction 375 double m_lenY; ///< element length in Y direction 376 double m_thickness; ///< total element thickness 377 std::vector<double> m_GaussZ; ///< layer separation z values (scaled to [-1,1]) 378 double m_GaussScaling; ///< scaling factor due to change of integration intervals 379 double m_Alpha; ///< structural damping 380 VectorN m_GravForceScale; ///< Gravity scaling matrix used to get the generalized force due to gravity 381 ChMatrixNM<double, 24, 24> m_MassMatrix; ///< mass matrix 382 ChMatrixNM<double, 24, 24> m_JacobianMatrix; ///< Jacobian matrix (Kfactor*[K] + Rfactor*[R]) 383 ChMatrixNM<double, 8, 3> m_d0; ///< initial nodal coordinates 384 ChMatrixNM<double, 8, 8> m_d0d0T; ///< matrix m_d0 * m_d0^T 385 ChMatrixNM<double, 8, 3> m_d; ///< current nodal coordinates 386 ChMatrixNM<double, 8, 8> m_ddT; ///< matrix m_d * m_d^T 387 ChVectorN<double, 24> m_d_dt; ///< current nodal velocities 388 ChVectorN<double, 8> m_strainANS; ///< ANS strain 389 ChMatrixNM<double, 8, 24> m_strainANS_D; ///< ANS strain derivatives 390 std::vector<ChVectorN<double, 5>> m_alphaEAS; ///< EAS parameters (5 per layer) 391 std::vector<ChMatrixNM<double, 5, 5>> m_KalphaEAS; ///< EAS Jacobians (a 5x5 matrix per layer) 392 static const double m_toleranceEAS; ///< tolerance for nonlinear EAS solver (on residual) 393 static const int m_maxIterationsEAS; ///< maximum number of nonlinear EAS iterations 394 395 public: 396 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 397 398 friend class ShellANCF_Mass; 399 friend class ShellANCF_Gravity; 400 friend class ShellANCF_Force; 401 friend class ShellANCF_Jacobian; 402 }; 403 404 /// @} fea_elements 405 406 } // end of namespace fea 407 } // end of namespace chrono 408 409 #endif 410