Merge pull request #240 from DongWuTUM/shell_singularity_test

Shell test for singularity in shell model

Author: Xiangyu-Hu

src/for_2D_build/particles/solid_particles_supplementary.cpp

@@ -51,4 +51,16 @@ namespace SPH
 			sqrt(sigmaxx * sigmaxx + sigmayy * sigmayy - sigmaxx * sigmayy + 3.0 * sigmaxy * sigmaxy);
 	}
 	//=============================================================================================//
+	void MidSurfaceVonMisesStressofShells::update(size_t index_i, Real dt)
+	{
+		Matd sigma = mid_surface_cauchy_stress_[index_i];
+
+		Real sigmaxx = sigma(0, 0);
+		Real sigmayy = sigma(1, 1);
+		Real sigmaxy = sigma(0, 1);
+
+		derived_variable_[index_i] =
+			sqrt(sigmaxx * sigmaxx + sigmayy * sigmayy - sigmaxx * sigmayy + 3.0 * sigmaxy * sigmaxy);
+	}
+	//=============================================================================================//
 }

src/for_3D_build/particles/solid_particles_supplementary.cpp

@@ -56,5 +56,22 @@ namespace SPH
 				 sigmaxx * sigmayy - sigmaxx * sigmazz - sigmayy * sigmazz +
 				 3.0 * (sigmaxy * sigmaxy + sigmaxz * sigmaxz + sigmayz * sigmayz));
 	}
+	//=============================================================================================//
+	void MidSurfaceVonMisesStressofShells::update(size_t index_i, Real dt)
+	{
+		Matd sigma = mid_surface_cauchy_stress_[index_i];
+
+		Real sigmaxx = sigma(0, 0);
+		Real sigmayy = sigma(1, 1);
+		Real sigmazz = sigma(2, 2);
+		Real sigmaxy = sigma(0, 1);
+		Real sigmaxz = sigma(0, 2);
+		Real sigmayz = sigma(1, 2);
+
+		derived_variable_[index_i] =
+			sqrt(sigmaxx * sigmaxx + sigmayy * sigmayy + sigmazz * sigmazz -
+				sigmaxx * sigmayy - sigmaxx * sigmazz - sigmayy * sigmazz +
+				3.0 * (sigmaxy * sigmaxy + sigmaxz * sigmaxz + sigmayz * sigmayz));
+	}
 	//=================================================================================================//
 }

src/shared/materials/elastic_solid.cpp

@@ -15,22 +15,6 @@ namespace SPH
 		cs0_ = sqrt(G0_ / rho0_);
 	};
 	//=================================================================================================//
-	Matd ElasticSolid::NumericalDampingRightCauchy(
-		Matd &F, Matd &dF_dt, Real smoothing_length, size_t particle_index_i)
-	{
-		Matd strain_rate = 0.5 * (dF_dt.transpose() * F + F.transpose() * dF_dt);
-		Matd normal_rate = getDiagonal(strain_rate);
-		return 0.5 * rho0_ * (cs0_ * (strain_rate - normal_rate) + c0_ * normal_rate) * smoothing_length;
-	}
-	//=================================================================================================//
-	Matd ElasticSolid::NumericalDampingLeftCauchy(
-		Matd &F, Matd &dF_dt, Real smoothing_length, size_t particle_index_i)
-	{
-		Matd strain_rate = 0.5 * (dF_dt * F.transpose() + F * dF_dt.transpose());
-		Matd normal_rate = getDiagonal(strain_rate);
-		return 0.5 * rho0_ * (cs0_ * (strain_rate - normal_rate) + c0_ * normal_rate) * smoothing_length;
-	}
-	//=================================================================================================//
 	Real ElasticSolid::PairNumericalDamping(Real dE_dt_ij, Real smoothing_length)
 	{
 		return 0.5 * rho0_ * c0_ * dE_dt_ij * smoothing_length;

src/shared/materials/elastic_solid.h

@@ -75,9 +75,21 @@ class ElasticSolid : public Solid
     /** Cauchy stress through Eulerian Almansi strain tensor. */
     virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) = 0;
     /** Numerical damping stress using right Cauchy tensor. */
