From 1494786a1721cb07878c559dc8fb434fcc4474f9 Mon Sep 17 00:00:00 2001 From: Winterless Date: Thu, 27 Jun 2024 12:25:43 +0200 Subject: [PATCH] Add global F, F_bending, and transformation matrix from initial local coordinates to current local coordinates to simplify the calculation. --- .../thin_structure_dynamics.cpp | 25 +++++++++++-------- .../solid_dynamics/thin_structure_dynamics.h | 12 +++------ 2 files changed, 18 insertions(+), 19 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 index bd96f3c675..f8337a70ae 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp @@ -108,6 +108,8 @@ ShellStressRelaxationFirstHalf:: global_moment_(*particles_->registerSharedVariable("GlobalMoment")), mid_surface_cauchy_stress_(*particles_->registerSharedVariable("MidSurfaceCauchyStress")), global_shear_stress_(*particles_->registerSharedVariable("GlobalShearStress")), + global_F_(*particles_->registerSharedVariable("GlobalDeformationGradient")), + global_F_bending_(*particles_->registerSharedVariable("GlobalBendingDeformationGradient")), E0_(elastic_solid_.YoungsModulus()), G0_(elastic_solid_.ShearModulus()), nu_(elastic_solid_.PoissonRatio()), @@ -144,13 +146,18 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) F_[index_i] += dF_dt_[index_i] * dt * 0.5; F_bending_[index_i] += dF_bending_dt_[index_i] * dt * 0.5; + global_F_[index_i] = transformation_matrix0_[index_i].transpose() * F_[index_i] * transformation_matrix0_[index_i]; + global_F_bending_[index_i] = transformation_matrix0_[index_i].transpose() * F_bending_[index_i] * transformation_matrix0_[index_i]; + Real J = F_[index_i].determinant(); - Matd inverse_F = F_[index_i].inverse(); + Matd inverse_transpose_global_F = global_F_[index_i].inverse().transpose(); rho_[index_i] = rho0_ / J; /** Get transformation matrix from global coordinates to current local coordinates. */ Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]); + /** Get transformation matrix from initial local coordinates to current local coordinates. */ + Matd transformation_matrix_0_to_current = current_transformation_matrix * transformation_matrix0_[index_i].transpose(); Matd resultant_stress = Matd::Zero(); Matd resultant_moment = Matd::Zero(); @@ -161,9 +168,9 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) Matd F_gaussian_point = F_[index_i] + gaussian_point_[i] * F_bending_[index_i] * thickness_[index_i] * 0.5; Matd dF_gaussian_point_dt = dF_dt_[index_i] + gaussian_point_[i] * dF_bending_dt_[index_i] * thickness_[index_i] * 0.5; Matd inverse_F_gaussian_point = F_gaussian_point.inverse(); - Matd current_local_almansi_strain = current_transformation_matrix * transformation_matrix0_[index_i].transpose() * 0.5 * + Matd current_local_almansi_strain = transformation_matrix_0_to_current * 0.5 * (Matd::Identity() - inverse_F_gaussian_point.transpose() * inverse_F_gaussian_point) * - transformation_matrix0_[index_i] * current_transformation_matrix.transpose(); + transformation_matrix_0_to_current.transpose(); /** correct Almansi strain tensor according to plane stress problem. */ current_local_almansi_strain = getCorrectedAlmansiStrain(current_local_almansi_strain, nu_); @@ -171,9 +178,9 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) /** correct out-plane numerical damping. */ numerical_damping_scaling_matrix_(Dimensions - 1, Dimensions - 1) = thickness_[i] < smoothing_length_ ? thickness_[i] : smoothing_length_; Matd cauchy_stress = elastic_solid_.StressCauchy(current_local_almansi_strain, index_i) + - current_transformation_matrix * transformation_matrix0_[index_i].transpose() * F_gaussian_point * + transformation_matrix_0_to_current * F_gaussian_point * elastic_solid_.NumericalDampingRightCauchy(F_gaussian_point, dF_gaussian_point_dt, numerical_damping_scaling_matrix_, index_i) * - F_gaussian_point.transpose() * transformation_matrix0_[index_i] * current_transformation_matrix.transpose() / F_gaussian_point.determinant(); + F_gaussian_point.transpose() * transformation_matrix_0_to_current.transpose() / F_gaussian_point.determinant(); /** Impose modeling assumptions. */ cauchy_stress.col(Dimensions - 1) *= shear_correction_factor_; @@ -199,13 +206,9 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) /** stress and moment in global coordinates for pair interaction */ global_stress_[index_i] = J * current_transformation_matrix.transpose() * - resultant_stress * current_transformation_matrix * - transformation_matrix0_[index_i].transpose() * - inverse_F.transpose() * transformation_matrix0_[index_i]; + resultant_stress * current_transformation_matrix * inverse_transpose_global_F; global_moment_[index_i] = J * current_transformation_matrix.transpose() * - resultant_moment * current_transformation_matrix * - transformation_matrix0_[index_i].transpose() * inverse_F.transpose() * - transformation_matrix0_[index_i]; + resultant_moment * current_transformation_matrix * inverse_transpose_global_F; global_shear_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_shear_stress; } //=================================================================================================// 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 8f3bd3db5f..e4ed80ab10 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h @@ -205,20 +205,15 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation Vecd e_ij = inner_neighborhood.e_ij_[n]; Real r_ij = inner_neighborhood.r_ij_[n]; Real weight = inner_neighborhood.W_ij_[n] * inv_W0_; - Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i], - transformation_matrix0_[index_i].transpose() * F_[index_i] * transformation_matrix0_[index_i], - pos_[index_j], - transformation_matrix0_[index_i].transpose() * F_[index_j] * transformation_matrix0_[index_i]); + Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i], global_F_[index_i], pos_[index_j], global_F_[index_j]); Real limiter_pos = SMIN(2.0 * pos_jump.norm() / r_ij, 1.0); force += mass_[index_i] * hourglass_control_factor_ * weight * G0_ * pos_jump * Dimensions * inner_neighborhood.dW_ij_[n] * Vol_[index_j] * limiter_pos; 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_matrix0_[index_i].transpose() * F_bending_[index_i] * transformation_matrix0_[index_i], - pseudo_n_variation_j, - transformation_matrix0_[index_j].transpose() * F_bending_[index_j] * transformation_matrix0_[index_j]); + Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_variation_i, global_F_bending_[index_i], + pseudo_n_variation_j, global_F_bending_[index_j]); 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_ij_[n] * Vol_[index_j] * pow(thickness_[index_i], 2) * limiter_pseudo_n; @@ -246,6 +241,7 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation StdLargeVec &rho_, &mass_; StdLargeVec &global_stress_, &global_moment_, &mid_surface_cauchy_stress_; StdLargeVec &global_shear_stress_; + StdLargeVec &global_F_, &global_F_bending_; Real E0_, G0_, nu_, hourglass_control_factor_; bool hourglass_control_; const Real inv_W0_ = 1.0 / sph_body_.sph_adaptation_->getKernel()->W0(ZeroVecd);