ME751 Advanced Computational Multibody Dynamics November 23, 2016 - - PowerPoint PPT Presentation

▶
me751 advanced computational multibody dynamics
SMART_READER_LITE
LIVE PREVIEW

ME751 Advanced Computational Multibody Dynamics November 23, 2016 - - PowerPoint PPT Presentation

Mar 15, 2023 384 likes •725 views

ME751 Advanced Computational Multibody Dynamics November 23, 2016 Antonio Recuero University of Wisconsin-Madison Before we get started Last time, we learned: Plates and shells Introduce Chronos bilinear ANCF shell element

slide-1
SLIDE 1

ME751 Advanced Computational Multibody Dynamics

November 23, 2016

Antonio Recuero University of Wisconsin-Madison

slide-2
SLIDE 2

Before we get started…

Last time, we learned:

Plates and shells

Introduce Chrono’s bilinear ANCF shell element

Kinematics of ANCF shell element

Curvilinear initial configuration

Inertia forces

Internal forces

This lecture…

Homework

Convergence

Locking Issues

Assumed Natural Strain

Enhanced Assumed Strain

Implementation in Chrono

Examples

2

slide-3
SLIDE 3
  • 4. Shell strains: Orthotropic and

curvilinear reference

3

slide-4
SLIDE 4
  • 4. Shell strains: Orthotropic and

curvilinear reference

4

slide-5
SLIDE 5
  • 4. Shell strains: Orthotropic and

curvilinear reference

5

Current deformed reference

slide-6
SLIDE 6
  • 4. Shell strains: Orthotropic and

curvilinear reference

6

slide-7
SLIDE 7
  • 4. Shell strains: Orthotropic and

curvilinear reference

7

slide-8
SLIDE 8
  • 4. Shell strains: Orthotropic

and curvilinear reference

8

Final expression! These strains are related to constitutive equations

slide-9
SLIDE 9
  • 4. Shell strains: Orthotropic

and curvilinear reference

9

slide-10
SLIDE 10
  • 1. Locking

10

slide-11
SLIDE 11
  • 1. Locking

11

slide-12
SLIDE 12
  • 1. Locking issues: Convergence

12

slide-13
SLIDE 13

1.a Volumetric Locking

13

slide-14
SLIDE 14

1.a Volumetric Locking

14

Rigid body Strain modes Hourglass modes

A uniform strain hexahedron and quadrilateral with hourglass control: Flanagan and Belytschko, 1981

slide-15
SLIDE 15

1.b Shear Locking

15

slide-16
SLIDE 16
  • 2. ANCF bilinear shell element

16

slide-17
SLIDE 17
  • 2. ANCF bilinear shell. ANS

17

slide-18
SLIDE 18
  • 2. ANCF bilinear shell. ANS

18

Substitution

  • f these

compatible strains in equations

slide-19
SLIDE 19
  • 2. ANCF bilinear shell. ANS

19

slide-20
SLIDE 20
  • 2. ANCF bilinear shell. EAS

20

slide-21
SLIDE 21
  • 2. ANCF bilinear shell. EAS

21

slide-22
SLIDE 22
  • 2. ANCF bilinear shell. EAS

22

slide-23
SLIDE 23
  • 2. ANCF bilinear shell. EAS

23

Substitution of thickness strain

slide-24
SLIDE 24
  • 3. Implementation of bilinear

ANCF shell element in Chrono

