diff --git a/src/shared/materials/general_continuum.cpp b/src/shared/materials/general_continuum.cpp index 22afa89d4e..d799cdfec8 100644 --- a/src/shared/materials/general_continuum.cpp +++ b/src/shared/materials/general_continuum.cpp @@ -43,16 +43,15 @@ Mat3d PlasticContinuum::ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &st Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose()); Mat3d deviatoric_strain_rate = strain_rate - (1.0 / stress_dimension_) * strain_rate.trace() * Mat3d::Identity(); Mat3d stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate + K_ * strain_rate.trace() * Mat3d::Identity() + stress_tensor * (spin_rate.transpose()) + spin_rate * stress_tensor; - Real stress_tensor_I1 = stress_tensor.trace(); Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity(); Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); - Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor_I1 - k_c_; + Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor.trace() - k_c_; Real lambda_dot_ = 0; Mat3d g = Mat3d::Zero(); if (f >= TinyReal) { Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum(); - //non-associate flow rule + // non-associate flow rule lambda_dot_ = (3.0 * alpha_phi_ * K_ * strain_rate.trace() + (G_ / sqrt(stress_tensor_J2)) * deviatoric_stress_times_strain_rate) / (9.0 * alpha_phi_ * K_ * getDPConstantsA(psi_) + G_); g = lambda_dot_ * (3.0 * K_ * getDPConstantsA(psi_) * Mat3d::Identity() + G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2))); } @@ -64,9 +63,7 @@ Mat3d PlasticContinuum::ReturnMapping(Mat3d &stress_tensor) { Real stress_tensor_I1 = stress_tensor.trace(); if (-alpha_phi_ * stress_tensor_I1 + k_c_ < 0) - { stress_tensor -= (1.0 / stress_dimension_) * (stress_tensor_I1 - k_c_ / alpha_phi_) * Mat3d::Identity(); - } stress_tensor_I1 = stress_tensor.trace(); Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity(); Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); @@ -77,4 +74,46 @@ Mat3d PlasticContinuum::ReturnMapping(Mat3d &stress_tensor) } return stress_tensor; } +//=================================================================================================// +Matd J2Plasticity::ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress, Real &hardening_factor) +{ + Matd strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); + Matd deviatoric_strain_rate = strain_rate - (1.0 / (Real)Dimensions) * strain_rate.trace() * Matd::Identity(); + Matd shear_stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate; + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + Real lambda_dot_ = 0; + Matd g = Matd::Zero(); + if (f > TinyReal) + { + Real deviatoric_stress_times_strain_rate = (shear_stress.cwiseProduct(strain_rate)).sum(); + lambda_dot_ = deviatoric_stress_times_strain_rate / (sqrt(2.0 * stress_tensor_J2) * (1.0 + hardening_modulus_ / (3.0 * G_))); + g = lambda_dot_ * (sqrt(2.0) * G_ * shear_stress / (sqrt(stress_tensor_J2))); + } + return shear_stress_rate_elastic - g; +} +//=================================================================================================// +Matd J2Plasticity::ReturnMappingShearStress(Matd &shear_stress, Real &hardening_factor) +{ + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + Real r = 1.0; + if (f > TinyReal) + r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); + return r * shear_stress; +} +//=================================================================================================// +Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor) +{ + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + return (f > TinyReal) ? (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal) : 1.0; +} +//=================================================================================================// +Real J2Plasticity::HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor) +{ + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + return (f > TinyReal) ? 0.5 * f / (G_ + hardening_modulus_ / 3.0) : 0.0; +} } // namespace SPH \ No newline at end of file diff --git a/src/shared/materials/general_continuum.h b/src/shared/materials/general_continuum.h index 115b4eeba2..681d652324 100644 --- a/src/shared/materials/general_continuum.h +++ b/src/shared/materials/general_continuum.h @@ -22,7 +22,7 @@ * ------------------------------------------------------------------------- */ /** * @file general_continuum.h - * @brief Describe the linear elastic and Drucker-Prager's plastic model + * @brief Describe the linear elastic, J2 plasticity, and Drucker-Prager's plastic model * @author Shuaihao Zhang and Xiangyu Hu */ @@ -43,7 +43,7 @@ class GeneralContinuum : public WeaklyCompressibleFluid Real contact_stiffness_; /* contact-force stiffness related to bulk modulus*/ public: explicit GeneralContinuum(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio) - : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(c0 * c0) + : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(rho0_ * c0 * c0) { material_type_name_ = "GeneralContinuum"; E_ = youngs_modulus; @@ -98,5 +98,31 @@ class PlasticContinuum : public GeneralContinuum virtual GeneralContinuum *ThisObjectPtr() override { return this; }; }; + +class J2Plasticity : public GeneralContinuum +{ + protected: + Real yield_stress_; + Real hardening_modulus_; + const Real sqrt_2_over_3_ = sqrt(2.0 / 3.0); + + public: + explicit J2Plasticity(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio, Real yield_stress, Real hardening_modulus = 0.0) + : GeneralContinuum(rho0, c0, youngs_modulus, poisson_ratio), + yield_stress_(yield_stress), hardening_modulus_(hardening_modulus) + { + material_type_name_ = "J2Plasticity"; + }; + virtual ~J2Plasticity(){}; + + Real YieldStress() { return yield_stress_; }; + Real HardeningModulus() { return hardening_modulus_; }; + + virtual Matd ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress, Real &hardening_factor); + virtual Matd ReturnMappingShearStress(Matd &shear_stress, Real &hardening_factor); + virtual Real ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor); + virtual Real HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor); + virtual J2Plasticity *ThisObjectPtr() override { return this; }; +}; } // namespace SPH #endif // GENERAL_CONTINUUM_H diff --git a/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h b/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h index dbcc5ee06e..328c11e76a 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h +++ b/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h @@ -30,3 +30,4 @@ #include "base_continuum_dynamics.h" #include "continuum_integration.hpp" +#include "continuum_dynamics_variable.h" diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp new file mode 100644 index 0000000000..c084c36ba4 --- /dev/null +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp @@ -0,0 +1,57 @@ +#include "continuum_dynamics_variable.h" + +namespace SPH +{ +namespace continuum_dynamics +{ +//=============================================================================================// +VonMisesStress::VonMisesStress(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "VonMisesStress"), + p_(particles_->getVariableDataByName("Pressure")), + shear_stress_(particles_->getVariableDataByName("ShearStress")) {} +//=============================================================================================// +void VonMisesStress::update(size_t index_i, Real dt) +{ + Matd stress_tensor = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); + derived_variable_[index_i] = getVonMisesStressFromMatrix(stress_tensor); +} +//=============================================================================================// +VonMisesStrain::VonMisesStrain(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "VonMisesStrain"), + strain_tensor_(particles_->getVariableDataByName("StrainTensor")) {} +//=============================================================================================// +void VonMisesStrain::update(size_t index_i, Real dt) +{ + derived_variable_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); +} +//=============================================================================================// +VerticalStress::VerticalStress(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "VerticalStress"), + stress_tensor_3D_(particles_->getVariableDataByName("StressTensor3D")) {} +//=============================================================================================// +void VerticalStress::update(size_t index_i, Real dt) +{ + derived_variable_[index_i] = stress_tensor_3D_[index_i](1, 1); +} +//=============================================================================================// +AccDeviatoricPlasticStrain::AccDeviatoricPlasticStrain(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "AccDeviatoricPlasticStrain"), + plastic_continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), + stress_tensor_3D_(particles_->getVariableDataByName("StressTensor3D")), + strain_tensor_3D_(particles_->getVariableDataByName("StrainTensor3D")), + E_(plastic_continuum_.getYoungsModulus()), nu_(plastic_continuum_.getPoissonRatio()) {} +//=============================================================================================// +void AccDeviatoricPlasticStrain::update(size_t index_i, Real dt) +{ + Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1.0 / 3.0) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity(); + Real hydrostatic_pressure = (1.0 / 3.0) * stress_tensor_3D_[index_i].trace(); + Mat3d elastic_strain_tensor_3D = deviatoric_stress / (2.0 * plastic_continuum_.getShearModulus(E_, nu_)) + + hydrostatic_pressure * Mat3d::Identity() / (9.0 * plastic_continuum_.getBulkModulus(E_, nu_)); + Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D; + Mat3d deviatoric_strain_tensor = plastic_strain_tensor_3D - (1.0 / (Real)Dimensions) * plastic_strain_tensor_3D.trace() * Mat3d::Identity(); + Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum(); + derived_variable_[index_i] = sqrt(sum * 2.0 / 3.0); +} +//=================================================================================================// +} // namespace continuum_dynamics +} // namespace SPH diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h new file mode 100644 index 0000000000..54fd194b37 --- /dev/null +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h @@ -0,0 +1,99 @@ +/* ------------------------------------------------------------------------- * + * SPHinXsys * + * ------------------------------------------------------------------------- * + * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle * + * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * + * physical accurate simulation and aims to model coupled industrial dynamic * + * systems including fluid, solid, multi-body dynamics and beyond with SPH * + * (smoothed particle hydrodynamics), a meshless computational method using * + * particle discretization. * + * * + * SPHinXsys is partially funded by German Research Foundation * + * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * + * HU1527/12-1 and HU1527/12-4. * + * * + * Portions copyright (c) 2017-2023 Technical University of Munich and * + * the authors' affiliations. * + * * + * Licensed under the Apache License, Version 2.0 (the "License"); you may * + * not use this file except in compliance with the License. You may obtain a * + * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * + * * + * ------------------------------------------------------------------------- */ +/** + * @file continuum_dynamics_variable.h + * @brief Here, we define the algorithm classes for computing derived solid dynamics variables. + * @details These variable can be added into variable list for state output. + * @author Shuaihao Zhang and Xiangyu Hu + */ + +#ifndef CONTINUUM_DYNAMICS_VARIABLE_H +#define CONTINUUM_DYNAMICS_VARIABLE_H + +#include "base_general_dynamics.h" +#include "general_continuum.h" + +namespace SPH +{ +namespace continuum_dynamics +{ +/** + * @class VonMisesStress + * @brief computing von_Mises_stress + */ +class VonMisesStress : public BaseDerivedVariable +{ + public: + explicit VonMisesStress(SPHBody &sph_body); + virtual ~VonMisesStress(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + Real *p_; + Matd *shear_stress_; +}; +/** + * @class VonMisesStrain + * @brief computing von_Mises_strain + */ +class VonMisesStrain : public BaseDerivedVariable +{ + public: + explicit VonMisesStrain(SPHBody &sph_body); + virtual ~VonMisesStrain(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + Matd *strain_tensor_; +}; +/** + * @class VerticalStress + */ +class VerticalStress : public BaseDerivedVariable +{ + public: + explicit VerticalStress(SPHBody &sph_body); + virtual ~VerticalStress(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + Mat3d *stress_tensor_3D_; +}; +/** + * @class AccumulatedDeviatoricPlasticStrain + */ +class AccDeviatoricPlasticStrain : public BaseDerivedVariable +{ + public: + explicit AccDeviatoricPlasticStrain(SPHBody &sph_body); + virtual ~AccDeviatoricPlasticStrain(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + PlasticContinuum &plastic_continuum_; + Mat3d *stress_tensor_3D_, *strain_tensor_3D_; + Real E_, nu_; +}; +} // namespace continuum_dynamics +} // namespace SPH +#endif // CONTINUUM_DYNAMICS_VARIABLE_H \ No newline at end of file diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index 6df3bd571c..aac2151fde 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -11,119 +11,152 @@ ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body) vel_(particles_->registerStateVariable("Velocity")), stress_tensor_3D_(particles_->registerStateVariable("StressTensor3D")) {} //=================================================================================================// -ShearAccelerationRelaxation::ShearAccelerationRelaxation(BaseInnerRelation &inner_relation) - : fluid_dynamics::BaseIntegration(inner_relation), - continuum_(DynamicCast(this, particles_->getBaseMaterial())), - G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), - smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), - acc_shear_(particles_->registerStateVariable("AccelerationByShear")) +AcousticTimeStep::AcousticTimeStep(SPHBody &sph_body, Real acousticCFL) + : LocalDynamicsReduce(sph_body), + fluid_(DynamicCast(this, particles_->getBaseMaterial())), + rho_(particles_->getVariableDataByName("Density")), + p_(particles_->getVariableDataByName("Pressure")), + vel_(particles_->getVariableDataByName("Velocity")), + smoothing_length_min_(sph_body.sph_adaptation_->MinimumSmoothingLength()), + acousticCFL_(acousticCFL) {} +//=================================================================================================// +Real AcousticTimeStep::reduce(size_t index_i, Real dt) { - particles_->addVariableToSort("AccelerationByShear"); + return fluid_.getSoundSpeed(p_[index_i], rho_[index_i]) + vel_[index_i].norm(); } //=================================================================================================// -void ShearAccelerationRelaxation::interaction(size_t index_i, Real dt) +Real AcousticTimeStep::outputResult(Real reduced_value) { - Real rho_i = rho_[index_i]; - Vecd acceleration = Vecd::Zero(); + return acousticCFL_ * smoothing_length_min_ / (reduced_value + TinyReal); +} +//====================================================================================// +StressDiffusion::StressDiffusion(BaseInnerRelation &inner_relation) + : BasePlasticIntegration(inner_relation), + phi_(DynamicCast(this, plastic_continuum_).getFrictionAngle()), + smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), + sound_speed_(plastic_continuum_.ReferenceSoundSpeed()) {} +//====================================================================================// +void StressDiffusion::interaction(size_t index_i, Real dt) +{ + Vecd acc_prior_i = force_prior_[index_i] / mass_[index_i]; + Real gravity = abs(acc_prior_i(1, 0)); + Real density = plastic_continuum_.getDensity(); + Mat3d diffusion_stress_rate = Mat3d::Zero(); + Mat3d diffusion_stress = Mat3d::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; Real r_ij = inner_neighborhood.r_ij_[n]; Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; - Vecd &e_ij = inner_neighborhood.e_ij_[n]; - Real eta_ij = 2 * (0.7 * (Real)Dimensions + 2.1) * (vel_[index_i] - vel_[index_j]).dot(e_ij) / (r_ij + TinyReal); - acceleration += eta_ij * dW_ijV_j * e_ij; + Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0); + diffusion_stress = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; + diffusion_stress(0, 0) -= (1 - sin(phi_)) * density * gravity * y_ij; + diffusion_stress(1, 1) -= density * gravity * y_ij; + diffusion_stress(2, 2) -= (1 - sin(phi_)) * density * gravity * y_ij; + diffusion_stress_rate += 2 * zeta_ * smoothing_length_ * sound_speed_ * + diffusion_stress * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); } - acc_shear_[index_i] += G_ * acceleration * dt / rho_i; + stress_rate_3D_[index_i] = diffusion_stress_rate; } -//=================================================================================================// -ShearStressRelaxation::ShearStressRelaxation(BaseInnerRelation &inner_relation) +//====================================================================================// +ShearStressRelaxationHourglassControl1stHalf :: + ShearStressRelaxationHourglassControl1stHalf(BaseInnerRelation &inner_relation, Real xi) : fluid_dynamics::BaseIntegration(inner_relation), continuum_(DynamicCast(this, particles_->getBaseMaterial())), shear_stress_(particles_->registerStateVariable("ShearStress")), - shear_stress_rate_(particles_->registerStateVariable("ShearStressRate")), velocity_gradient_(particles_->registerStateVariable("VelocityGradient")), strain_tensor_(particles_->registerStateVariable("StrainTensor")), - strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), - von_mises_stress_(particles_->registerStateVariable("VonMisesStress")), - von_mises_strain_(particles_->registerStateVariable("VonMisesStrain")), - Vol_(particles_->getVariableDataByName("VolumetricMeasure")), - B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")) + B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), + scale_penalty_force_(particles_->registerStateVariable("ScalePenaltyForce")), xi_(xi) { particles_->addVariableToSort("ShearStress"); - particles_->addVariableToSort("ShearStressRate"); - particles_->addVariableToSort("VonMisesStress"); - particles_->addVariableToSort("VonMisesStrain"); particles_->addVariableToSort("VelocityGradient"); particles_->addVariableToSort("StrainTensor"); - particles_->addVariableToSort("StrainTensorRate"); -} -//====================================================================================// -void ShearStressRelaxation::initialization(size_t index_i, Real dt) -{ - strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt; - shear_stress_[index_i] += shear_stress_rate_[index_i] * 0.5 * dt; + particles_->addVariableToSort("ScalePenaltyForce"); } //====================================================================================// -void ShearStressRelaxation::interaction(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl1stHalf::interaction(size_t index_i, Real dt) { Matd velocity_gradient = Matd::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + Vecd vel_i = vel_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_i]; + Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; Vecd &e_ij = inner_neighborhood.e_ij_[n]; - Vecd v_ij = vel_[index_i] - vel_[index_j]; - velocity_gradient -= v_ij * (B_[index_i] * e_ij * dW_ijV_j).transpose(); + Matd velocity_gradient_ij = -(vel_i - vel_[index_j]) * (B_[index_i] * e_ij * dW_ijV_j).transpose(); + velocity_gradient += velocity_gradient_ij; } velocity_gradient_[index_i] = velocity_gradient; - /*calculate strain*/ - Matd strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); - strain_tensor_rate_[index_i] = strain_rate; - strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt; - Matd strain_i = strain_tensor_[index_i]; - von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_i); } //====================================================================================// -void ShearStressRelaxation::update(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl1stHalf::update(size_t index_i, Real dt) { - shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); - shear_stress_[index_i] += shear_stress_rate_[index_i] * dt * 0.5; - Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); - von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); + Matd shear_stress_rate = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); + shear_stress_[index_i] += shear_stress_rate * dt; + scale_penalty_force_[index_i] = xi_; + Matd strain_tensor_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); + strain_tensor_[index_i] += strain_tensor_rate * dt; } //====================================================================================// -StressDiffusion::StressDiffusion(BaseInnerRelation &inner_relation) - : BasePlasticIntegration(inner_relation), - fai_(DynamicCast(this, plastic_continuum_).getFrictionAngle()), - smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), - sound_speed_(plastic_continuum_.ReferenceSoundSpeed()) {} +ShearStressRelaxationHourglassControl2ndHalf :: + ShearStressRelaxationHourglassControl2ndHalf(BaseInnerRelation &inner_relation) + : fluid_dynamics::BaseIntegration(inner_relation), + continuum_(DynamicCast(this, particles_->getBaseMaterial())), + shear_stress_(particles_->getVariableDataByName("ShearStress")), + velocity_gradient_(particles_->getVariableDataByName("VelocityGradient")), + acc_shear_(particles_->registerStateVariable("AccelerationByShear")), + acc_hourglass_(particles_->registerStateVariable("AccelerationHourglass")), + scale_penalty_force_(particles_->getVariableDataByName("ScalePenaltyForce")), + G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())) +{ + particles_->addVariableToSort("AccelerationByShear"); + particles_->addVariableToSort("AccelerationHourglass"); +} //====================================================================================// -void StressDiffusion::interaction(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl2ndHalf::interaction(size_t index_i, Real dt) { - Vecd acc_prior_i = force_prior_[index_i] / mass_[index_i]; - Real gravity = abs(acc_prior_i(1, 0)); - Real density = plastic_continuum_.getDensity(); - Mat3d diffusion_stress_rate_ = Mat3d::Zero(); - Mat3d diffusion_stress_ = Mat3d::Zero(); + Real rho_i = rho_[index_i]; + Matd shear_stress_i = shear_stress_[index_i]; + Vecd acceleration = Vecd::Zero(); + Vecd acceleration_hourglass = Vecd::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { size_t index_j = inner_neighborhood.j_[n]; - Real r_ij = inner_neighborhood.r_ij_[n]; Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; - Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0); - diffusion_stress_ = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; - diffusion_stress_(0, 0) -= (1 - sin(fai_)) * density * gravity * y_ij; - diffusion_stress_(1, 1) -= density * gravity * y_ij; - diffusion_stress_(2, 2) -= (1 - sin(fai_)) * density * gravity * y_ij; - diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * - diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); + Vecd &e_ij = inner_neighborhood.e_ij_[n]; + acceleration += ((shear_stress_i + shear_stress_[index_j]) / rho_i) * dW_ijV_j * e_ij; + Vecd v_ij = vel_[index_i] - vel_[index_j]; + Real r_ij = inner_neighborhood.r_ij_[n]; + Vecd v_ij_correction = v_ij - 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_j]) * r_ij * e_ij; + acceleration_hourglass += 0.5 * (scale_penalty_force_[index_i] + scale_penalty_force_[index_j]) * G_ * v_ij_correction * dW_ijV_j * dt / (rho_i * r_ij); } - stress_rate_3D_[index_i] = diffusion_stress_rate_; + acc_hourglass_[index_i] += acceleration_hourglass; + acc_shear_[index_i] = acceleration + acc_hourglass_[index_i]; +} +//====================================================================================// +ShearStressRelaxationHourglassControl1stHalfJ2Plasticity :: + ShearStressRelaxationHourglassControl1stHalfJ2Plasticity(BaseInnerRelation &inner_relation, Real xi) + : ShearStressRelaxationHourglassControl1stHalf(inner_relation, xi), + J2_plasticity_(DynamicCast(this, particles_->getBaseMaterial())), + hardening_factor_(particles_->registerStateVariable("HardeningFactor")) +{ + particles_->addVariableToSort("HardeningFactor"); } //====================================================================================// +void ShearStressRelaxationHourglassControl1stHalfJ2Plasticity::update(size_t index_i, Real dt) +{ + Matd shear_stress_rate = J2_plasticity_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i], hardening_factor_[index_i]); + Matd shear_stress_try = shear_stress_[index_i] + shear_stress_rate * dt; + Real hardening_factor_increment = J2_plasticity_.HardeningFactorRate(shear_stress_try, hardening_factor_[index_i]); + hardening_factor_[index_i] += sqrt(2.0 / 3.0) * hardening_factor_increment; + scale_penalty_force_[index_i] = xi_ * J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]); + shear_stress_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_try, hardening_factor_[index_i]); + Matd strain_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); + strain_tensor_[index_i] += strain_rate * dt; +} } // namespace continuum_dynamics -} // namespace SPH +} // namespace SPH \ No newline at end of file diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 9d7115699d..85362e5bbc 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -49,6 +49,22 @@ class ContinuumInitialCondition : public LocalDynamics Mat3d *stress_tensor_3D_; }; +class AcousticTimeStep : public LocalDynamicsReduce +{ + public: + explicit AcousticTimeStep(SPHBody &sph_body, Real acousticCFL = 0.6); + virtual ~AcousticTimeStep(){}; + Real reduce(size_t index_i, Real dt = 0.0); + virtual Real outputResult(Real reduced_value) override; + + protected: + Fluid &fluid_; + Real *rho_, *p_; + Vecd *vel_; + Real smoothing_length_min_; + Real acousticCFL_; +}; + template class BaseIntegration1stHalf : public FluidDynamicsType { @@ -63,35 +79,6 @@ class BaseIntegration1stHalf : public FluidDynamicsType using Integration1stHalf = BaseIntegration1stHalf; using Integration1stHalfRiemann = BaseIntegration1stHalf; -class ShearAccelerationRelaxation : public fluid_dynamics::BaseIntegration -{ - public: - explicit ShearAccelerationRelaxation(BaseInnerRelation &inner_relation); - virtual ~ShearAccelerationRelaxation(){}; - void interaction(size_t index_i, Real dt = 0.0); - - protected: - GeneralContinuum &continuum_; - Real G_, smoothing_length_; - Vecd *acc_shear_; -}; - -class ShearStressRelaxation : public fluid_dynamics::BaseIntegration -{ - public: - explicit ShearStressRelaxation(BaseInnerRelation &inner_relation); - virtual ~ShearStressRelaxation(){}; - void initialization(size_t index_i, Real dt = 0.0); - void interaction(size_t index_i, Real dt = 0.0); - void update(size_t index_i, Real dt = 0.0); - - protected: - GeneralContinuum &continuum_; - Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_; - Real *von_mises_stress_, *von_mises_strain_, *Vol_; - Matd *B_; -}; - template class BasePlasticIntegration : public fluid_dynamics::BaseIntegration { @@ -103,7 +90,6 @@ class BasePlasticIntegration : public fluid_dynamics::BaseIntegration, RiemannSolverType> protected: RiemannSolverType riemann_solver_; - Real *acc_deviatoric_plastic_strain_, *vertical_stress_; Real *Vol_, *mass_; - Real E_, nu_; - - Real getDeviatoricPlasticStrain(Mat3d &strain_tensor); }; using PlasticIntegration2ndHalfInnerNoRiemann = PlasticIntegration2ndHalf, NoRiemannSolver>; using PlasticIntegration2ndHalfInnerRiemann = PlasticIntegration2ndHalf, AcousticRiemannSolver>; @@ -200,9 +182,49 @@ class StressDiffusion : public BasePlasticIntegration void interaction(size_t index_i, Real dt = 0.0); protected: - Real zeta_ = 0.1, fai_; /*diffusion coefficient*/ + Real zeta_ = 0.1, phi_; /*diffusion coefficient*/ Real smoothing_length_, sound_speed_; }; + +class ShearStressRelaxationHourglassControl1stHalf : public fluid_dynamics::BaseIntegration +{ + public: + explicit ShearStressRelaxationHourglassControl1stHalf(BaseInnerRelation &inner_relation, Real xi = 4.0); + virtual ~ShearStressRelaxationHourglassControl1stHalf(){}; + void interaction(size_t index_i, Real dt = 0.0); + void update(size_t index_i, Real dt = 0.0); + + protected: + GeneralContinuum &continuum_; + Matd *shear_stress_, *velocity_gradient_, *strain_tensor_, *B_; + Real *scale_penalty_force_, xi_; +}; + +class ShearStressRelaxationHourglassControl2ndHalf : public fluid_dynamics::BaseIntegration +{ + public: + explicit ShearStressRelaxationHourglassControl2ndHalf(BaseInnerRelation &inner_relation); + virtual ~ShearStressRelaxationHourglassControl2ndHalf(){}; + void interaction(size_t index_i, Real dt = 0.0); + + protected: + GeneralContinuum &continuum_; + Matd *shear_stress_, *velocity_gradient_; + Vecd *acc_shear_, *acc_hourglass_; + Real *scale_penalty_force_, G_; +}; + +class ShearStressRelaxationHourglassControl1stHalfJ2Plasticity : public ShearStressRelaxationHourglassControl1stHalf +{ + public: + explicit ShearStressRelaxationHourglassControl1stHalfJ2Plasticity(BaseInnerRelation &inner_relation, Real xi = 0.2); + virtual ~ShearStressRelaxationHourglassControl1stHalfJ2Plasticity(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + J2Plasticity &J2_plasticity_; + Real *hardening_factor_; +}; } // namespace continuum_dynamics } // namespace SPH #endif // CONTINUUM_INTEGRATION_H \ No newline at end of file diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp index fa9ca9243b..511c526c74 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp @@ -1,7 +1,5 @@ #pragma once - #include "continuum_integration.h" - #include "base_particles.hpp" namespace SPH { @@ -11,7 +9,7 @@ namespace continuum_dynamics template BaseIntegration1stHalf::BaseIntegration1stHalf(BaseInnerRelation &inner_relation) : FluidDynamicsType(inner_relation), - acc_shear_(this->particles_->template getVariableDataByName("AccelerationByShear")) {} + acc_shear_(this->particles_->template registerStateVariable("AccelerationByShear")) {} //=================================================================================================// template void BaseIntegration1stHalf::update(size_t index_i, Real dt) @@ -28,12 +26,8 @@ BasePlasticIntegration::BasePlasticIntegration(BaseRelationT strain_tensor_3D_(this->particles_->template registerStateVariable("StrainTensor3D")), stress_rate_3D_(this->particles_->template registerStateVariable("StressRate3D")), strain_rate_3D_(this->particles_->template registerStateVariable("StrainRate3D")), - elastic_strain_tensor_3D_(this->particles_->template registerStateVariable("ElasticStrainTensor3D")), - elastic_strain_rate_3D_(this->particles_->template registerStateVariable("ElasticStrainRate3D")), velocity_gradient_(this->particles_->template registerStateVariable("VelocityGradient")) { - this->particles_->template addVariableToSort("ElasticStrainTensor3D"); - this->particles_->template addVariableToSort("ElasticStrainRate3D"); this->particles_->template addVariableToSort("StrainTensor3D"); this->particles_->template addVariableToSort("StressTensor3D"); this->particles_->template addVariableToSort("StrainRate3D"); @@ -142,15 +136,8 @@ template PlasticIntegration2ndHalf, RiemannSolverType>::PlasticIntegration2ndHalf(BaseInnerRelation &inner_relation) : BasePlasticIntegration(inner_relation), riemann_solver_(plastic_continuum_, plastic_continuum_, 20.0 * (Real)Dimensions), - acc_deviatoric_plastic_strain_(particles_->registerStateVariable("AccDeviatoricPlasticStrain")), - vertical_stress_(particles_->registerStateVariable("VerticalStress")), Vol_(particles_->getVariableDataByName("VolumetricMeasure")), - mass_(particles_->getVariableDataByName("Mass")), - E_(plastic_continuum_.getYoungsModulus()), nu_(plastic_continuum_.getPoissonRatio()) -{ - particles_->addVariableToSort("AccDeviatoricPlasticStrain"); - particles_->addVariableToSort("VerticalStress"); -} + mass_(particles_->getVariableDataByName("Mass")) {} //=================================================================================================// template void PlasticIntegration2ndHalf, RiemannSolverType>::initialization(size_t index_i, Real dt) @@ -181,14 +168,6 @@ void PlasticIntegration2ndHalf, RiemannSolverType>::interaction(size_t i } //=================================================================================================// template -Real PlasticIntegration2ndHalf, RiemannSolverType>::getDeviatoricPlasticStrain(Mat3d &strain_tensor) -{ - Mat3d deviatoric_strain_tensor = strain_tensor - (1.0 / (Real)Dimensions) * strain_tensor.trace() * Mat3d::Identity(); - Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum(); - return sqrt(sum * 2.0 / 3.0); -} -//=================================================================================================// -template void PlasticIntegration2ndHalf, RiemannSolverType>::update(size_t index_i, Real dt) { rho_[index_i] += drho_dt_[index_i] * dt * 0.5; @@ -198,18 +177,9 @@ void PlasticIntegration2ndHalf, RiemannSolverType>::update(size_t index_ stress_rate_3D_[index_i] += stress_tensor_rate_3D_; stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt; /*return mapping*/ - Mat3d stress_tensor_ = plastic_continuum_.ReturnMapping(stress_tensor_3D_[index_i]); - stress_tensor_3D_[index_i] = stress_tensor_; - vertical_stress_[index_i] = stress_tensor_3D_[index_i](1, 1); + stress_tensor_3D_[index_i] = plastic_continuum_.ReturnMapping(stress_tensor_3D_[index_i]); strain_rate_3D_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); strain_tensor_3D_[index_i] += strain_rate_3D_[index_i] * dt; - /*calculate elastic strain*/ - Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1.0 / 3.0) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity(); - Real hydrostatic_pressure = (1.0 / 3.0) * stress_tensor_3D_[index_i].trace(); - elastic_strain_tensor_3D_[index_i] = deviatoric_stress / (2.0 * plastic_continuum_.getShearModulus(E_, nu_)) + - hydrostatic_pressure * Mat3d::Identity() / (9.0 * plastic_continuum_.getBulkModulus(E_, nu_)); - Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D_[index_i]; - acc_deviatoric_plastic_strain_[index_i] = getDeviatoricPlasticStrain(plastic_strain_tensor_3D); } //=================================================================================================// template diff --git a/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp b/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp index 650ed6d87c..a99589fc47 100644 --- a/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp +++ b/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp @@ -99,12 +99,10 @@ int main(int ac, char *av[]) Gravity gravity(Vecd(0.0, -gravity_g)); SimpleDynamics> constant_gravity(soil_block, gravity); SimpleDynamics wall_boundary_normal_direction(wall_boundary); - Dynamics1Level granular_stress_relaxation(soil_block_inner, soil_block_contact); Dynamics1Level granular_density_relaxation(soil_block_inner, soil_block_contact); InteractionWithUpdate soil_density_by_summation(soil_block_inner, soil_block_contact); InteractionDynamics stress_diffusion(soil_block_inner); - ReduceDynamics soil_acoustic_time_step(soil_block, 0.4); //---------------------------------------------------------------------- // Define the methods for I/O operations, observations @@ -113,9 +111,10 @@ int main(int ac, char *av[]) BodyStatesRecordingToVtp body_states_recording(sph_system); body_states_recording.addToWrite(soil_block, "Pressure"); body_states_recording.addToWrite(soil_block, "Density"); + SimpleDynamics vertical_stress(soil_block); body_states_recording.addToWrite(soil_block, "VerticalStress"); + SimpleDynamics accumulated_deviatoric_plastic_strain(soil_block); body_states_recording.addToWrite(soil_block, "AccDeviatoricPlasticStrain"); - body_states_recording.addToWrite(wall_boundary, "NormalDirection"); RestartIO restart_io(sph_system); RegressionTestDynamicTimeWarping> write_mechanical_energy(soil_block, gravity); @@ -135,8 +134,8 @@ int main(int ac, char *av[]) int screen_output_interval = 500; int observation_sample_interval = screen_output_interval * 2; int restart_output_interval = screen_output_interval * 10; - Real End_Time = 1.0; /**< End time. */ - Real D_Time = End_Time / 50; /**< Time stamps for output of body states. */ + Real End_Time = 0.8; /**< End time. */ + Real D_Time = End_Time / 40; /**< Time stamps for output of body states. */ Real Dt = 0.1 * D_Time; //---------------------------------------------------------------------- // Statistics for CPU time @@ -205,6 +204,8 @@ int main(int ac, char *av[]) interval_updating_configuration += TickCount::now() - time_instance; } TickCount t2 = TickCount::now(); + vertical_stress.exec(); + accumulated_deviatoric_plastic_strain.exec(); body_states_recording.writeToFile(); TickCount t3 = TickCount::now(); interval += t3 - t2; diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml index df2b1d25ba..461e092951 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml @@ -1,51 +1,42 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml index 296f1f7722..b376707536 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml deleted file mode 100644 index 49c242b5e7..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml deleted file mode 100644 index fac5f2d438..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml deleted file mode 100644 index 8a1dc7daa1..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml index 9cacfb4cc1..f63c0bb733 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml index 20390edd95..a85d9d3ba4 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml index 9e018fae74..0e388bbb73 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml index fe62912d6f..e0347772fe 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml index ffb3d2a9ac..238c385f80 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml index fd67bd9b30..9e22ce14bf 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml index f71552ddb8..4facd994cc 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml deleted file mode 100644 index 8297225a57..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml deleted file mode 100644 index f02d5b8586..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml index f3f38a290a..9f850afadb 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat index 19c6fe5aad..ffe6691bfd 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat @@ -1,3 +1,3 @@ true -13 +8 4 \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp index f28c809311..aca9b8b2e9 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -1,8 +1,9 @@ /* ----------------------------------------------------------------------------* - * SPHinXsys: 2D oscillation beam--updated Lagrangian method * + * SPHinXsys: 2D oscillation beam * * ----------------------------------------------------------------------------* * This is the one of the basic test cases for understanding SPH method for * - * solid simulation based on updated Lagrangian method * + * solid simulation based on updated Lagrangian method. * + * A generalized hourglass control method is used here. * * In this case, the constraint of the beam is implemented with * * internal constrained subregion. * * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu * @@ -26,9 +27,7 @@ BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus Real poisson = 0.3975; // Poisson ratio -// Real poisson = 0.4; //Poisson ratio Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -Real gravity_g = 0.0; //---------------------------------------------------------------------- // Parameters for initial condition on velocity //---------------------------------------------------------------------- @@ -109,7 +108,6 @@ int main(int ac, char *av[]) // Creating body, materials and particles. //---------------------------------------------------------------------- RealBody beam_body(sph_system, makeShared("BeamBody")); - beam_body.sph_adaptation_->resetKernel(); beam_body.defineMaterial(rho0_s, c0, Youngs_modulus, poisson); beam_body.generateParticles(); @@ -130,17 +128,15 @@ int main(int ac, char *av[]) // this section define all numerical methods will be used in this case //----------------------------------------------------------------------------- /** initial condition */ - InteractionWithUpdate correction_matrix(beam_body_inner); SimpleDynamics beam_initial_velocity(beam_body); - - InteractionDynamics beam_shear_acceleration(beam_body_inner); - Dynamics1Level beam_shear_stress_relaxation(beam_body_inner); + InteractionWithUpdate correction_matrix(beam_body_inner); Dynamics1Level beam_pressure_relaxation(beam_body_inner); Dynamics1Level beam_density_relaxation(beam_body_inner); + InteractionWithUpdate beam_shear_stress(beam_body_inner); + InteractionDynamics beam_shear_acceleration(beam_body_inner); SimpleDynamics beam_volume_update(beam_body); - - ReduceDynamics fluid_advection_time_step(beam_body, U_ref, 0.2); - ReduceDynamics fluid_acoustic_time_step(beam_body, 0.4); + ReduceDynamics advection_time_step(beam_body, U_ref, 0.2); + ReduceDynamics acoustic_time_step(beam_body, 0.4); // clamping a solid body part. BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); SimpleDynamics constraint_beam_base(beam_base); @@ -148,10 +144,10 @@ int main(int ac, char *av[]) // outputs //----------------------------------------------------------------------------- BodyStatesRecordingToVtp write_beam_states(beam_body); - write_beam_states.addToWrite(beam_body, "VonMisesStress"); - write_beam_states.addToWrite(beam_body, "VonMisesStrain"); write_beam_states.addToWrite(beam_body, "Density"); write_beam_states.addToWrite(beam_body, "Pressure"); + SimpleDynamics beam_von_mises_stress(beam_body); + write_beam_states.addToWrite(beam_body, "VonMisesStress"); ObservedQuantityRecording write_beam_tip_displacement("Position", beam_observer_contact); RegressionTestDynamicTimeWarping> write_beam_kinetic_energy(beam_body); //---------------------------------------------------------------------- @@ -184,17 +180,16 @@ int main(int ac, char *av[]) while (integration_time < D_Time) { Real relaxation_time = 0.0; - Real advection_dt = fluid_advection_time_step.exec(); + Real advection_dt = advection_time_step.exec(); beam_volume_update.exec(); - while (relaxation_time < advection_dt) { - Real acoustic_dt = fluid_acoustic_time_step.exec(); - beam_shear_stress_relaxation.exec(acoustic_dt); + Real acoustic_dt = acoustic_time_step.exec(); beam_pressure_relaxation.exec(acoustic_dt); constraint_beam_base.exec(); - beam_density_relaxation.exec(acoustic_dt); + beam_shear_stress.exec(acoustic_dt); beam_shear_acceleration.exec(acoustic_dt); + beam_density_relaxation.exec(acoustic_dt); ite++; relaxation_time += acoustic_dt; integration_time += acoustic_dt; @@ -211,6 +206,7 @@ int main(int ac, char *av[]) beam_body_inner.updateConfiguration(); correction_matrix.exec(); } + beam_von_mises_stress.exec(); write_beam_tip_displacement.writeToFile(ite); write_beam_kinetic_energy.writeToFile(ite); TickCount t2 = TickCount::now(); @@ -231,6 +227,5 @@ int main(int ac, char *av[]) { write_beam_kinetic_energy.testResult(); } - return 0; } \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy.xml new file mode 100644 index 0000000000..c6b0a3834f --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy.xml @@ -0,0 +1,103 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml index 1f10003773..7a8cacca58 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml new file mode 100644 index 0000000000..4d9716f718 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml new file mode 100644 index 0000000000..8de5e83ed5 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml index bb2d859705..4f75fb4c37 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml new file mode 100644 index 0000000000..9ada01ea11 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml index 5789bcd8a4..0313a9973c 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml index 8be7b3f132..e59f0fa6a5 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt b/tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt deleted file mode 100644 index 23ebe87350..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) -PROJECT("${CURRENT_FOLDER}") - -SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) -SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") -SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") -SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") - -file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) -file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ - DESTINATION ${BUILD_INPUT_PATH}) - -add_executable(${PROJECT_NAME}) -aux_source_directory(. DIR_SRCS) -target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) -target_link_libraries(${PROJECT_NAME} sphinxsys_3d) -set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") - -add_test(NAME ${PROJECT_NAME} - COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp deleted file mode 100644 index dceabddb3d..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp +++ /dev/null @@ -1,261 +0,0 @@ -/** - * @file oscillating_plate_UL.cpp - * @brief This is the test case for the hourglass manuscript. - * @details We consider vibration deformation of a square plate under initial vertical velocity field. - * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu - */ -#include "sphinxsys.h" -using namespace SPH; -//---------------------------------------------------------------------- -// Basic geometry parameters and numerical setup. -//---------------------------------------------------------------------- -Real PL = 0.4; /** Length of the square plate. */ -Real PH = 0.4; /** Width of the square plate. */ -Real PT = 0.01; /** Thickness of the square plate. */ -int particle_number = 3; /** Particle number in the direction of the thickness. */ -Real particle_spacing_ref = PT / (Real)particle_number; /** Initial reference particle spacing. */ -int particle_number_PL = PL / particle_spacing_ref; -int particle_number_PH = PH / particle_spacing_ref; -int BWD = 1; -Real BW = particle_spacing_ref * (Real)BWD; /** Boundary width, determined by specific layer of boundary particles. */ -BoundingBox system_domain_bounds(Vec3d(-BW, -BW, -PT / 2), Vec3d(PL + BW, PH + BW, PT / 2)); -StdVec observation_location = {Vecd(0.5 * PL, 0.5 * PH, 0.0), Vecd(-BW, -BW, 0.0)}; -//---------------------------------------------------------------------- -// Material parameters. -//---------------------------------------------------------------------- -Real rho0_s = 1000.0; /** Normalized density. */ -Real Youngs_modulus = 100.0e6; /** Normalized Youngs Modulus. */ -Real poisson = 0.3; /** Poisson ratio. */ -Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -Real gravity_g = 0.0; - -Real governing_vibration_integer_x = 2.0; -Real governing_vibration_integer_y = 2.0; -Real U_ref = 1.0; // Maximum velocity - -namespace SPH -{ -//---------------------------------------------------------------------- -// Complex shape for wall boundary, note that no partial overlap is allowed -// for the shapes in a complex shape. -//---------------------------------------------------------------------- -class Plate; -template <> -class ParticleGenerator : public ParticleGenerator -{ - public: - explicit ParticleGenerator(SPHBody &sph_body, BaseParticles &base_particles) - : ParticleGenerator(sph_body, base_particles){}; - virtual void prepareGeometricData() override - { - for (int k = 0; k < particle_number; k++) - { - for (int i = 0; i < particle_number_PL + 2 * BWD; i++) - { - for (int j = 0; j < particle_number_PH + 2 * BWD; j++) - { - Real x = particle_spacing_ref * i - BW + particle_spacing_ref * 0.5; - Real y = particle_spacing_ref * j - BW + particle_spacing_ref * 0.5; - Real z = particle_spacing_ref * (k - ((particle_number - 1.0) / 2.0)); - addPositionAndVolumetricMeasure( - Vecd(x, y, z), particle_spacing_ref * particle_spacing_ref * particle_spacing_ref); - } - } - } - } -}; - -/** Define the boundary geometry. */ -class BoundaryGeometry : public BodyPartByParticle -{ - public: - BoundaryGeometry(SPHBody &body, const std::string &body_part_name) - : BodyPartByParticle(body, body_part_name) - { - TaggingParticleMethod tagging_particle_method = std::bind(&BoundaryGeometry::tagManually, this, _1); - tagParticles(tagging_particle_method); - }; - virtual ~BoundaryGeometry(){}; - - private: - void tagManually(size_t index_i) - { - if ((base_particles_.ParticlePositions()[index_i][2] < 0.5 * particle_spacing_ref) && - (base_particles_.ParticlePositions()[index_i][2] > -0.5 * particle_spacing_ref) && - (base_particles_.ParticlePositions()[index_i][0] < 0.0 || - base_particles_.ParticlePositions()[index_i][0] > PL || - base_particles_.ParticlePositions()[index_i][1] < 0.0 || - base_particles_.ParticlePositions()[index_i][1] > PH)) - { - body_part_particles_.push_back(index_i); - } - }; -}; - -/** Define the initial condition. */ -class BeamInitialCondition - : public fluid_dynamics::FluidInitialCondition -{ - public: - explicit BeamInitialCondition(SPHBody &sph_body) - : fluid_dynamics::FluidInitialCondition(sph_body){}; - - void update(size_t index_i, Real dt) - { - /** initial velocity profile */ - vel_[index_i][2] = sin(governing_vibration_integer_x * Pi * pos_[index_i][0] / PL) * - sin(governing_vibration_integer_y * Pi * pos_[index_i][1] / PH); - }; -}; -} // namespace SPH - -//---------------------------------------------------------------------- -// Main program starts here. -//---------------------------------------------------------------------- -int main(int ac, char *av[]) -{ - //---------------------------------------------------------------------- - // Build up an SPHSystem and IO environment. - //---------------------------------------------------------------------- - SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); - sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment(); - //---------------------------------------------------------------------- - // Creating bodies with corresponding materials and particles. - //---------------------------------------------------------------------- - SolidBody plate_body(sph_system, makeShared("PlateBody")); - plate_body.defineMaterial(rho0_s, c0, Youngs_modulus, poisson); - plate_body.generateParticles(); - - ObserverBody plate_observer(sph_system, "PlateObserver"); - plate_observer.generateParticles(observation_location); - //---------------------------------------------------------------------- - // Define body relation map. - // The contact map gives the topological connections between the bodies. - // Basically the the range of bodies to build neighbor particle lists. - // Generally, we first define all the inner relations, then the contact relations. - //---------------------------------------------------------------------- - InnerRelation plate_body_inner(plate_body); - ContactRelation plate_observer_contact(plate_observer, {&plate_body}); - //---------------------------------------------------------------------- - // Define the numerical methods used in the simulation. - // Note that there may be data dependence on the sequence of constructions. - // Generally, the geometric models or simple objects without data dependencies, - // such as gravity, should be initiated first. - // Then the major physical particle dynamics model should be introduced. - // Finally, the auxillary models such as time step estimator, initial condition, - // boundary condition and other constraints should be defined. - //---------------------------------------------------------------------- - SimpleDynamics initial_velocity(plate_body); - InteractionWithUpdate corrected_configuration(plate_body_inner); - - InteractionDynamics plate_shear_acceleration(plate_body_inner); - Dynamics1Level plate_pressure_relaxation(plate_body_inner); - Dynamics1Level plate_density_relaxation(plate_body_inner); - Dynamics1Level plate_shear_stress_relaxation(plate_body_inner); - - ReduceDynamics fluid_advection_time_step(plate_body, U_ref, 0.2); - ReduceDynamics fluid_acoustic_time_step(plate_body, 0.4); - BoundaryGeometry boundary_geometry(plate_body, "BoundaryGeometry"); - SimpleDynamics constrain_holder(boundary_geometry, Vecd(1.0, 1.0, 0.0)); - SimpleDynamics constrain_mass_center(plate_body, Vecd(1.0, 1.0, 0.0)); - //---------------------------------------------------------------------- - // Define the methods for I/O operations, observations - // and regression tests of the simulation. - //---------------------------------------------------------------------- - BodyStatesRecordingToVtp write_states(sph_system); - write_states.addToWrite(plate_body, "VonMisesStress"); - write_states.addToWrite(plate_body, "VonMisesStress"); - write_states.addToWrite(plate_body, "Pressure"); - write_states.addToWrite(plate_body, "Density"); - RestartIO restart_io(sph_system); - ObservedQuantityRecording write_plate_displacement("Position", plate_observer_contact); - RegressionTestDynamicTimeWarping> write_kinetic_energy(plate_body); - //---------------------------------------------------------------------- - // Prepare the simulation with cell linked list, configuration - // and case specified initial condition if necessary. - //---------------------------------------------------------------------- - sph_system.initializeSystemCellLinkedLists(); - sph_system.initializeSystemConfigurations(); - initial_velocity.exec(); - constrain_holder.exec(); - corrected_configuration.exec(); - //---------------------------------------------------------------------- - // Setup for time-stepping control - //---------------------------------------------------------------------- - Real &physical_time = *sph_system.getSystemVariableDataByName("PhysicalTime"); - size_t number_of_iterations = 0; - int screen_output_interval = 500; - int restart_output_interval = screen_output_interval * 10; - Real end_time = 0.02; - Real output_period = end_time / 50.0; - //---------------------------------------------------------------------- - // Statistics for CPU time - //---------------------------------------------------------------------- - TickCount t1 = TickCount::now(); - TimeInterval interval; - //---------------------------------------------------------------------- - // First output before the main loop. - //---------------------------------------------------------------------- - write_states.writeToFile(0); - write_plate_displacement.writeToFile(0); - write_kinetic_energy.writeToFile(0); - //---------------------------------------------------------------------- - // Main loop starts here. - //---------------------------------------------------------------------- - while (physical_time < end_time) - { - Real integration_time = 0.0; - while (integration_time < output_period) - { - Real relaxation_time = 0.0; - Real advection_dt = fluid_advection_time_step.exec(); - - while (relaxation_time < advection_dt) - { - Real acoustic_dt = fluid_acoustic_time_step.exec(); - plate_shear_stress_relaxation.exec(acoustic_dt); - plate_pressure_relaxation.exec(acoustic_dt); - constrain_holder.exec(acoustic_dt); - plate_density_relaxation.exec(acoustic_dt); - plate_shear_acceleration.exec(acoustic_dt); - number_of_iterations++; - relaxation_time += acoustic_dt; - integration_time += acoustic_dt; - physical_time += acoustic_dt; - if (number_of_iterations % screen_output_interval == 0) - { - std::cout << "N=" << number_of_iterations << " Time: " - << physical_time << " advection_dt: " - << advection_dt << " acoustic_dt: " - << acoustic_dt << "\n"; - if (number_of_iterations % restart_output_interval == 0 && number_of_iterations != sph_system.RestartStep()) - restart_io.writeToFile(Real(number_of_iterations)); - } - } - plate_body.updateCellLinkedList(); - plate_body_inner.updateConfiguration(); - corrected_configuration.exec(); - } - write_plate_displacement.writeToFile(number_of_iterations); - write_kinetic_energy.writeToFile(number_of_iterations); - TickCount t2 = TickCount::now(); - write_states.writeToFile(); - TickCount t3 = TickCount::now(); - interval += t3 - t2; - } - - TickCount t4 = TickCount::now(); - TickCount::interval_t tt; - tt = t4 - t1 - interval; - std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; - - if (sph_system.GenerateRegressionData()) - { - write_kinetic_energy.generateDataBase(1.0e-3); - } - else - { - write_kinetic_energy.testResult(); - } - return 0; -} diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml deleted file mode 100644 index 0595b41a58..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml deleted file mode 100644 index eedda0a708..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml deleted file mode 100644 index 89c7502971..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp b/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp index c80282c26c..dae594254e 100644 --- a/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp +++ b/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp @@ -171,7 +171,9 @@ int main(int ac, char *av[]) BodyStatesRecordingToVtp body_states_recording(sph_system); body_states_recording.addToWrite(soil_block, "Density"); body_states_recording.addToWrite(soil_block, "Pressure"); + SimpleDynamics vertical_stress(soil_block); body_states_recording.addToWrite(soil_block, "VerticalStress"); + SimpleDynamics accumulated_deviatoric_plastic_strain(soil_block); body_states_recording.addToWrite(soil_block, "AccDeviatoricPlasticStrain"); RestartIO restart_io(sph_system); RegressionTestDynamicTimeWarping> write_soil_mechanical_energy(soil_block, gravity); @@ -263,6 +265,8 @@ int main(int ac, char *av[]) interval_updating_configuration += TickCount::now() - time_instance; } TickCount t2 = TickCount::now(); + vertical_stress.exec(); + accumulated_deviatoric_plastic_strain.exec(); body_states_recording.writeToFile(); TickCount t3 = TickCount::now(); interval += t3 - t2; diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt new file mode 100644 index 0000000000..acb72b326d --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt @@ -0,0 +1,26 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_3d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + +set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") + +set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "continuum_dynamics") diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml new file mode 100644 index 0000000000..da81bfd22e --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml new file mode 100644 index 0000000000..60bdb4418f --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml new file mode 100644 index 0000000000..d23ad8ee5c --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml new file mode 100644 index 0000000000..4764515cdb --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml new file mode 100644 index 0000000000..0d57d75632 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml new file mode 100644 index 0000000000..3886ec8fdb --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml new file mode 100644 index 0000000000..2e8750dac5 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml similarity index 55% rename from tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml rename to tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml index 9497963594..bcf6137195 100644 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_runtimes.dat b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat similarity index 100% rename from tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_runtimes.dat rename to tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py similarity index 80% rename from tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py rename to tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py index 52acb2cabc..2612da0e6a 100644 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py @@ -7,12 +7,12 @@ from regression_test_base_tool import SphinxsysRegressionTest """ -case name: test_3d_oscillating_plate_UL +case name: test_3d_taylor_bar """ -case_name = "test_3d_oscillating_plate_UL" -body_name = "PlateBody" -parameter_name = "TotalKineticEnergy" +case_name = "test_3d_taylor_bar_UL" +body_name = "MyObserver" +parameter_name = "Position" number_of_run_times = 0 converged = 0 @@ -21,12 +21,13 @@ while True: print("Now start a new run......") - sphinxsys.run_case() + sphinxsys.run_particle_relaxation() + sphinxsys.run_case_with_reload() number_of_run_times += 1 converged = sphinxsys.read_dat_file() print("Please note: This is the", number_of_run_times, "run!") if number_of_run_times <= 200: - if converged == "true": + if (converged == "true"): print("The tested parameters of all variables are converged, and the run will stop here!") break elif converged != "true": diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp new file mode 100644 index 0000000000..943e53bc6f --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -0,0 +1,179 @@ +/** + * @file taylor_bar.cpp + * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu + */ +#include "taylor_bar_UL.h" +int main(int ac, char *av[]) +{ + /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ + SPHSystem system(system_domain_bounds, particle_spacing_ref); + system.setRunParticleRelaxation(false); + system.setReloadParticles(true); + system.handleCommandlineOptions(ac, av)->setIOEnvironment(); + + RealBody column(system, makeShared("Column")); + column.defineBodyLevelSetShape()->writeLevelSet(system); + column.defineMaterial(rho0_s, c0, Youngs_modulus, poisson, yield_stress); + (!system.RunParticleRelaxation() && system.ReloadParticles()) + ? column.generateParticles(column.getName()) + : column.generateParticles(); + + SolidBody wall_boundary(system, makeShared("Wall")); + wall_boundary.defineMaterial(rho0_s, Youngs_modulus, poisson); + wall_boundary.generateParticles(); + + ObserverBody my_observer(system, "MyObserver"); + StdVec observation_location = {Vecd(0.0, 0.0, PW)}; + my_observer.generateParticles(observation_location); + + /**body relation topology */ + InnerRelation column_inner(column); + ContactRelation my_observer_contact(my_observer, {&column}); + SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); + //---------------------------------------------------------------------- + // Run particle relaxation for body-fitted distribution if chosen. + //---------------------------------------------------------------------- + if (system.RunParticleRelaxation()) + { + using namespace relax_dynamics; + /** Random reset the insert body particle position. */ + SimpleDynamics random_column_particles(column); + /** Write the body state to Vtp file. */ + BodyStatesRecordingToVtp write_column_to_vtp(system); + /** Write the particle reload files. */ + ReloadParticleIO write_particle_reload_files(column); + /** A Physics relaxation step. */ + RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); + /** + * @brief Particle relaxation starts here. + */ + random_column_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + write_column_to_vtp.writeToFile(0); + /** relax particles of the insert body. */ + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 200 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the column body N = " << ite_p << "\n"; + write_column_to_vtp.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process of cylinder body finish !" << std::endl; + /** Output results. */ + write_particle_reload_files.writeToFile(0.0); + return 0; + } + //---------------------------------------------------------------------- + // All numerical methods will be used in this case. + //---------------------------------------------------------------------- + SimpleDynamics initial_condition(column); + SimpleDynamics wall_normal_direction(wall_boundary); + InteractionWithUpdate corrected_configuration(column_inner); + Dynamics1Level column_pressure_relaxation(column_inner); + Dynamics1Level column_density_relaxation(column_inner); + InteractionWithUpdate column_shear_stress(column_inner); + InteractionDynamics column_shear_acceleration(column_inner); + SimpleDynamics column_volume_update(column); + ReduceDynamics advection_time_step(column, U_max, 0.2); + ReduceDynamics acoustic_time_step(column, 0.4); + InteractionDynamics column_wall_contact_force(column_wall_contact); + //---------------------------------------------------------------------- + // Output + //---------------------------------------------------------------------- + /**define simple data file input and outputs functions. */ + BodyStatesRecordingToVtp write_states(system); + write_states.addToWrite(column, "Pressure"); + write_states.addToWrite(column, "Density"); + SimpleDynamics column_von_mises_stress(column); + write_states.addToWrite(column, "VonMisesStress"); + RegressionTestDynamicTimeWarping> write_displacement("Position", my_observer_contact); + //---------------------------------------------------------------------- + // From here the time stepping begins. + //---------------------------------------------------------------------- + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + wall_normal_direction.exec(); + corrected_configuration.exec(); + initial_condition.exec(); + //---------------------------------------------------------------------- + // Setup time-stepping related simulation parameters. + //---------------------------------------------------------------------- + Real &physical_time = *system.getSystemVariableDataByName("PhysicalTime"); + int ite = 0; + Real end_time = 6.0e-5; + int screen_output_interval = 100; + Real output_period = end_time / 60; + /** Statistics for computing time. */ + TickCount t1 = TickCount::now(); + TimeInterval interval; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + write_states.writeToFile(); + write_displacement.writeToFile(0); + //---------------------------------------------------------------------- + // Main time-stepping loop. + //---------------------------------------------------------------------- + while (physical_time < end_time) + { + Real integration_time = 0.0; + while (integration_time < output_period) + { + Real relaxation_time = 0.0; + Real advection_dt = advection_time_step.exec(); + column_volume_update.exec(); + while (relaxation_time < advection_dt) + { + Real acoustic_dt = acoustic_time_step.exec(); + if (ite % screen_output_interval == 0) + { + std::cout << "N=" << ite << " Time: " + << physical_time << " advection_dt: " + << advection_dt << " acoustic_dt: " + << acoustic_dt << "\n"; + } + column_wall_contact_force.exec(acoustic_dt); + + column_pressure_relaxation.exec(acoustic_dt); + column_shear_stress.exec(acoustic_dt); + column_shear_acceleration.exec(acoustic_dt); + column_density_relaxation.exec(acoustic_dt); + + ite++; + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + physical_time += acoustic_dt; + } + column.updateCellLinkedList(); + column_inner.updateConfiguration(); + column_wall_contact.updateConfiguration(); + corrected_configuration.exec(); + } + TickCount t2 = TickCount::now(); + column_von_mises_stress.exec(); + write_states.writeToFile(ite); + write_displacement.writeToFile(ite); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall_boundary time for computation: " << tt.seconds() << " seconds." << std::endl; + + if (system.GenerateRegressionData()) + { + write_displacement.generateDataBase(0.1); + } + else + { + write_displacement.testResult(); + } + + return 0; +} diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h new file mode 100644 index 0000000000..2e3e3df898 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h @@ -0,0 +1,132 @@ +/** + * @file taylor_bar.cpp + * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu + */ +#pragma once +#include "sphinxsys.h" +using namespace SPH; +//---------------------------------------------------------------------- +// Global geometry parameters. +//---------------------------------------------------------------------- +Real PL = 0.00391; /**< X-direction domain. */ +Real PW = 0.02346; /**< Z-direction domain. */ +Real particle_spacing_ref = PL / 12.0; +Real SL = particle_spacing_ref * 4.0; +Real inner_circle_radius = PL; +Vec3d domain_lower_bound(-4.0 * PL, -4.0 * PL, -SL); +Vec3d domain_upper_bound(4.0 * PL, 4.0 * PL, 2.0 * PW); +BoundingBox system_domain_bounds(domain_lower_bound, domain_upper_bound); +int resolution(20); +//---------------------------------------------------------------------- +// Material properties and global parameters +//---------------------------------------------------------------------- +Real rho0_s = 2700.0; /**< Reference density. */ +Real poisson = 0.3; /**< Poisson ratio. */ +Real Youngs_modulus = 78.2e9; +Real yield_stress = 0.29e9; +Real vel_0 = 373.0; +Real U_max = vel_0; +Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); +/** Define the wall. */ +class WallBoundary : public ComplexShape +{ + public: + explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) + { + Vecd halfsize_holder(3.0 * PL, 3.0 * PL, 0.5 * SL); + Vecd translation_holder(0.0, 0.0, -0.5 * SL); + add(halfsize_holder, resolution, translation_holder); + } +}; +/** Define the body. */ +class Column : public ComplexShape +{ +public: + explicit Column(const std::string& shape_name) : ComplexShape(shape_name) + { + Vecd translation_column(0.0, 0.0, 0.5 * PW + particle_spacing_ref); + add(SimTK::UnitVec3(0, 0, 1.0), inner_circle_radius, + 0.5 * PW, resolution, translation_column); + } +}; +/** + * application dependent initial condition + */ +class InitialCondition + : public fluid_dynamics::FluidInitialCondition +{ + public: + explicit InitialCondition(SPHBody &sph_body) + : fluid_dynamics::FluidInitialCondition(sph_body){}; + + void update(size_t index_i, Real dt) + { + vel_[index_i][2] = -vel_0; + } +}; +/** Contact force. */ +class DynamicContactForceWithWall : public LocalDynamics, + public DataDelegateContact +{ + public: + explicit DynamicContactForceWithWall(SurfaceContactRelation &solid_body_contact_relation, Real penalty_strength = 1.0) + : LocalDynamics(solid_body_contact_relation.getSPHBody()), + DataDelegateContact(solid_body_contact_relation), + continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), + Vol_(particles_->getVariableDataByName("VolumetricMeasure")), + vel_(particles_->getVariableDataByName("Velocity")), + force_prior_(particles_->getVariableDataByName("ForcePrior")), + penalty_strength_(penalty_strength) + { + impedance_ = continuum_.ReferenceDensity() * sqrt(continuum_.ContactStiffness()); + reference_pressure_ = continuum_.ReferenceDensity() * continuum_.ContactStiffness(); + for (size_t k = 0; k != contact_particles_.size(); ++k) + { + contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName("VolumetricMeasure")); + contact_vel_.push_back(contact_particles_[k]->registerStateVariable("Velocity")); + contact_n_.push_back(contact_particles_[k]->template getVariableDataByName("NormalDirection")); + } + }; + virtual ~DynamicContactForceWithWall(){}; + void interaction(size_t index_i, Real dt = 0.0) + { + Vecd force = Vecd::Zero(); + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Real particle_spacing_j1 = 1.0 / contact_bodies_[k]->sph_adaptation_->ReferenceSpacing(); + Real particle_spacing_ratio2 = + 1.0 / (sph_body_.sph_adaptation_->ReferenceSpacing() * particle_spacing_j1); + particle_spacing_ratio2 *= 0.1 * particle_spacing_ratio2; + + Vecd *n_k = contact_n_[k]; + Vecd *vel_n_k = contact_vel_[k]; + Real *Vol_k = contact_Vol_[k]; + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + Vecd e_ij = contact_neighborhood.e_ij_[n]; + Vecd n_k_j = n_k[index_j]; + Real impedance_p = 0.5 * impedance_ * (vel_[index_i] - vel_n_k[index_j]).dot(-n_k_j); + Real overlap = contact_neighborhood.r_ij_[n] * n_k_j.dot(e_ij); + Real delta = 2.0 * overlap * particle_spacing_j1; + Real beta = delta < 1.0 ? (1.0 - delta) * (1.0 - delta) * particle_spacing_ratio2 : 0.0; + Real penalty_p = penalty_strength_ * beta * fabs(overlap) * reference_pressure_; + // force due to pressure + force -= 2.0 * (impedance_p + penalty_p) * e_ij.dot(n_k_j) * + n_k_j * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; + } + } + force_prior_[index_i] += force * Vol_[index_i]; + }; + + protected: + GeneralContinuum &continuum_; + Real *Vol_; + Vecd *vel_, *force_prior_; // note that prior force directly used here + StdVec contact_Vol_; + StdVec contact_vel_, contact_n_; + Real penalty_strength_; + Real impedance_, reference_pressure_; +};