-    virtual Matd NumericalDampingRightCauchy(Matd &deformation, Matd &deformation_rate, Real smoothing_length, size_t particle_index_i);
-    /** Numerical damping stress using left Cauchy tensor. */
-    virtual Matd NumericalDampingLeftCauchy(Matd &deformation, Matd &deformation_rate, Real smoothing_length, size_t particle_index_i);
+		template <typename ScalingType>
+		Matd NumericalDampingRightCauchy(const Matd &deformation, const Matd &deformation_rate, const ScalingType &scaling, size_t particle_index_i)
+		{
+			Matd strain_rate = 0.5 * (deformation_rate.transpose() * deformation + deformation.transpose() * deformation_rate);
+			Matd normal_rate = getDiagonal(strain_rate);
+			return 0.5 * rho0_ * (cs0_ * (strain_rate - normal_rate) + c0_ * normal_rate) * scaling;
+		}
+		/** Numerical damping stress using left Cauchy tensor. */
+		template <typename ScalingType>
+		Matd NumericalDampingLeftCauchy(const Matd &deformation, const Matd &deformation_rate, const ScalingType& scaling, size_t particle_index_i)
+		{
+			Matd strain_rate = 0.5 * (deformation_rate * deformation.transpose() + deformation * deformation_rate.transpose());
+			Matd normal_rate = getDiagonal(strain_rate);
+			return 0.5 * rho0_ * (cs0_ * (strain_rate - normal_rate) + c0_ * normal_rate) * scaling;
+		}
     /** Numerical damping is computed between particles i and j */
     virtual Real PairNumericalDamping(Real dE_dt_ij, Real smoothing_length);
 

src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp

@@ -72,6 +72,8 @@ namespace SPH
 			  elastic_solid_(particles_->elastic_solid_),
 			  global_stress_(particles_->global_stress_),
 			  global_moment_(particles_->global_moment_),
+			  mid_surface_cauchy_stress_(particles_->mid_surface_cauchy_stress_),
+			  numerical_damping_scaling_(particles_->numerical_damping_scaling_),
 			  global_shear_stress_(particles_->global_shear_stress_),
 			  n_(particles_->n_),
 			  rho0_(elastic_solid_.ReferenceDensity()),
@@ -98,15 +100,8 @@ namespace SPH
 				gaussian_point_ = three_gaussian_points_;
 				gaussian_weight_ = three_gaussian_weights_;
 			}
-			/** Define the factor of hourglass control algorithm according to the dimension. */
-			if (Dimensions == 2)
-			{
-				hourglass_control_factor_ = 0.05;
-			}
-			else
-			{
-				hourglass_control_factor_ = 0.01;
-			}
+			/** Define the factor of hourglass control algorithm. */
+			hourglass_control_factor_ = 0.005;
 		}
 		//=================================================================================================//
 		void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)
@@ -120,12 +115,16 @@ namespace SPH
 
 			F_[index_i] += dF_dt_[index_i] * dt * 0.5;
 			F_bending_[index_i] += dF_bending_dt_[index_i] * dt * 0.5;
-			rho_[index_i] = rho0_ / F_[index_i].determinant();
+
+			Real J = F_[index_i].determinant();
+			Matd inverse_F = F_[index_i].inverse();
+
+			rho_[index_i] = rho0_ / J;
 
 			/** Calculate the current normal direction of mid-surface. */
 			n_[index_i] = transformation_matrix_[index_i].transpose() * getNormalFromDeformationGradientTensor(F_[index_i]);
 			/** Get transformation matrix from global coordinates to current local coordinates. */
-			Matd current_transformation_matrix = getTransformationMatrix(n_[index_i]);
+            Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]);
 
 			Matd resultant_stress = Matd::Zero();
 			Matd resultant_moment = Matd::Zero();
@@ -144,39 +143,40 @@ namespace SPH
 				current_local_almansi_strain = getCorrectedAlmansiStrain(current_local_almansi_strain, nu_);
 
 				/** correct out-plane numerical damping. */
-				Matd numerical_damping = current_transformation_matrix * transformation_matrix_[index_i].transpose() * F_gaussian_point *
-										 elastic_solid_.NumericalDampingRightCauchy(F_gaussian_point, dF_gaussian_point_dt, smoothing_length_, index_i) * F_gaussian_point.transpose() * transformation_matrix_[index_i] * current_transformation_matrix.transpose() / F_gaussian_point.determinant();
-				numerical_damping.col(Dimensions - 1) *= thickness_[index_i] / smoothing_length_;
-
-				Matd cauchy_stress = elastic_solid_.StressCauchy(current_local_almansi_strain, F_gaussian_point, index_i) + numerical_damping;
+				Matd cauchy_stress = elastic_solid_.StressCauchy(current_local_almansi_strain, F_gaussian_point, index_i) 
+									 + current_transformation_matrix * transformation_matrix_[index_i].transpose() * F_gaussian_point *
+									 elastic_solid_.NumericalDampingRightCauchy(F_gaussian_point, dF_gaussian_point_dt, numerical_damping_scaling_[index_i], index_i)
+									 * F_gaussian_point.transpose() * transformation_matrix_[index_i] * current_transformation_matrix.transpose()
+									 / F_gaussian_point.determinant();
 
 				/** Impose modeling assumptions. */
 				cauchy_stress.col(Dimensions - 1) *= shear_correction_factor_;
 				cauchy_stress.row(Dimensions - 1) *= shear_correction_factor_;
 				cauchy_stress(Dimensions - 1, Dimensions - 1) = 0.0;
