Modify shell hourglass control algorithm.

Xiangyu-Hu · Jul 23, 2023 · 21785dc · 21785dc
1 parent 416323e
commit 21785dc
Show file tree

Hide file tree

Showing 4 changed files with 11 additions and 30 deletions.
diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp
@@ -101,7 +101,7 @@ ShellStressRelaxationFirstHalf::
         gaussian_weight_ = three_gaussian_weights_;
     }
     /** Define the factor of hourglass control algorithm. */
-    hourglass_control_factor_ = 0.005;
+    hourglass_control_factor_ = 0.002;
 }
 //=================================================================================================//
 void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)

diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h
@@ -202,22 +202,24 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation
             {
                 Vecd e_ij = inner_neighborhood.e_ij_[n];
                 Real r_ij = inner_neighborhood.r_ij_[n];
-                Real dim_inv_r_ij = Dimensions / r_ij;
                 Real weight = inner_neighborhood.W_ij_[n] * inv_W0_;
                 Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i],
                                                       transformation_matrix_[index_i].transpose() * F_[index_i] * transformation_matrix_[index_i],
                                                       pos_[index_j],
                                                       transformation_matrix_[index_i].transpose() * F_[index_j] * transformation_matrix_[index_i]);
-                acceleration += hourglass_control_factor_ * weight * G0_ * pos_jump * dim_inv_r_ij *
-                                inner_neighborhood.dW_ijV_j_[n] * thickness_[index_i];
+                Real limiter_pos = SMIN(2.0 * pos_jump.norm() / r_ij, 1.0);
+                acceleration += hourglass_control_factor_ * weight * G0_ * pos_jump * Dimensions *
+                                inner_neighborhood.dW_ijV_j_[n] * limiter_pos;
 
-                Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_[index_i] - n0_[index_i],
+                Vecd pseudo_n_variation_i = pseudo_n_[index_i] - n0_[index_i];
+                Vecd pseudo_n_variation_j = pseudo_n_[index_j] - n0_[index_j];
+                Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_variation_i,
                                                            transformation_matrix_[index_i].transpose() * F_bending_[index_i] * transformation_matrix_[index_i],
-                                                           pseudo_n_[index_j] - n0_[index_j],
+                                                           pseudo_n_variation_j,
                                                            transformation_matrix_[index_j].transpose() * F_bending_[index_j] * transformation_matrix_[index_j]);
-                Vecd rotation_jump = getRotationJump(pseudo_n_jump, transformation_matrix_[index_i]);
-                pseudo_normal_acceleration += hourglass_control_factor_ * weight * G0_ * rotation_jump * dim_inv_r_ij *
-                                              inner_neighborhood.dW_ijV_j_[n] * pow(thickness_[index_i], 3);
+                Real limiter_pseudo_n = SMIN(2.0 * pseudo_n_jump.norm() / ((pseudo_n_variation_i- pseudo_n_variation_j).norm() + Eps), 1.0);
+                pseudo_normal_acceleration += hourglass_control_factor_ * weight * G0_ * pseudo_n_jump * Dimensions *
+                                              inner_neighborhood.dW_ijV_j_[n] * pow(thickness_[index_i], 2) * limiter_pseudo_n;
             }
 
             acceleration += (global_stress_i + global_stress_[index_j]) * inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n];

diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp
@@ -173,23 +173,6 @@ Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_
     return getWENOStateWithStencilPoints(v1, v2, v3, v4);
 }
 //=================================================================================================//
-Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix)
-{
-    Vec2d local_rotation_jump = Vec2d::Zero();
-    Vec2d local_pseuodo_n_jump = transformation_matrix * pseudo_n_jump;
-    local_rotation_jump[0] = local_pseuodo_n_jump[0];
-    return transformation_matrix.transpose() * local_rotation_jump;
-}
-//=================================================================================================//
-Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix)
-{
-    Vec3d local_rotation_jump = Vec3d::Zero();
-    Vec3d local_pseuodo_n_jump = transformation_matrix * pseudo_n_jump;
-    local_rotation_jump[0] = local_pseuodo_n_jump[0];
-    local_rotation_jump[1] = local_pseuodo_n_jump[1];
-    return transformation_matrix.transpose() * local_rotation_jump;
-}
-//=================================================================================================//
 Mat2d getCorrectedAlmansiStrain(const Mat2d &current_local_almansi_strain, const Real &nu_)
 {
     Mat2d corrected_almansi_strain = current_local_almansi_strain;

diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h
@@ -76,10 +76,6 @@ Vecd getWENOLeftState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i
 Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value,
                        const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value);
 
-/** get the artificial rotation from the pseudo-normal jump. */
-Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix);
-Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix);
-
 /** get the corrected Eulerian Almansi strain tensor according to plane stress problem. */
 Mat2d getCorrectedAlmansiStrain(const Mat2d &current_local_almansi_strain, const Real &nu_);
 Mat3d getCorrectedAlmansiStrain(const Mat3d &current_local_almansi_strain, const Real &nu_);