From 2a4913fe2aa5ac585b23827b6af2827b10cefde5 Mon Sep 17 00:00:00 2001 From: xipenglyu Date: Wed, 2 Aug 2023 14:35:35 +0800 Subject: [PATCH] add getRotationFromPseudoNormalForSmallDeformation. --- .../slender_structure_dynamics.cpp | 5 ++-- .../slender_structure_dynamics.h | 6 ++--- .../solid_dynamics/slender_structure_math.cpp | 4 +-- .../thin_structure_dynamics.cpp | 18 ++++++++++--- .../solid_dynamics/thin_structure_dynamics.h | 5 ++-- .../solid_dynamics/thin_structure_math.cpp | 27 ++++++++----------- 6 files changed, 36 insertions(+), 29 deletions(-) diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp index c082cb18e1..fca680fb22 100644 --- a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp @@ -231,8 +231,9 @@ void BarStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) void BarStressRelaxationFirstHalf::update(size_t index_i, Real dt) { vel_[index_i] += (acc_prior_[index_i] + acc_[index_i]) * dt; - angular_vel_[index_i] += (dangular_vel_dt_[index_i]) * dt; - angular_b_vel_[index_i] += (dangular_b_vel_dt_[index_i]) * dt; + Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i], pseudo_b_n_[index_i]); + angular_vel_[index_i] += (transformation_matrix_[index_i] * (current_transformation_matrix.transpose() * dangular_vel_dt_[index_i])) * dt; + angular_b_vel_[index_i] += (transformation_matrix_[index_i] * (current_transformation_matrix.transpose() * dangular_b_vel_dt_[index_i])) * dt; } //=================================================================================================// void BarStressRelaxationSecondHalf::initialization(size_t index_i, Real dt) diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h index d9bbd593c4..b673335e1e 100644 --- a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h @@ -232,9 +232,9 @@ class BarStressRelaxationFirstHalf : public BaseBarRelaxation acc_[index_i] = acceleration * inv_rho0_ / (thickness_[index_i] * width_[index_i]); dpseudo_n_d2t_[index_i] = pseudo_normal_acceleration * inv_rho0_ * 12.0 / pow(thickness_[index_i], 4); dpseudo_b_n_d2t_[index_i] = -pseudo_b_normal_acceleration * inv_rho0_ * 12.0 / pow(thickness_[index_i], 4); - - Vecd local_dpseudo_n_d2t = transformation_matrix_[index_i] * dpseudo_n_d2t_[index_i]; - Vecd local_dpseudo_b_n_d2t = transformation_matrix_[index_i] * dpseudo_b_n_d2t_[index_i]; + Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i], pseudo_b_n_[index_i]); + Vecd local_dpseudo_n_d2t = current_transformation_matrix * dpseudo_n_d2t_[index_i]; + Vecd local_dpseudo_b_n_d2t = current_transformation_matrix * dpseudo_b_n_d2t_[index_i]; dangular_b_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation_b( Vec3d(local_dpseudo_b_n_d2t), Vec3d(local_dpseudo_n_d2t), Vec3d(rotation_b_[index_i]), Vec3d(angular_b_vel_[index_i]), dt); dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation( diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp index d54976c6a0..0bf8c0cd73 100644 --- a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp @@ -77,14 +77,14 @@ Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t Vec3d getRotationFromPseudoNormalForSmallDeformation( const Vec3d &dpseudo_b_n_d2t, const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt) { - return Vec3d(0.0, dpseudo_n_d2t[0], 0.0); + return Vec3d(-dpseudo_n_d2t[1] - dpseudo_b_n_d2t[2], dpseudo_n_d2t[0], dpseudo_b_n_d2t[0]); } //=================================================================================================// Vec3d getRotationFromPseudoNormalForSmallDeformation_b( const Vec3d &dpseudo_b_n_d2t, const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt) { - return Vec3d(0.0, 0.0, dpseudo_b_n_d2t[0]); + return Vec3d(-dpseudo_n_d2t[1] - dpseudo_b_n_d2t[2], dpseudo_n_d2t[0], dpseudo_b_n_d2t[0]); } //=================================================================================================// diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp index c618f8b0b1..b5f47bf5ed 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp @@ -177,16 +177,26 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) void ShellStressRelaxationFirstHalf::update(size_t index_i, Real dt) { vel_[index_i] += (acc_prior_[index_i] + acc_[index_i]) * dt; - angular_vel_[index_i] += dangular_vel_dt_[index_i] * dt; + Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]); + angular_vel_[index_i] += (transformation_matrix_[index_i] * (current_transformation_matrix.transpose() * dangular_vel_dt_[index_i])) * dt; } //=================================================================================================// void ShellStressRelaxationSecondHalf::initialization(size_t index_i, Real dt) { pos_[index_i] += vel_[index_i] * dt * 0.5; rotation_[index_i] += angular_vel_[index_i] * dt * 0.5; - dpseudo_n_dt_[index_i] = transformation_matrix_[index_i].transpose() * - getVectorChangeRateAfterThinStructureRotation(local_pseudo_n_0, rotation_[index_i], angular_vel_[index_i]); - pseudo_n_[index_i] += dpseudo_n_dt_[index_i] * dt * 0.5; + Vecd pseudo_n_temp = pseudo_n_[index_i]; + + pseudo_n_[index_i] = transformation_matrix_[index_i].transpose() * getVectorAfterThinStructureRotation(local_pseudo_n_0, rotation_[index_i]); + + if (dt < 1e-10) + { + dpseudo_n_dt_[index_i] = Vecd::Zero(); + } + else + { + dpseudo_n_dt_[index_i] = (pseudo_n_[index_i] - pseudo_n_temp) / (0.5 * dt); + } } //=================================================================================================// void ShellStressRelaxationSecondHalf::update(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 index 61dbfd100b..aa6a4005c1 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h @@ -230,8 +230,9 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation dpseudo_n_d2t_[index_i] = pseudo_normal_acceleration * inv_rho0_ * 12.0 / pow(thickness_[index_i], 3); /** the relation between pseudo-normal and rotations */ - Vecd local_dpseudo_n_d2t = transformation_matrix_[index_i] * dpseudo_n_d2t_[index_i]; - dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForFiniteDeformation(local_dpseudo_n_d2t, rotation_[index_i], angular_vel_[index_i], dt); + Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]); + Vecd local_dpseudo_n_d2t = current_transformation_matrix * dpseudo_n_d2t_[index_i]; + dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation(local_dpseudo_n_d2t, rotation_[index_i], angular_vel_[index_i], dt); }; void update(size_t index_i, Real dt = 0.0); diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp index fd1c84537a..fe380df118 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp @@ -23,22 +23,17 @@ Vec2d getVectorAfterThinStructureRotation(const Vec2d &initial_vector, const Vec Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles) { /**The rotation matrix is the rotation around Y-axis \times rotation around X-axis. */ - Real sin_angle_x = sin(rotation_angles[0]); - Real cos_angle_x = cos(rotation_angles[0]); - - Real sin_angle_y = sin(rotation_angles[1]); - Real cos_angle_y = cos(rotation_angles[1]); - - Mat3d rotation_matrix = Mat3d{ - {cos_angle_y, 0.0, sin_angle_y}, - {0.0, 1.0, 0.0}, - {-sin_angle_y, 0.0, cos_angle_y}, - } * - Mat3d{ - {1.0, 0.0, 0.0}, - {0.0, cos_angle_x, -sin_angle_x}, - {0.0, sin_angle_x, cos_angle_x}, - }; + + Real theta = sqrt(rotation_angles[0] * rotation_angles[0] + rotation_angles[1] * rotation_angles[1] + rotation_angles[2] * rotation_angles[2]); + Mat3d R1 = Mat3d::Zero(); + R1(0, 1) = -rotation_angles[2]; + R1(0, 2) = rotation_angles[1]; + R1(1, 0) = rotation_angles[2]; + R1(1, 2) = -rotation_angles[0]; + R1(2, 0) = -rotation_angles[1]; + R1(2, 1) = rotation_angles[0]; + + Mat3d rotation_matrix = Mat3d::Identity() + sin(theta) / (theta + Eps) * R1 + (1 - cos(theta)) / (theta * theta + Eps) * R1 * R1; return rotation_matrix * initial_vector; }