+				
+				if (i == 0)
+				{
+					mid_surface_cauchy_stress_[index_i] = cauchy_stress;
+				}
 
-				Matd stress_PK2_gaussian_point = F_gaussian_point.determinant() * F_gaussian_point.inverse() * transformation_matrix_[index_i] *
-												 current_transformation_matrix.transpose() * cauchy_stress * current_transformation_matrix * transformation_matrix_[index_i].transpose() * F_gaussian_point.inverse().transpose();
-
-				Vecd shear_stress_PK2_gaussian_point = -stress_PK2_gaussian_point.col(Dimensions - 1);
-				Matd moment_PK2_gaussian_point = stress_PK2_gaussian_point * gaussian_point_[i] * thickness_[index_i] * 0.5;
-
+				/** Integrate Cauchy stress along thickness. */
 				resultant_stress +=
-					0.5 * thickness_[index_i] * gaussian_weight_[i] * F_gaussian_point * stress_PK2_gaussian_point;
+					0.5 * thickness_[index_i] * gaussian_weight_[i] * cauchy_stress;
 				resultant_moment +=
-					0.5 * thickness_[index_i] * gaussian_weight_[i] * F_gaussian_point * moment_PK2_gaussian_point;
-				resultant_shear_stress +=
-					0.5 * thickness_[index_i] * gaussian_weight_[i] * F_gaussian_point * shear_stress_PK2_gaussian_point;
+					0.5 * thickness_[index_i] * gaussian_weight_[i] * (cauchy_stress * gaussian_point_[i] * thickness_[index_i] * 0.5);
+				resultant_shear_stress -=
+					0.5 * thickness_[index_i] * gaussian_weight_[i] * cauchy_stress.col(Dimensions - 1);
+
+				resultant_stress.col(Dimensions - 1) = Vecd::Zero();
+				resultant_moment.col(Dimensions - 1) = Vecd::Zero();
 			}
-			/** Only one (for 2D) or two (for 3D) angular momentum equations left. */
-			resultant_moment.col(Dimensions - 1) = Vecd::Zero();
-			resultant_moment.row(Dimensions - 1) = Vecd::Zero().transpose();
-			resultant_shear_stress[Dimensions - 1] = 0.0;
 
 			/** stress and moment in global coordinates for pair interaction */
-			global_stress_[index_i] = transformation_matrix_[index_i].transpose() * resultant_stress * transformation_matrix_[index_i];
-			global_moment_[index_i] = transformation_matrix_[index_i].transpose() * resultant_moment * transformation_matrix_[index_i];
-			global_shear_stress_[index_i] = transformation_matrix_[index_i].transpose() * resultant_shear_stress;
+			global_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_stress * current_transformation_matrix
+				* transformation_matrix_[index_i].transpose() * inverse_F.transpose() * transformation_matrix_[index_i];
+			global_moment_[index_i] = J * current_transformation_matrix.transpose() * resultant_moment * current_transformation_matrix
+				* transformation_matrix_[index_i].transpose() * inverse_F.transpose() * transformation_matrix_[index_i];
+			global_shear_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_shear_stress;
 		}
 		//=================================================================================================//
 		void ShellStressRelaxationFirstHalf::update(size_t index_i, Real dt)
@@ -190,7 +190,7 @@ namespace SPH
 			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]);
+                                     getVectorChangeRateAfterThinStructureRotation(local_pseudo_n_0, rotation_[index_i], angular_vel_[index_i]);
 			pseudo_n_[index_i] += dpseudo_n_dt_[index_i] * dt * 0.5;
 		}
 		//=================================================================================================//

src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h

@@ -208,7 +208,7 @@ namespace SPH
 															  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 * E0_ * pos_jump * dim_inv_r_ij *
+						acceleration += hourglass_control_factor_ * weight * G0_ * pos_jump * dim_inv_r_ij *
 										inner_neighborhood.dW_ijV_j_[n] * thickness_[index_i];
 
 						Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_[index_i] - n0_[index_i],
