Skip to content

Fix orthotropic material #684

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
WeiyiVirtonomy wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Fix orthotropic material #684

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions src/for_3D_build/materials/elastic_solid_3d.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,18 +7,18 @@ namespace SPH
void OrthotropicSolid::CalculateAllMu()
{

Mu_[0] = 1 / G_[0] + 1 / G_[2] - 1 / G_[1];
Mu_[1] = 1 / G_[1] + 1 / G_[0] - 1 / G_[2];
Mu_[2] = 1 / G_[2] + 1 / G_[1] - 1 / G_[0];
Mu_[0] = G_[0] + G_[2] - G_[1];
Mu_[1] = G_[1] + G_[0] - G_[2];
Mu_[2] = G_[2] + G_[1] - G_[0];
}
//=================================================================================================//
void OrthotropicSolid::CalculateAllLambda()
{
// first we calculate the upper left part, a 3x3 matrix of the full compliance matrix
Matd Compliance = Matd::Zero();
Compliance.col(0) = Vecd(1 / E_[0], -poisson_[0] / E_[1], -poisson_[1] / E_[2]);
Compliance.col(0) = Vecd(1 / E_[0], -poisson_[0] / E_[0], -poisson_[1] / E_[2]);
Compliance.col(1) = Vecd(-poisson_[0] / E_[0], 1 / E_[1], -poisson_[2] / E_[1]);
Compliance.col(2) = Vecd(-poisson_[1] / E_[0], -poisson_[2] / E_[1], 1 / E_[2]);
Compliance.col(2) = Vecd(-poisson_[1] / E_[2], -poisson_[2] / E_[1], 1 / E_[2]);
// we calculate the inverse of the Compliance matrix, and calculate the lambdas elementwise
Matd Compliance_inv = Compliance.inverse();
// Lambda_ is a 3x3 matrix
Expand Down
28 changes: 24 additions & 4 deletions src/shared/materials/elastic_solid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,9 +147,8 @@ Real NeoHookeanSolidIncompressible::VolumetricKirchhoff(Real J)
//=================================================================================================//
OrthotropicSolid::OrthotropicSolid(Real rho_0, std::array<Vecd, Dimensions> a, std::array<Real, Dimensions> E,
std::array<Real, Dimensions> G, std::array<Real, Dimensions> poisson)
// set parameters for parent class: LinearElasticSolid
// we take the max. E and max. poisson to approximate the maximum of
// the Bulk modulus --> for time step size calculation
// since Poisson ratio can be larger than 0.5 for orthotropic materials,
// the sound speed needs to be reset after parameters are computed
: LinearElasticSolid(rho_0, *std::max_element(E.cbegin(), E.cend()),
*std::max_element(poisson.cbegin(), poisson.cend())),
a_(a), E_(E), G_(G), poisson_(poisson)
Expand All @@ -159,6 +158,9 @@ OrthotropicSolid::OrthotropicSolid(Real rho_0, std::array<Vecd, Dimensions> a, s
CalculateA0();
CalculateAllMu();
CalculateAllLambda();
// we take the max. K in three pinciple directions to approximate the maximum of
// the Bulk modulus --> for time step size calculation
reset_sound_speeds();
};
//=================================================================================================//
Matd OrthotropicSolid::StressPK2(Matd &F, size_t index_i)
Expand All @@ -175,11 +177,17 @@ Matd OrthotropicSolid::StressPK2(Matd &F, size_t index_i)
Summa2 += Lambda_(i, j) * (CalculateBiDotProduct(A_[i], strain) * A_[j] +
CalculateBiDotProduct(A_[j], strain) * A_[i]);
}
stress_PK2 += Mu_[i] * (((A_[i] * strain) + (strain * A_[i])) + 1 / 2 * (Summa2));
stress_PK2 += Mu_[i] * (A_[i] * strain + strain * A_[i]) + 0.5 * Summa2;
}
return stress_PK2;
}
//=================================================================================================//
Matd StressCauchy(Matd &almansi_strain, size_t particle_index_i)
{
// TODO: implement
throw(std::runtime_error("Cauchy stress of orthotropic material not yet implemented"));
}
//=================================================================================================//
Real OrthotropicSolid::VolumetricKirchhoff(Real J)
{
return K0_ * J * (J - 1);
Expand All @@ -191,6 +199,18 @@ void OrthotropicSolid::CalculateA0()
A_[i] = a_[i] * a_[i].transpose();
}
//=================================================================================================//
void OrthotropicSolid::reset_sound_speeds()
{
std::array<Real, Dimensions> K;
for (int i = 0; i < Dimensions; ++i)
K[i] = (Lambda_.row(i).sum() + 2.0 * Mu_[i]) / 3.0;
K0_ = *std::max_element(K.cbegin(), K.cend());
G0_ = *std::max_element(G_.cbegin(), G_.cend());
lambda0_ = Lambda_.maxCoeff();
setSoundSpeeds();
setContactStiffness(c0_);
}
//=================================================================================================//
Matd FeneNeoHookeanSolid::StressPK2(Matd &F, size_t index_i)
{
Matd right_cauchy = F.transpose() * F;
Expand Down
24 changes: 13 additions & 11 deletions src/shared/materials/elastic_solid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -213,10 +213,10 @@ class NeoHookeanSolidIncompressible : public LinearElasticSolid
* @class OrthotropicSolid
* @brief Generic definition with 3 orthogonal directions + 9 independent parameters,
* ONLY for 3D applications
* @param "a" --> 3 principal direction vectors
* @param "E" --> 3 principal Young's moduli
* @param "G" --> 3 principal shear moduli
* @param "poisson" --> 3 principal Poisson's ratios
* @param "a" --> 3 principal direction vectors in the order of a_x, a_y, a_z
* @param "E" --> 3 principal Young's moduli in the order of E_x, E_y, E_z
* @param "G" --> 3 principal shear moduli in the order of G_xy, G_yz, G_zx
* @param "poisson" --> 3 principal Poisson's ratios in the order of nu_xy, ny_zx, nu_yz
*/
class OrthotropicSolid : public LinearElasticSolid
{
Expand All @@ -225,9 +225,10 @@ class OrthotropicSolid : public LinearElasticSolid
std::array<Real, Dimensions> G, std::array<Real, Dimensions> poisson);

/** second Piola-Kirchhoff stress related with green-lagrangian deformation tensor */
virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
Matd StressPK2(Matd &deformation, size_t particle_index_i) override;
Matd StressCauchy(Matd &almansi_strain, size_t particle_index_i) override;
/** Volumetric Kirchhoff stress determinate */
virtual Real VolumetricKirchhoff(Real J) override;
Real VolumetricKirchhoff(Real J) override;

protected:
// input data
Expand All @@ -236,13 +237,14 @@ class OrthotropicSolid : public LinearElasticSolid
std::array<Real, Dimensions> G_;
std::array<Real, Dimensions> poisson_;
// calculated data
Real Mu_[Dimensions];
std::array<Real, Dimensions> Mu_;
Matd Lambda_;
Matd A_[Dimensions];
std::array<Matd, Dimensions> A_;

virtual void CalculateAllMu();
virtual void CalculateAllLambda();
virtual void CalculateA0();
void CalculateAllMu();
void CalculateAllLambda();
void CalculateA0();
void reset_sound_speeds();
};

/**
Expand Down