Skip to content

Commit

Permalink
Merge pull request #651 from Shuaihao-Zhang/SPH-GNOG
Browse files Browse the repository at this point in the history
Generalized non-hourglass formulation for ULSPH solid dynamics
  • Loading branch information
Xiangyu-Hu authored Sep 19, 2024
2 parents 8b33bac + c302adc commit 599c900
Show file tree
Hide file tree
Showing 53 changed files with 1,067 additions and 598 deletions.
49 changes: 44 additions & 5 deletions src/shared/materials/general_continuum.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)));
}
Expand All @@ -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();
Expand All @@ -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
30 changes: 28 additions & 2 deletions src/shared/materials/general_continuum.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/

Expand All @@ -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;
Expand Down Expand Up @@ -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
Original file line number Diff line number Diff line change
Expand Up @@ -30,3 +30,4 @@

#include "base_continuum_dynamics.h"
#include "continuum_integration.hpp"
#include "continuum_dynamics_variable.h"
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
#include "continuum_dynamics_variable.h"

namespace SPH
{
namespace continuum_dynamics
{
//=============================================================================================//
VonMisesStress::VonMisesStress(SPHBody &sph_body)
: BaseDerivedVariable<Real>(sph_body, "VonMisesStress"),
p_(particles_->getVariableDataByName<Real>("Pressure")),
shear_stress_(particles_->getVariableDataByName<Matd>("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<Real>(sph_body, "VonMisesStrain"),
strain_tensor_(particles_->getVariableDataByName<Matd>("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<Real>(sph_body, "VerticalStress"),
stress_tensor_3D_(particles_->getVariableDataByName<Mat3d>("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<Real>(sph_body, "AccDeviatoricPlasticStrain"),
plastic_continuum_(DynamicCast<PlasticContinuum>(this, sph_body_.getBaseMaterial())),
stress_tensor_3D_(particles_->getVariableDataByName<Mat3d>("StressTensor3D")),
strain_tensor_3D_(particles_->getVariableDataByName<Mat3d>("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
Original file line number Diff line number Diff line change
@@ -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<Real>
{
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<Real>
{
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<Real>
{
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<Real>
{
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
Loading

0 comments on commit 599c900

Please sign in to comment.