@@ -236,7 +236,7 @@ namespace SPH
 
 		protected:
 			ElasticSolid &elastic_solid_;
-			StdLargeVec<Matd> &global_stress_, &global_moment_;
+			StdLargeVec<Matd> &global_stress_, &global_moment_, &mid_surface_cauchy_stress_, &numerical_damping_scaling_;
 			StdLargeVec<Vecd> &global_shear_stress_, &n_;
 			Real rho0_, inv_rho0_;
 			Real smoothing_length_, E0_, G0_, nu_, hourglass_control_factor_;

src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp

@@ -44,85 +44,64 @@ namespace SPH
 		//=================================================================================================//
 		Vec2d getVectorChangeRateAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles, const Vec2d &angular_vel)
 		{
-			/**The derivative of the rotation matrix. */
-			Real sin_angle = sin(rotation_angles[0]);
-			Real cos_angle = cos(rotation_angles[0]);
-			Mat2d drotation_matrix_dt{
-				{-sin_angle * angular_vel[0],  cos_angle * angular_vel[0]},
-				{-cos_angle * angular_vel[0], -sin_angle * angular_vel[0]},
-			};
-
-			return drotation_matrix_dt * initial_vector;
+			return Vec2d(cos(rotation_angles[0]) * angular_vel[0], -sin(rotation_angles[0]) * angular_vel[0]);
 		}
 		//=================================================================================================//
 		Vec3d getVectorChangeRateAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles, const Vec3d &angular_vel)
 		{
-			/**The rotation matrix about the 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]);
+			Real sin_rotation_0 = sin(rotation_angles[0]);
+			Real cos_rotation_0 = cos(rotation_angles[0]);
 
-			Real angular_vel_x = angular_vel[0];
-			Real angular_vel_y = angular_vel[1];
+			Real sin_rotation_1 = sin(rotation_angles[1]);
+			Real cos_rotation_1 = cos(rotation_angles[1]);
 
-			Mat3d rotation_matrix_x{
-				{1.0, 0.0, 0.0},
-				{0.0, cos_angle_x, -sin_angle_x},
-				{0.0, sin_angle_x, cos_angle_x},
-			};
-			/**The rotation matrix about the Y-axis. */
-			Mat3d rotation_matrix_y{
-				{cos_angle_y, 0.0, sin_angle_y},
-				{0.0, 1.0, 0.0},
-				{-sin_angle_y, 0.0, cos_angle_y},
-			};
-			/**The derivative of the rotation matrix of the X-axis. */
-			Mat3d drotation_matrix_x_dt{
-				{0.0, 0.0, 0.0},
-				{0.0,-sin_angle_x * angular_vel_x, -cos_angle_x * angular_vel_x},
-				{0.0, cos_angle_x * angular_vel_x, -sin_angle_x * angular_vel_x},
-			};
-			/**The derivative of the rotation matrix of the Y-axis. */
-			Mat3d drotation_matrix_y_dt{
-				{-sin_angle_y * angular_vel_y, 0.0, cos_angle_y * angular_vel_y},
-				{0.0, 0.0, 0.0},
-				{-cos_angle_y * angular_vel_y, 0.0,-sin_angle_y * angular_vel_y},
-			};
+			Real dpseudo_n_dt_0 = -sin_rotation_0 * sin_rotation_1 * angular_vel[0] + cos_rotation_0 * cos_rotation_1 * angular_vel[1];
+			Real dpseudo_n_dt_1 = -cos_rotation_0 * angular_vel[0];
+			Real dpseudo_n_dt_2 = -sin_rotation_0 * cos_rotation_1 * angular_vel[0] - cos_rotation_0 * sin_rotation_1 * angular_vel[1];
 
-			return (drotation_matrix_y_dt * rotation_matrix_x + rotation_matrix_y * drotation_matrix_x_dt)* initial_vector;
+			return Vec3d(dpseudo_n_dt_0, dpseudo_n_dt_1, dpseudo_n_dt_2);
 		}
 		//=================================================================================================//
 		Vec2d getRotationFromPseudoNormalForFiniteDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt)