24 void MyForce::Evaluate(ChMatrixNM<double, 54, 1>& result, const double x, const double y, const double z) { // Element shape function ChMatrixNM<double, 1, 8> N; m_element‐>ShapeFunctions(N, x, y, z); // Determinant of position vector gradient matrix: Initial configuration ChMatrixNM<double, 1, 8> Nx; ChMatrixNM<double, 1, 8> Ny; ChMatrixNM<double, 1, 8> Nz; ChMatrixNM<double, 1, 3> Nx_d0; ChMatrixNM<double, 1, 3> Ny_d0; ChMatrixNM<double, 1, 3> Nz_d0; double detJ0 = m_element‐>Calc_detJ0(x, y, z, Nx, Ny, Nz, Nx_d0, Ny_d0, Nz_d0); // ANS shape function ChMatrixNM<double, 1, 4> S_ANS; // Shape function vector for Assumed Natural Strain ChMatrixNM<double, 6, 5> M; // Shape function vector for Enhanced Assumed Strain m_element‐>ShapeFunctionANSbilinearShell(S_ANS, x, y); m_element‐>Basis_M(M, x, y, z); // Transformation : Orthogonal transformation (A and J) ChVector<double> G1xG2; // Cross product of first and second column of double G1dotG1; // Dot product of first column of position vector gradient G1xG2.x = Nx_d0[0][1] * Ny_d0[0][2] ‐ Nx_d0[0][2] * Ny_d0[0][1]; G1xG2.y = Nx_d0[0][2] * Ny_d0[0][0] ‐ Nx_d0[0][0] * Ny_d0[0][2]; G1xG2.z = Nx_d0[0][0] * Ny_d0[0][1] ‐ Nx_d0[0][1] * Ny_d0[0][0]; G1dotG1 = Nx_d0[0][0] * Nx_d0[0][0] + Nx_d0[0][1] * Nx_d0[0][1] + Nx_d0[0][2] * Nx_d0[0][2]; // Tangent Frame ChVector<double> A1; ChVector<double> A2; ChVector<double> A3; A1.x = Nx_d0[0][0]; A1.y = Nx_d0[0][1]; A1.z = Nx_d0[0][2]; A1 = A1 / sqrt(G1dotG1); A3 = G1xG2.GetNormalized(); A2.Cross(A3, A1); Evaluation of 3D shape functions Shape function derivatives Vectors of initial deformation gradient New interpolation for ANS Interpolation for 5 internal parameters Covariant vector base, initial config. Normalization

  • f such a frame
slide-25
SLIDE 25
  • 3. Internal forces

25 // Direction for orthotropic material double theta = m_element‐>GetLayer(m_kl).Get_theta(); // Fiber angle ChVector<double> AA1; ChVector<double> AA2; ChVector<double> AA3; AA1 = A1 * cos(theta) + A2 * sin(theta); AA2 = ‐A1 * sin(theta) + A2 * cos(theta); AA3 = A3; /// Beta ChMatrixNM<double, 3, 3> j0; ChVector<double> j01; ChVector<double> j02; ChVector<double> j03; ChMatrixNM<double, 9, 1> beta; // Calculates inverse of rd0 (j0) (position vector gradient: Initial Configuration) j0(0, 0) = Ny_d0[0][1] * Nz_d0[0][2] ‐ Nz_d0[0][1] * Ny_d0[0][2]; j0(0, 1) = Ny_d0[0][2] * Nz_d0[0][0] ‐ Ny_d0[0][0] * Nz_d0[0][2]; j0(0, 2) = Ny_d0[0][0] * Nz_d0[0][1] ‐ Nz_d0[0][0] * Ny_d0[0][1]; j0(1, 0) = Nz_d0[0][1] * Nx_d0[0][2] ‐ Nx_d0[0][1] * . . . j0.MatrDivScale(detJ0); j01(0) = j0(0, 0); j02(0) = j0(1, 0); j03(0) = j0(2, 0); . . . j02(2) = j0(1, 2); j03(2) = j0(2, 2); // Coefficients of contravariant transformation beta(0, 0) = Vdot(AA1, j01); beta(1, 0) = Vdot(AA2, j01); beta(2, 0) = Vdot(AA3, j01); . . . beta(7, 0) = Vdot(AA2, j03); beta(8, 0) = Vdot(AA3, j03); // Transformation matrix, function of fiber angle const ChMatrixNM<double, 6, 6>& T0 = m_element‐ >GetLayer(m_kl).Get_T0(); // Determinant of the initial position vector gradient at the element center double detJ0C = m_element‐>GetLayer(m_kl).Get_detJ0C(); // Enhanced Assumed Strain ChMatrixNM<double, 6, 5> G = T0 * M * (detJ0C / detJ0); ChMatrixNM<double, 6, 1> strain_EAS = G * (*m_alpha_eas); ChMatrixNM<double, 8, 1> ddNx; ChMatrixNM<double, 8, 1> ddNy; ChMatrixNM<double, 8, 1> ddNz; ddNx.MatrMultiplyT(m_element‐>m_ddT, Nx); ddNy.MatrMultiplyT(m_element‐>m_ddT, Ny); ddNz.MatrMultiplyT(m_element‐>m_ddT, Nz); Transform orthogonal system with direction of fiber (orthotropic material) Calculate beta coefficients Constant transformation to account for initial configuration EAS calculation

slide-26
SLIDE 26

26 ChMatrixNM<double, 8, 1> ddNx; ChMatrixNM<double, 8, 1> ddNy; ChMatrixNM<double, 8, 1> ddNz; ddNx.MatrMultiplyT(m_element‐>m_ddT, Nx); ddNy.MatrMultiplyT(m_element‐>m_ddT, Ny); ddNz.MatrMultiplyT(m_element‐>m_ddT, Nz); ChMatrixNM<double, 8, 1> d0d0Nx; ChMatrixNM<double, 8, 1> d0d0Ny; ChMatrixNM<double, 8, 1> d0d0Nz; d0d0Nx.MatrMultiplyT(m_element‐>m_d0d0T, Nx); d0d0Ny.MatrMultiplyT(m_element‐>m_d0d0T, Ny); d0d0Nz.MatrMultiplyT(m_element‐>m_d0d0T, Nz); Derivatives w.r.t. coordinates

  • 3. Internal forces
slide-27
SLIDE 27

27

// Strain component ChMatrixNM<double, 6, 1> strain_til; strain_til(0, 0) = 0.5 * ((Nx * ddNx)(0, 0) ‐ (Nx * d0d0Nx)(0, 0)); strain_til(1, 0) = 0.5 * ((Ny * ddNy)(0, 0) ‐ (Ny * d0d0Ny)(0, 0)); strain_til(2, 0) = (Nx * ddNy)(0, 0) ‐ (Nx * d0d0Ny)(0, 0); strain_til(3, 0) = N(0, 0) * m_element‐>m_strainANS(0, 0) + N(0, 2) * m_element‐>m_strainANS(1, 0) + N(0, 4) * m_element‐>m_strainANS(2, 0) + N(0, 6) * m_element‐>m_strainANS(3, 0); strain_til(4, 0) = S_ANS(0, 2) * m_element‐>m_strainANS(6, 0) + S_ANS(0, 3) * m_element‐>m_strainANS(7, 0); strain_til(5, 0) = S_ANS(0, 0) * m_element‐>m_strainANS(4, 0) + S_ANS(0, 1) * m_element‐>m_strainANS(5, 0); // For orthotropic material ChMatrixNM<double, 6, 1> strain; strain(0, 0) = strain_til(0, 0) * beta(0) * beta(0) + strain_til(1, 0) * beta(3) * beta(3) + strain_til(2, 0) * beta(0) * beta(3) + strain_til(3, 0) * beta(6) * beta(6) + strain_til(4, 0) * beta(0) * beta(6) + strain_til(5, 0) * beta(3) * beta(6); . . . strain(5, 0) = strain_til(0, 0) * 2.0 * beta(1) * beta(2) + strain_til(1, 0) * 2.0 * beta(4) * beta(5) + strain_til(2, 0) * (beta(2) * beta(4) + beta(1) * beta(5)) + strain_til(3, 0) * 2.0 * beta(7) * beta(8) + strain_til(4, 0) * (beta(2) * beta(7) + beta(1) * beta(8)) + strain_til(5, 0) * (beta(5) * beta(7) + beta(4) * beta(8));

Calculation of covariant compatible strains ANS: Transverse shear and thickness strain Beta matrix/vector multiplication

  • 3. Internal forces
slide-28
SLIDE 28

28 // Strain derivative component ChMatrixNM<double, 6, 24> strainD_til; ChMatrixNM<double, 1, 24> tempB; ChMatrixNM<double, 1, 24> tempBB; ChMatrixNM<double, 1, 3> tempB3; ChMatrixNM<double, 1, 3> tempB31; strainD_til.Reset(); tempB3.MatrMultiply(Nx, m_element‐>m_d); for (int i = 0; i < 8; i++) { for (int j = 0; j < 3; j++) { tempB(0, i * 3 + j) = tempB3(0, j) * Nx(0, i); } } strainD_til.PasteClippedMatrix(&tempB, 0, 0, 1, 24, 0, 0); tempB3.MatrMultiply(Ny, m_element‐>m_d); for (int i = 0; i < 8; i++) { for (int j = 0; j < 3; j++) { tempB(0, i * 3 + j) = tempB3(0, j) * Ny(0, i); } } strainD_til.PasteClippedMatrix(&tempB, 0, 0, 1, 24, 1, 0); tempB31.MatrMultiply(Nx, m_element‐>m_d); for (int i = 0; i < 8; i++) { for (int j = 0; j < 3; j++) { tempB(0, i * 3 + j) = tempB3(0, j) * Nx(0, i) + tempB31(0, j) * Ny(0, i); } } . . . strainD_til.PasteClippedMatrix(&tempB, 0, 0, 1, 24, 5, 0);

  • 3. Internal forces

// For orthotropic material ChMatrixNM<double, 6, 24> strainD; // Derivative

  • f the strains w.r.t. the coordinates. Includes
  • rthotropy

for (int ii = 0; ii < 24; ii++) { strainD(0, ii) = strainD_til(0, ii) * beta(0) * beta(0) + strainD_til(1, ii) * beta(3) * beta(3) + strainD_til(2, ii) * beta(0) * beta(3) + strainD_til(3, ii) * beta(6) * beta(6) + strainD_til(4, ii) * beta(0) * beta(6) + strainD_til(5, ii) * beta(3) * beta(6); . . . strainD(5, ii) = strainD_til(0, ii) * 2.0 * beta(1) * beta(2) + strainD_til(1, ii) * 2.0 * beta(4) * beta(5) + strainD_til(2, ii) * (beta(2) * beta(4) + beta(1) * beta(5)) + strainD_til(3, ii) * 2.0 * beta(7) * beta(8) + strainD_til(4, ii) * (beta(2) * beta(7) + beta(1) * beta(8)) + strainD_til(5, ii) * (beta(5) * beta(7) + beta(4) * beta(8)); } Derivative of covariant strains w.r.t. coordinates: Needed for generalized forces and Jacobians Calculations Pre- and post- multiplication by beta vector: Derivate of engineering strain

slide-29
SLIDE 29
  • 3. Internal forces

29 // Enhanced Assumed Strain 2nd strain += strain_EAS; // Strain time derivative for structural damping ChMatrixNM<double, 6, 1> DEPS; DEPS.Reset(); for (int ii = 0; ii < 24; ii++) { DEPS(0, 0) = DEPS(0, 0) + strainD(0, ii) * m_element‐>m_d_dt(ii, 0); DEPS(1, 0) = DEPS(1, 0) + strainD(1, ii) * m_element‐>m_d_dt(ii, 0); DEPS(2, 0) = DEPS(2, 0) + strainD(2, ii) * m_element‐>m_d_dt(ii, 0); DEPS(3, 0) = DEPS(3, 0) + strainD(3, ii) * m_element‐>m_d_dt(ii, 0); DEPS(4, 0) = DEPS(4, 0) + strainD(4, ii) * m_element‐>m_d_dt(ii, 0); DEPS(5, 0) = DEPS(5, 0) + strainD(5, ii) * m_element‐>m_d_dt(ii, 0); } // Add structural damping strain += DEPS * m_element‐>m_Alpha; // Matrix of elastic coefficients: the input assumes the material *could* be orthotropic const ChMatrixNM<double, 6, 6>& E_eps = m_element‐>GetLayer(m_kl).GetMaterial()‐>Get_E_eps(); // Internal force calculation ChMatrixNM<double, 24, 6> tempC; tempC.MatrTMultiply(strainD, E_eps); ChMatrixNM<double, 24, 1> Fint = (tempC * strain) * (detJ0 * m_element‐>m_GaussScaling); // EAS terms ChMatrixNM<double, 5, 6> temp56; temp56.MatrTMultiply(G, E_eps); ChMatrixNM<double, 5, 1> HE = (temp56 * strain) * (detJ0 * m_element‐>m_GaussScaling); // EAS residual ChMatrixNM<double, 5, 5> KALPHA = (temp56 * G) * (detJ0 * m_element‐>m_GaussScaling); // EAS Jacobian /// Total result vector result.PasteClippedMatrix(&Fint, 0, 0, 24, 1, 0, 0); Sum enhanced assumed strain Add structural damping (not covered, but straightforward) Generalized force per Gauss point EAS generalized force and Jacobian

slide-30
SLIDE 30

30

  • 3. Internal forces
slide-31
SLIDE 31
  • 3. Internal forces

31

slide-32
SLIDE 32
  • 4. Applications

32

slide-33
SLIDE 33

33