[cig-commits] r6954 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu May 24 09:56:34 PDT 2007
Author: willic3
Date: 2007-05-24 09:56:27 -0700 (Thu, 24 May 2007)
New Revision: 6954
Modified:
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
Log:
Finished putting in stuff for calStress (not tested yet).
I still need to put some of the time-dependent stress computation stuff
in a separate routine, put in an updateState function, and make use of
logical flags.
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-05-24 14:27:53 UTC (rev 6953)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-05-24 16:56:27 UTC (rev 6954)
@@ -209,9 +209,9 @@
const double e23 = totalStrain[4];
const double e13 = totalStrain[5];
- const double etraceTpdt = e11 + e22 + e33;
- const double emeanTpdt = etraceTpdt/3.0;
- const double s123 = lambda * etraceTpdt;
+ const double traceStrainTpdt = e11 + e22 + e33;
+ const double meanStrainTpdt = traceStrainTpdt/3.0;
+ const double s123 = lambda * traceStrainTpdt;
if (useElasticBehavior()) {
(*stress)[0] = s123 + 2.0*mu*e11;
@@ -222,7 +222,11 @@
(*stress)[5] = 2.0 * mu * e13;
} else {
const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 };
- const double smeanTpdt = bulkmodulus * etraceTpdt;
+ const double meanStressTpdt = bulkmodulus * traceStrainTpdt;
+ const double meanStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][0] +
+ parameters[_MaxwellIsotropic3D::pidStrainT][1] +
+ parameters[_MaxwellIsotropic3D::pidStrainT][2];
+
// The code below should probably be in a separate function since it
// is used more than once. I should also probably cover the possibility
// that Maxwell time is zero (although this should never happen).
@@ -234,14 +238,30 @@
double factorial = 1.0;
double fraction = 1.0;
dq = 1.0;
- for (int iTerm=0; iTerm < numTerms; ++iTerm) {
+ for (int iTerm=2; iTerm <= numTerms; ++iTerm) {
factorial *= iTerm;
fSign *= -1.0;
fraction *= _dt/maxwelltime;
dq += fSign*fraction/factorial;
} // for
} else
- dq = maxwelltime*(1.0-exp(-dt/maxwelltime))/dt;
+ dq = maxwelltime*(1.0-exp(-_dt/maxwelltime))/_dt;
+ const double expFac = exp(-_dt/maxwelltime);
+ const double elasFac = 2.0*mu;
+ double devStrainTpdt = 0.0;
+ double devStrainT = 0.0;
+ double devStressTpdt = 0.0;
+ for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp) {
+ devStrainTpdt = totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
+ devStrainT = parameters[_MaxwellIsotropic3D::pidStrainT][iComp] -
+ diag[iComp]*meanStrainT;
+ visStrain = expFac*parameters[_MaxwellIsotropic3D::pidVisStrain][iComp] +
+ dq*(devStrainTpdt - devStrainT);
+ devStressTpdt = elasFac*visStrain;
+ // Later I will want to put in initial stresses.
+ (*stress)[iComp] =diag[iComp]*smeanTpdt+devStressTpdt;
+ } // for
+ } //else
} // _calcStress
// ----------------------------------------------------------------------
More information about the cig-commits
mailing list