-		{
-			Vec2d dangular_vel_dt = Vec2d::Zero();
-			dangular_vel_dt[0] = -(dpseudo_n_d2t[0] + sin(rotation[0]) * pow(angular_vel[0], 2))
-								 / (2 * sin(rotation[0]) * angular_vel[0] * dt - cos(rotation[0]));
-			return dangular_vel_dt;
+        {
+			Real cos_rotation_0 = cos(rotation[0]);
+			Real sin_rotation_0 = sin(rotation[0]);
+
+			Real angle_vel_dt_0 = cos_rotation_0 * (dpseudo_n_d2t[0] + sin_rotation_0 * angular_vel[0] * angular_vel[0]) 
+								  - sin_rotation_0 * (dpseudo_n_d2t[1] + cos_rotation_0 * angular_vel[0] * angular_vel[0]);
+
+            return Vec2d(angle_vel_dt_0, 0.0);
 		}
 		//=================================================================================================//
 		Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt)
 		{
-			Real sin_rotation_x = sin(rotation[0]);
-			Real cos_rotation_x = cos(rotation[0]);
-			Real sin_rotation_y = sin(rotation[1]);
-			Real cos_rotation_y = cos(rotation[1]);
+			Real sin_rotation_0 = sin(rotation[0]);
+			Real cos_rotation_0 = cos(rotation[0]);
+			Real sin_rotation_1 = sin(rotation[1]);
+			Real cos_rotation_1 = cos(rotation[1]);
 
-			Real angle_vel_dt_x = (dpseudo_n_d2t[1] - sin_rotation_x * pow(angular_vel[0], 2))
-								 / (2 * sin_rotation_x * angular_vel[0] * dt - cos_rotation_x);
+			Real rotation_0_a = -(dpseudo_n_d2t[2] * cos_rotation_1 + dpseudo_n_d2t[0] * sin_rotation_1
+								  + angular_vel[1] * angular_vel[1] * cos_rotation_0
+								  + angular_vel[0] * angular_vel[0] * cos_rotation_0);
+			Real rotation_0_b = sin_rotation_0 * angular_vel[0] * angular_vel[0] - dpseudo_n_d2t[1];
+			Real angle_vel_dt_0 = sin_rotation_0 * rotation_0_a + cos_rotation_0 * rotation_0_b;
 
-			Real angle_vel_dt_y = (dpseudo_n_d2t[0] + cos_rotation_x * sin_rotation_y
-								  * (pow(angular_vel[0], 2) + pow(angular_vel[1], 2))
-								  + 2 * sin_rotation_x * cos_rotation_y * angular_vel[0] * angular_vel[1]
-								  + (2 * cos_rotation_x * sin_rotation_y * angular_vel[0] * dt
-								  + 2 * sin_rotation_x * cos_rotation_y * angular_vel[1] * dt
-								  + sin_rotation_x * cos_rotation_y) * angle_vel_dt_x)
-								 / (-2 * sin_rotation_x * cos_rotation_y * angular_vel[0] * dt
-									- 2 * cos_rotation_x * sin_rotation_y * angular_vel[1] * dt
-									+ cos_rotation_x * cos_rotation_y);
-									
-			return Vec3d(angle_vel_dt_x, angle_vel_dt_y, 0.0);
+			Real rotation_1_a = dpseudo_n_d2t[0] * cos_rotation_1 - dpseudo_n_d2t[2] * sin_rotation_1
+								+ 2.0 * angular_vel[1] * angular_vel[0] * sin_rotation_0;
+			Real rotation_1_b1 = dpseudo_n_d2t[0] * cos_rotation_0
+								 + angular_vel[1] * angular_vel[1] * cos_rotation_0 * cos_rotation_0 * sin_rotation_1
+								 + angular_vel[0] * angular_vel[0] * sin_rotation_1
+								 - dpseudo_n_d2t[1] * sin_rotation_1 * sin_rotation_0
+								 + 2.0 * angular_vel[1] * angular_vel[0] * cos_rotation_1 * cos_rotation_0 * sin_rotation_0;
+			Real rotation_1_b2 = -(dpseudo_n_d2t[2] * cos_rotation_0
+								   + angular_vel[1] * angular_vel[1] * cos_rotation_1 * cos_rotation_0 * cos_rotation_0
+								   + angular_vel[0] * angular_vel[0] * cos_rotation_1
+								   - dpseudo_n_d2t[1] * cos_rotation_1 * sin_rotation_0
+								   - 2.0 * angular_vel[1] * angular_vel[0] * cos_rotation_0 * sin_rotation_1 * sin_rotation_0);
+			Real angle_vel_dt_1 = rotation_1_a * rotation_1_a * (rotation_1_b1 * cos_rotation_1 + rotation_1_b2 * sin_rotation_1)
+								  / (rotation_1_b1 * rotation_1_b1 + rotation_1_b2 * rotation_1_b2 + Eps);
+					
+			return Vec3d(angle_vel_dt_0, angle_vel_dt_1, 0.0);
 		}
 		//=================================================================================================//
 		Vec2d getRotationFromPseudoNormalForSmallDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt)