#include "modelDeemo.h"
#include "version.txt"
// excerpt-export-start
int __stdcall DllMain(void *,unsigned, void *) {
return 1;
}
extern "C" EXPORT_TAG const char *getName() {
#ifdef MODELDEBUG
return "cmodelDeemod";
#else
return "cmodelDeemo";
#endif
}
extern "C" EXPORT_TAG unsigned getMajorVersion() {
return MAJOR_VERSION;
}
extern "C" EXPORT_TAG unsigned getMinorVersion() {
return UPDATE_VERSION;
}
extern "C" EXPORT_TAG void *createInstance() {
models::ModelDeemo *m = new models::ModelDeemo();
return (void *)m;
}
// excerpt-export-end
namespace models {
// excerpt-con-start
ModelDeemo::ModelDeemo() : bulk_(0.0), Mshear_(0.0), Mviscosity1_(0.0), Kshear1_(0.0), Kviscosity1_(0.0), Kshear2_(0.0), Kviscosity2_(0.0), Mviscosity2_(0.0), A_(0.0) {
//kelvinStrain1
Mekd_[0] = 0.0;
Mekd_[1] = 0.0;
Mekd_[2] = 0.0;
Mekd_[3] = 0.0;
Mekd_[4] = 0.0;
Mekd_[5] = 0.0;
//kelvinStrain2
Mekd_[6] = 0.0;
Mekd_[7] = 0.0;
Mekd_[8] = 0.0;
Mekd_[9] = 0.0;
Mekd_[10] = 0.0;
Mekd_[11] = 0.0;
}
// excerpt-con-end
String ModelDeemo::getName() const {
#ifdef MODELDEBUG
return L"Deemo-debug";
#else
return L"Deemo";
#endif
}
String ModelDeemo::getFullName() const {
#ifdef MODELDEBUG
return L"Deemo Debug";
#else
return L"Deemo";
#endif
}
UInt ModelDeemo::getMinorVersion() const {
return UPDATE_VERSION;
}
String ModelDeemo::getProperties() const { //获取参数
return L"bulk_, Kshear1_, Mshear_, Kviscosity1_, Mviscosity1_,"//burgers的参数
L"Kshear2_, Kviscosity2_, Mviscosity2_, A_," //塑性部分的参数
L"strain-kelvin-xx,strain-kelvin-yy,strain-kelvin-zz,strain-kelvin-xy,strain-kelvin-xz,strain-kelvin-yz,"//黏弹性体应变
L"cohesion,friction,dilation,tension,strain-shear-plastic,strain-tensile-plastic -strain-tension-plastic";//MC准则
}
String ModelDeemo::getStates() const {
return L"shear-n,tension-n,shear-p,tension-p";
}
Variant ModelDeemo::getProperty(UInt index) const {
switch (index) {
case 1: return bulk_;
case 2: return Kshear1_;
case 3: return Mshear_;
case 4: return Kviscosity1_;
case 5: return Mviscosity1_;
case 6: return Kshear2_;
case 7: return Kviscosity2_;
case 8: return Mviscosity2_;
case 9: return A_;//模型参数
case 10: return Mekd_[0];
case 11: return Mekd_[1];
case 12: return Mekd_[2];
case 13: return Mekd_[3];
case 14: return Mekd_[4];
case 15: return Mekd_[5];
case 16: return Mekd_[6];
case 17: return Mekd_[7];
case 18: return Mekd_[8];
case 19: return Mekd_[9];
case 20: return Mekd_[10];
case 21: return Mekd_[11];
case 22: return cohesion_;
case 23: return friction_;
case 24: return dilation_;
case 25: return tension_;
case 26: return sHP_;//剪切应变
case 27: return tHP_;//张拉应变
}
return(0.0);
}
void ModelDeemo::setProperty(UInt index,const Variant &p,UInt restoreVersion) {
ConstitutiveModel::setProperty(index, p, restoreVersion);
switch (index) {
case 1: bulk_ = p.toDouble(); break;
case 2: Kshear1_ = p.toDouble(); break;
case 3: Mshear_ = p.toDouble(); break;
case 4: Kviscosity1_ = p.toDouble(); break;
case 5: Mviscosity1_ = p.toDouble(); break;
case 6: Kshear2_ = p.toDouble(); break;
case 7: Kviscosity2_ = p.toDouble(); break;
case 8: Mviscosity2_ = p.toDouble(); break;
case 9: A_ = p.toDouble(); break;
case 10: Mekd_[0] = p.toDouble(); break;
case 11: Mekd_[1] = p.toDouble(); break;
case 12: Mekd_[2] = p.toDouble(); break;
case 13: Mekd_[3] = p.toDouble(); break;
case 14: Mekd_[4] = p.toDouble(); break;
case 15: Mekd_[5] = p.toDouble(); break;
case 16: Mekd_[6] = p.toDouble(); break;
case 17: Mekd_[7] = p.toDouble(); break;
case 18: Mekd_[8] = p.toDouble(); break;
case 19: Mekd_[9] = p.toDouble(); break;
case 20: Mekd_[10] = p.toDouble(); break;
case 21: Mekd_[11] = p.toDouble(); break;
case 22: cohesion_ = p.toDouble(); break;
case 23: friction_ = p.toDouble(); break;
case 24: dilation_ = p.toDouble(); break;
case 25: tension_ = p.toDouble(); break;
case 26: sHP_ = p.toDouble(); break;
case 27: tHP_ = p.toDouble(); break;
}
}
/*bool ModelDeemo::isPropertyAdvanced(UInt i) const {
if (i <= 4)
return ModelDeemo::isPropertyAdvanced(i);
else if (i==9)
return true;
return false;
}*/
//下面这段cvisc中是vm,这里是mm
void ModelDeemo::copy(const ConstitutiveModel *m) {
const ModelDeemo *mm = dynamic_cast<const ModelDeemo *>(m);
if (!mm) throw std::runtime_error("Internal error: constitutive model dynamic cast failed.");
//
ConstitutiveModel::copy(m);
//
bulk_ = mm->bulk_;
Kshear1_ = mm->Kshear1_;
Mshear_ = mm->Mshear_;
Kviscosity1_ = mm->Kviscosity1_;
Mviscosity1_ = mm->Mviscosity1_;
Kshear2_ = mm->Kshear2_;
Kviscosity2_ = mm->Kviscosity2_;
Mviscosity2_ = mm->Mviscosity2_;
A_ = mm->A_;
Mekd_[0] = mm->Mekd_[0];
Mekd_[1] = mm->Mekd_[1];
Mekd_[2] = mm->Mekd_[2];
Mekd_[3] = mm->Mekd_[3];
Mekd_[4] = mm->Mekd_[4];
Mekd_[5] = mm->Mekd_[5];
Mekd_[6] = mm->Mekd_[6];
Mekd_[7] = mm->Mekd_[7];
Mekd_[8] = mm->Mekd_[8];
Mekd_[9] = mm->Mekd_[9];
Mekd_[10] = mm->Mekd_[10];
Mekd_[11] = mm->Mekd_[11];
cohesion_ = mm->cohesion_;
friction_ = mm->friction_;
dilation_ = mm->dilation_;
tension_ = mm->tension_;
sHP_ = mm->sHP_;
tHP_ = mm->tHP_;
}
// excerpt-run-start参数初始化,计算模型中间参数值;
void ModelDeemo::initialize(UByte d,State *s) {
ConstitutiveModel::initialize(d,s);
//参数非零化
if (Mshear_ <= 0.0) Mshear_ = 1e-20;
if (Kshear1_ <= 0.0) Kshear1_ = 0.0;
if (Kshear2_ <= 0.0) Kshear2_ = 0.0;
if (Kviscosity1_ <= 0.0) {
Kviscosity1_ = 0.0;
Kshear1_ = 0.0;
}
if (Kviscosity2_ <= 0.0) {
Kviscosity2_ = 0.0;
Kshear2_ = 0.0;
}
if (Mviscosity1_ <= 0.0) Mviscosity1_ = 0.0;
if (Mviscosity2_ <= 0.0) Mviscosity2_ = 0.0;
if (A_ <= 0.0) A_ = 0.0;
// 3. 初始化Mohr-Coulomb准则相关参数(参考CVisc模型)
Double rsin = std::sin(friction_ * degrad); // degrad是角度转弧度的系数
nph_ = (1.0 + rsin) / (1.0 - rsin); // 塑性硬化参数
csn_ = 2.0 * cohesion_ * sqrt(nph_); // 修正后的内聚力
// 限制抗拉强度不超过内聚力和摩擦角决定的顶点值
if (friction_ > 0.0) {
Double apex = cohesion_ / std::tan(friction_ * degrad);
tension_ = std::min(tension_, apex);
}
// 剪胀角相关参数
rsin = std::sin(dilation_ * degrad);
nps_ = (1.0 + rsin) / (1.0 - rsin); // 剪胀参数
// 计算应力比参数
rc_ = std::sqrt(1.0 + nph_ * nph_);
}
static const UInt Dqs = 12;
static const UInt Dqt = 13;
void ModelDeemo::run(UByte d,State *s) {
ConstitutiveModel::run(d,s);
// excerpt-state-start
if (s->state_ & shear_now) s->state_ |= shear_past;
s->state_ &= ~shear_now;
if (s->state_ & tension_now) s->state_ |= tension_past;
s->state_ &= ~tension_now;
// excerpt-state-end
UInt iPlas = 0;//未屈服(黏弹性)
// dEkd values now stored in s->working_[] array (necessary for thread safety)
Double tempk1 = 0, tempk2 = 0, tempm1 = 0, tempm2 = 0.0;
Double dCrtdel = (s->isCreep() ? s->getTimeStep() : 0.0);
Double dSubZoneVolume = s->getSubZoneVolume();
if (!s->sub_zone_) {
s->working_[0] = 0.0;
s->working_[1] = 0.0;
s->working_[2] = 0.0;
s->working_[3] = 0.0;
s->working_[4] = 0.0;
s->working_[5] = 0.0;
s->working_[6] = 0.0;
s->working_[7] = 0.0;
s->working_[8] = 0.0;
s->working_[9] = 0.0;
s->working_[10] = 0.0;
s->working_[11] = 0.0;
s->working_[Dqs] = 0.0;
s->working_[Dqt] = 0.0;
}
if (Kviscosity1_ <= 0.0) tempk1 = 0.0;
else tempk1 = 1.0 / Kviscosity1_;
//
if (Mviscosity1_ <= 0.0) tempm1 = 0.0;
else tempm1 = 1.0 / Mviscosity1_;
if (Kviscosity2_ <= 0.0) tempk2 = 0.0;
else tempk2 = 1.0 / Kviscosity2_;
if (Mviscosity2_ <= 0.0) tempm2 = 0.0;
else tempm2 = 1.0 / Mviscosity2_;
Double temp1 = 0.5 * Kshear1_ * dCrtdel * tempk1;
Double a1_con = 1.0 + temp1;//A1
Double b1_con = 1.0 - temp1;//B1
Double ba1 = b1_con / a1_con;
Double bac1 = ba1 - 1.0;//(B1/A1-1)
Double temp2 = 0.5 * Kshear2_ * dCrtdel * tempk2;
Double a2_con = 1.0 + temp2;//A2
Double b2_con = 1.0 - temp2;//B2
Double ba2 = b2_con / a2_con;
Double bac2 = ba2 - 1.0;//(B2/A2-1)
Double tempb = (tempm1 + tempk1 / a1_con+ tempk2 / a2_con) * dCrtdel * 0.25;
Double tempa = 0.5 / Mshear_;
Double x_con = tempa + tempb;//a
Double y_con = tempa - tempb;//b
Double z1_con = dCrtdel * tempk1 / (4.0 * a1_con);
Double z2_con = dCrtdel * tempk2 / (4.0 * a2_con);
//Double z_con = z1_con + z2_con;//
Double c1dxc = 1.0 / x_con;//1/a
//;--- partition strains ---
Double dev = s->stnE_.s11() + s->stnE_.s22() + s->stnE_.s33();
Double dev3 = d1d3 * dev;
Double de11d = s->stnE_.s11() - dev3;
Double de22d = s->stnE_.s22() - dev3;
Double de33d = s->stnE_.s33() - dev3;//应变计算
//;--- partition stresses---
Double s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33());
Double s11d = s->stnS_.s11() - s0;
Double s22d = s->stnS_.s22() - s0;
Double s33d = s->stnS_.s33() - s0;//应力计算
//;--- remember old stresses ---
Double s11old = s11d;
Double s22old = s22d;
Double s33old = s33d;
Double s12old = s->stnS_.s12();
Double s13old = s->stnS_.s13();
Double s23old = s->stnS_.s23();
//;--- new trial deviator stresses assuming viscoelastic increments ---
s11d = (de11d + s11d * y_con - Mekd_[0] * bac1) * c1dxc;
s22d = (de22d + s22d * y_con - Mekd_[1] * bac1) * c1dxc;
s33d = (de33d + s33d * y_con - Mekd_[2] * bac1) * c1dxc;
Double s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con - Mekd_[3] * bac1) * c1dxc;
Double s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con - Mekd_[4] * bac1) * c1dxc;
Double s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con - Mekd_[5] * bac1) * c1dxc;
//;--- new trial isotropic stress assuming elastic increment ---
s0 += bulk_ * dev;
Double s11i = s11d + s0;
Double s22i = s22d + s0;
Double s33i = s33d + s0;
//; --- trial stresses ---
s->stnS_.rs11() = s11i;
s->stnS_.rs22() = s22i;
s->stnS_.rs33() = s33i;
s->stnS_.rs12() = s12i;
s->stnS_.rs13() = s13i;
s->stnS_.rs23() = s23i;
//先标定了没有进入塑性的部分
//下面是在做塑性状态的判断
if (canFail()) {
SymTensorInfo info;
DVect3 prin = s->stnS_.getEigenInfo(&info);
Double e1_ = bulk_ + d2d3 * c1dxc;//α1
Double e2_ = bulk_ - d1d3 * c1dxc;//α2
Double ra = e1_ - nps_ * e2_;
//Double e21_ = e2_ / e1_;
Double rb = e2_ * (1.0 - nps_);
Double rd = e2_ - nps_ * e1_;//应力更新公式的λ的系数
//定义塑性部分相关参数
Double G3 = Kshear2_; // 开尔文体弹性模量 (对应Kshear2_)
Double A = A_; // 非线性黏壶参数
Double fs = -prin.x() + nph_ * prin.z() - csn_;
Double ftz = prin.z() - tension_;
Double fty = prin.y() - tension_;
Double ftx = prin.x() - tension_;//MC破坏准则
Double fsd = fs / rc_;//(归一化但是还没搞清楚是啥)
Double lambda_s_1 = fs / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel)); //λs1
Double lambda_s_2 = fs * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1); //λs2
Double lambda_s = lambda_s_1+ lambda_s_2;//λs
Double beta_s_xt = ftx / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel));//βx
Double gamma_s_xt = ftx * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);//γx
Double lambda_s_xt = beta_s_xt + gamma_s_xt;//λsxt
Double beta_s_yt = fty / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel));//βy
Double gamma_s_yt = fty * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);//γy
Double lambda_s_yt = beta_s_yt + gamma_s_yt;//λsyt
Double beta_s_zt = ftz / (2 * G3)*(1.0 - exp(-G3 * tempk2 * dCrtdel));//βz
Double gamma_s_zt = ftz * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);//γz
Double lambda_s_zt = beta_s_zt + gamma_s_zt;//λszt
if (iPlas == 0) {
// 塑性未激活:仅用第一个 Kelvin 元件(原逻辑)
s11d = (de11d + s11d * y_con - Mekd_[0] * bac1) * c1dxc;
s22d = (de22d + s22d * y_con - Mekd_[1] * bac1) * c1dxc;
s33d = (de33d + s33d * y_con - Mekd_[2] * bac1) * c1dxc;
s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con - Mekd_[3] * bac1) * c1dxc;
s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con - Mekd_[4] * bac1) * c1dxc;
s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con - Mekd_[5] * bac1) * c1dxc;
}
else {
//Double s11_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s11d)) / (2 * m1_);
//Double s22_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s22d)) / (2 * m1_);
//Double s33_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s33d)) / (2 * m1_);
//Double s12_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s->stnS_.s12())) / (2 *m1_);
//Double s13_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s->stnS_.s13())) / (2 *m1_);
//Double s23_2 = (-m2_ + std::sqrt((m2_)*(m2_)+4 * m1_*s->stnS_.s23())) / (2 *m1_);
//非线性牛顿体
Double s11_2 = s11d * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);
Double s22_2 = s22d * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);
Double s33_2 = s33d * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);
Double s12_2 = s->stnS_.s12() * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);
Double s13_2 = s->stnS_.s13() * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);
Double s23_2 = s->stnS_.s23() * tempm2 / (2 * A)*(exp(A* dCrtdel) - 1);
// 塑性已激活:使用两个 Kelvin 元件和一个非线性牛顿体
s11d = (de11d + s11d * y_con - s11_2 - Mekd_[0] * bac1 - Mekd_[6] * bac2) * c1dxc;
s22d = (de22d + s22d * y_con - s22_2 - Mekd_[1] * bac1 - Mekd_[7] * bac2) * c1dxc;
s33d = (de33d + s33d * y_con - s33_2 - Mekd_[2] * bac1 - Mekd_[8] * bac2) * c1dxc;
s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con - s12_2 - Mekd_[3] * bac1 - Mekd_[9] * bac2) * c1dxc;
s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con - s13_2 - Mekd_[4] * bac1 - Mekd_[10] * bac2) * c1dxc;
s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con - s23_2 - Mekd_[5] * bac1 - Mekd_[11] * bac2) * c1dxc;
}
if (fsd > 0.0 && fsd >= ftz) {//剪切塑性修正
iPlas = 1;
s->state_ |= shear_now;
prin.rx() += lambda_s * ra;
prin.ry() += lambda_s * rb;
prin.rz() += lambda_s * rd;//已检查
}
else if (ftz > 0.0 && ftz >= fsd) {
s->state_ |= tension_now;//拉伸塑性修正
if (ftx > 0.0) {
iPlas = 4;
prin.rx() -= lambda_s_xt * e1_ + (lambda_s_yt + lambda_s_zt) * e2_;
prin.ry() -= lambda_s_yt * e1_ + (lambda_s_xt + lambda_s_zt) * e2_;
prin.rz() -= lambda_s_zt * e1_ + (lambda_s_xt + lambda_s_yt) * e2_;//三个主应力失效
}
else if (fty > 0.0) {
iPlas = 3;
prin.rx() -= (lambda_s_yt + lambda_s_zt) * e2_;
prin.ry() -= (lambda_s_zt * e2_ + lambda_s_yt * e1_);
prin.rz() -= (lambda_s_zt * e1_ + lambda_s_yt * e2_);//两个主应力失效
}
else {
iPlas = 2;
//Double tco = ftz * e21_;//应该是用不上了
prin.rx() -= lambda_s_zt * e2_ ;
prin.ry() -= lambda_s_zt * e2_ ;
prin.rz() -= lambda_s_zt * e1_ ;//一个主应力失效
}
}
if (friction_ > 0.0) {
Double apex = cohesion_ / tan(friction_*degrad);
if (prin.x() >= apex || prin.y() >= apex || prin.z() >= apex) {//如果任意方向上的应力超过这个阈值,进入拉伸状态
iPlas = 4;
s->state_ |= tension_now;
prin.rx() = apex;
prin.ry() = apex;
prin.rz() = apex;
}
}
if (iPlas) {
s->stnS_ = info.resolve(prin);//用当前的prin应力来解析并更新应力状态
if (iPlas == 1) {
Double dDe1p = -lambda_s;//1方向的塑性应变增量
Double dDe3p = -1 * lambda_s * nps_;//3方向的塑性应变增量
Double dDepa = d1d3 * (dDe1p + dDe3p);
dDe1p -= dDepa;
dDe3p -= dDepa;
s->working_[Dqs] += sqrt(0.5 * (dDe1p*dDe1p + dDepa * dDepa + dDe3p * dDe3p)) * dSubZoneVolume;
}
if (iPlas == 2) {
Double dAux = lambda_s_zt;
s->working_[Dqt] += dAux * dSubZoneVolume;
}
if (iPlas == 3) {
Double dAux = lambda_s_zt + lambda_s_yt;
s->working_[Dqt] += dAux * dSubZoneVolume;
}
if (iPlas == 4) {
Double dAux = lambda_s_zt + lambda_s_yt + lambda_s_xt;
s->working_[Dqt] += dAux * dSubZoneVolume;
}
}
}
//--- sub-zone contribution to Kelvin1-strains ---
s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33());
s->working_[0] += (Mekd_[0] * ba1 + (s->stnS_.s11() - s0 + s11old) * z1_con) * s->getSubZoneVolume();
s->working_[1] += (Mekd_[1] * ba1 + (s->stnS_.s22() - s0 + s22old) * z1_con) * s->getSubZoneVolume();
s->working_[2] += (Mekd_[2] * ba1 + (s->stnS_.s33() - s0 + s33old) * z1_con) * s->getSubZoneVolume();
s->working_[3] += (Mekd_[3] * ba1 + (s->stnS_.s12() + s12old) * z1_con) * s->getSubZoneVolume();
s->working_[4] += (Mekd_[4] * ba1 + (s->stnS_.s13() + s13old) * z1_con) * s->getSubZoneVolume();
s->working_[5] += (Mekd_[5] * ba1 + (s->stnS_.s23() + s23old) * z1_con) * s->getSubZoneVolume();
//--- sub-zone contribution to Kelvin2-strains ---
s->working_[6] += (Mekd_[6] * ba2 + (s->stnS_.s11() - s0 + s11old) * z2_con) * s->getSubZoneVolume();
s->working_[7] += (Mekd_[7] * ba2 + (s->stnS_.s22() - s0 + s22old) * z2_con) * s->getSubZoneVolume();
s->working_[8] += (Mekd_[8] * ba2 + (s->stnS_.s33() - s0 + s33old) * z2_con) * s->getSubZoneVolume();
s->working_[9] += (Mekd_[9] * ba2 + (s->stnS_.s12() + s12old) * z2_con) * s->getSubZoneVolume();
s->working_[10] += (Mekd_[10] * ba2 + (s->stnS_.s13() + s13old) * z2_con) * s->getSubZoneVolume();
s->working_[11] += (Mekd_[11] * ba2 + (s->stnS_.s23() + s23old) * z2_con) * s->getSubZoneVolume();
//--- update stored Kelvin-strains and plastic strain ---
if (s->sub_zone_ == s->total_sub_zones_ - 1) {//这行代码检查是否在最后一个子区域内,如果是,则执行后续的更新操作
Double Aux = 1. / s->getZoneVolume();//Aux变量是区域体积的倒数
if (s->overlay_ == 2) Aux *= 0.5;//如果overlay_等于2,则Aux乘以0.5.
Mekd_[0] = s->working_[0] * Aux;
Mekd_[1] = s->working_[1] * Aux;
Mekd_[2] = s->working_[2] * Aux;
Mekd_[3] = s->working_[3] * Aux;
Mekd_[4] = s->working_[4] * Aux;
Mekd_[5] = s->working_[5] * Aux;
Mekd_[6] = s->working_[6] * Aux;
Mekd_[7] = s->working_[7] * Aux;
Mekd_[8] = s->working_[8] * Aux;
Mekd_[9] = s->working_[9] * Aux;
Mekd_[10] = s->working_[10] * Aux;
Mekd_[11] = s->working_[11] * Aux;
if (canFail()) {
sHP_ += s->working_[Dqs] * Aux;
tHP_ += s->working_[Dqt] * Aux;
}
}
// excerpt-run-end
}
} // namespace models
// EOF这个是我的本构,不管怎么跑都会进入加速蠕变阶段,下面是burgers-mohr的你对比来看看#include "modelcvisc.h"
#include "base/src/version.h"
namespace models {
ModelCVisc::ModelCVisc() : ModelBurger(),
cohesion_(0.0),friction_(0.0),dilation_(0.0),tension_(0.0),
sHP_(0.0),tHP_(0.0),nph_(0.0),csn_(0.0),rc_(0.0),nps_(0.0) {
}
UInt ModelCVisc::getMinorVersion() const{
return version::update;
}
String ModelCVisc::getProperties(void) const {
return ModelBurger::getProperties() +
L",cohesion,friction,dilation,tension,strain-shear-plastic,strain-tensile-plastic -strain-tension-plastic";
}
String ModelCVisc::getStates(void) const {
return L"shear-n,tension-n,shear-p,tension-p";
}
Variant ModelCVisc::getProperty(UInt index) const {
if (index <= 11)
return ModelBurger::getProperty(index);
else {
switch (index) {
case 12: return cohesion_;
case 13: return friction_;
case 14: return dilation_;
case 15: return tension_;
case 16: return sHP_;
case 17: return tHP_;
}
}
return(0.0);
}
void ModelCVisc::setProperty(UInt index, const Variant &p,UInt restoreVersion) {
if (index <= 11)
ModelBurger::setProperty(index,p,restoreVersion);
else {
switch (index) {
case 12: cohesion_ = p.toDouble(); break;
case 13: friction_ = p.toDouble(); break;
case 14: dilation_ = p.toDouble(); break;
case 15: tension_ = p.toDouble(); break;
case 16: sHP_ = p.toDouble(); break;
case 17: tHP_ = p.toDouble(); break;
}
}
}
void ModelCVisc::copy(const ConstitutiveModel *m) {
const ModelCVisc *vm = dynamic_cast<const ModelCVisc *>(m);
if (!vm) throw std::runtime_error("Internal error: constitutive model dynamic cast failed.");
ModelBurger::copy(m);
cohesion_ = vm->cohesion_;
friction_ = vm->friction_;
dilation_ = vm->dilation_;
tension_ = vm->tension_;
sHP_ = vm->sHP_;
tHP_ = vm->tHP_;
}
void ModelCVisc::initialize(UByte d,State *s) {
ModelBurger::initialize(d, s);
Double rsin = std::sin(friction_ * degrad);
nph_ = (1.0 + rsin) / (1.0 - rsin);
csn_ = 2.0 * cohesion_ * sqrt(nph_);
if (friction_) {
Double apex = cohesion_ / std::tan(friction_ * degrad);
tension_ = std::min(tension_,apex);
}
rsin = std::sin(dilation_ * degrad);
nps_ = (1.0 + rsin) / (1.0 - rsin);
rc_ = std::sqrt(1.0 + nph_*nph_);
}
static const UInt Dqs = 6;
static const UInt Dqt = 7;
void ModelCVisc::run(UByte d,State *s) {
ConstitutiveModel::run(d,s);
if (s->state_ & shear_now) s->state_ |= shear_past;
s->state_ &= ~shear_now;
if (s->state_ & tension_now) s->state_ |= tension_past;
s->state_ &= ~tension_now;
UInt iPlas = 0;
Double tempk=0, tempm=0;
Double dCrtdel = (s->isCreep() ? s->getTimeStep() : 0.0);
Double dSubZoneVolume = s->getSubZoneVolume();
if (!s->sub_zone_) {
s->working_[0] = 0.0;
s->working_[1] = 0.0;
s->working_[2] = 0.0;
s->working_[3] = 0.0;
s->working_[4] = 0.0;
s->working_[5] = 0.0;
s->working_[Dqs] =0.0;
s->working_[Dqt] =0.0;
}
if(Kviscosity_ <= 0.0) tempk = 0.0;
else tempk = 1.0 / Kviscosity_ ;
//
if (Mviscosity_ <= 0.0) tempm = 0.0;
else tempm = 1.0 / Mviscosity_ ;
Double temp = 0.5 * Kshear_ * dCrtdel * tempk;
Double a_con = 1.0 + temp;
Double b_con = 1.0 - temp;
Double ba = b_con / a_con;
Double bal = ba - 1.0;
temp = (tempm + tempk / a_con) * dCrtdel * 0.25 ;
Double temp1 = 0.5 / Mshear_;
Double x_con = temp1 + temp;
Double y_con = temp1 - temp;
Double z_con = dCrtdel * tempk / (4.0 * a_con);
Double c1dxc = 1.0 / x_con;
//;--- partition strains ---
Double dev = s->stnE_.s11() + s->stnE_.s22() + s->stnE_.s33() ;
Double dev3 = d1d3 * dev;
Double de11d = s->stnE_.s11() - dev3;
Double de22d = s->stnE_.s22() - dev3;
Double de33d = s->stnE_.s33() - dev3;
//;--- partition stresses---
Double s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33());
Double s11d = s->stnS_.s11() - s0;
Double s22d = s->stnS_.s22() - s0;
Double s33d = s->stnS_.s33() - s0;
//;--- remember old stresses ---
Double s11old = s11d;
Double s22old = s22d;
Double s33old = s33d;
Double s12old = s->stnS_.s12();
Double s13old = s->stnS_.s13();
Double s23old = s->stnS_.s23();
//;--- new trial deviator stresses assuming viscoelastic increments ---
s11d = (de11d + s11d * y_con - Mekd_[0] * bal) * c1dxc ;
s22d = (de22d + s22d * y_con - Mekd_[1]* bal) * c1dxc ;
s33d = (de33d + s33d * y_con - Mekd_[2]* bal) * c1dxc ;
Double s12i = (s->stnE_.s12() + s->stnS_.s12() * y_con - Mekd_[3] * bal) * c1dxc ;
Double s13i = (s->stnE_.s13() + s->stnS_.s13() * y_con - Mekd_[4] * bal) * c1dxc ;
Double s23i = (s->stnE_.s23() + s->stnS_.s23() * y_con - Mekd_[5] * bal) * c1dxc ;
//;--- new trial isotropic stress assuming elastic increment ---
s0 += bulk_ * dev;
Double s11i = s11d + s0;
Double s22i = s22d + s0;
Double s33i = s33d + s0;
//; --- trial stresses ---
s->stnS_.rs11() = s11i;
s->stnS_.rs22() = s22i;
s->stnS_.rs33() = s33i;
s->stnS_.rs12() = s12i;
s->stnS_.rs13() = s13i;
s->stnS_.rs23() = s23i;
if (canFail()) {
SymTensorInfo info;
DVect3 prin = s->stnS_.getEigenInfo(&info);
Double e1_ = bulk_ + d2d3*c1dxc;
Double e2_ = bulk_ - d1d3*c1dxc;
Double ra = e1_ - nps_*e2_;
Double rb = e2_ - nps_*e1_;
Double rd = ra - rb*nph_;
Double sc1_ = ra / rd;
Double sc3_ = rb / rd;
Double sc2_ = e2_*(1.0 - nps_)/rd;
Double e21_ = e2_/e1_;
Double sf1_ = 1.0/rd;
Double sf3_ = - nps_/rd;
Double fs = - prin.x() + nph_ * prin.z() - csn_;
Double fsd = fs / rc_;
Double ftz = prin.z() - tension_;
Double fty = prin.y() - tension_;
Double ftx = prin.x() - tension_;
if (fsd > 0.0 && fsd >= ftz) {
iPlas = 1;
s->state_ |= shear_now;
prin.rx() += fs * sc1_;
prin.ry() += fs * sc2_;
prin.rz() += fs * sc3_;
}
else if (ftz > 0.0 && ftz >= fsd) {
s->state_ |= tension_now;
if (ftx > 0.0) {
iPlas = 4;
prin.rx() = tension_;
prin.ry() = tension_;
prin.rz() = tension_;
} else if (fty > 0.0) {
iPlas = 3;
prin.rx() -= (fty + ftz) * e2_/(e1_+e2_);
prin.ry() = tension_;
prin.rz() = tension_;
} else {
iPlas = 2;
Double tco = ftz * e21_;
prin.rx() -= tco;
prin.ry() -= tco;
prin.rz() = tension_;
}
}
if (friction_ > 0.0) {
Double apex = cohesion_ / tan(friction_*degrad);
if (prin.x()>=apex || prin.y()>=apex || prin.z()>=apex) {
iPlas = 4;
s->state_ |= tension_now;
prin.rx() = apex;
prin.ry() = apex;
prin.rz() = apex;
}
}
if (iPlas) {
s->stnS_ = info.resolve(prin);
if (iPlas == 1) {
Double dDe1p = fs * sf1_;
Double dDe3p = fs * sf3_;
Double dDepa = d1d3 * (dDe1p + dDe3p);
dDe1p -= dDepa;
dDe3p -= dDepa;
s->working_[Dqs] += sqrt(0.5 * (dDe1p*dDe1p+dDepa*dDepa+dDe3p*dDe3p)) * dSubZoneVolume;
}
if (iPlas == 2) {
Double dAux = ftz / e1_;
s->working_[Dqt] += dAux * dSubZoneVolume;
}
if (iPlas == 3) {
Double dAux = (ftz+fty) / (e1_+e2_);
s->working_[Dqt] += dAux * dSubZoneVolume;
}
if (iPlas == 4) {
Double dAux = (ftz+fty+ftx) / bulk_;
s->working_[Dqt] += dAux * dSubZoneVolume;
}
}
}
//--- sub-zone contribution to Kelvin-strains ---
s0 = d1d3 * (s->stnS_.s11() + s->stnS_.s22() + s->stnS_.s33());
s->working_[0] += (Mekd_[0] * ba + (s->stnS_.s11() - s0 + s11old) * z_con) * s->getSubZoneVolume();
s->working_[1] += (Mekd_[1] * ba + (s->stnS_.s22() - s0 + s22old) * z_con) * s->getSubZoneVolume();
s->working_[2] += (Mekd_[2] * ba + (s->stnS_.s33() - s0 + s33old) * z_con) * s->getSubZoneVolume();
s->working_[3] += (Mekd_[3] * ba + (s->stnS_.s12() + s12old) * z_con) * s->getSubZoneVolume();
s->working_[4] += (Mekd_[4] * ba + (s->stnS_.s13() + s13old) * z_con) * s->getSubZoneVolume();
s->working_[5] += (Mekd_[5] * ba + (s->stnS_.s23() + s23old) * z_con) * s->getSubZoneVolume();
//--- update stored Kelvin-strains and plastic strain ---
if (s->sub_zone_ == s->total_sub_zones_-1) {
Double Aux = 1./s->getZoneVolume();
if (s->overlay_==2) Aux *= 0.5;
Mekd_[0]= s->working_[0] * Aux;
Mekd_[1]= s->working_[1] * Aux;
Mekd_[2]= s->working_[2] * Aux;
Mekd_[3]= s->working_[3] * Aux;
Mekd_[4]= s->working_[4] * Aux;
Mekd_[5]= s->working_[5] * Aux;
if (canFail()) {
sHP_ += s->working_[Dqs]*Aux;
tHP_ += s->working_[Dqt]*Aux;
}
}
}
}
// EOF
最新发布