[cig-commits] r11844 - in long/3D/Gale/trunk: . input/benchmarks src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Tue Apr 22 11:23:03 PDT 2008
Author: walter
Date: 2008-04-22 11:23:03 -0700 (Tue, 22 Apr 2008)
New Revision: 11844
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/input/benchmarks/extension.xml
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
Log:
r2104 at earth: boo | 2008-04-22 11:23:15 -0700
Clean up MohrCoulomb and make it correct in 3D
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2103
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2104
Modified: long/3D/Gale/trunk/input/benchmarks/extension.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/extension.xml 2008-04-22 18:23:00 UTC (rev 11843)
+++ long/3D/Gale/trunk/input/benchmarks/extension.xml 2008-04-22 18:23:03 UTC (rev 11844)
@@ -345,11 +345,9 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
-<!-- <param name="Type">DruckerPrager</param> -->
<param name="Type">MohrCoulomb</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
- <param name="StrainRateField">StrainRateField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
<param name="Context">context</param>
<param name="StrainWeakening">strainWeakening</param>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-04-22 18:23:00 UTC (rev 11843)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-04-22 18:23:03 UTC (rev 11844)
@@ -111,7 +111,6 @@
MohrCoulomb* self,
FeVariable* pressureField,
FeVariable* velocityGradientsField,
- FeVariable* strainRateField,
MaterialPointsSwarm* materialPointsSwarm,
double cohesion,
double cohesionAfterSoftening,
@@ -119,13 +118,9 @@
double frictionCoefficientAfterSoftening,
double minimumYieldStress)
{
- StandardParticle materialPoint;
- Dimension_Index dim = materialPointsSwarm->dim;
-
self->materialPointsSwarm = materialPointsSwarm;
self->pressureField = pressureField;
self->velocityGradientsField = velocityGradientsField;
- self->strainRateField = strainRateField;
self->cohesion = cohesion;
self->frictionCoefficient = frictionCoefficient;
@@ -167,7 +162,6 @@
FeVariable* pressureField;
MaterialPointsSwarm* materialPointsSwarm;
FeVariable* velocityGradientsField;
- FeVariable* strainRateField;
/* Construct Parent */
_YieldRheology_Construct( self, cf, data );
@@ -185,14 +179,11 @@
"PressureField", FeVariable, True, data );
velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
"VelocityGradientsField", FeVariable, True, data );
- strainRateField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "StrainRateField", FeVariable, True, data );
_MohrCoulomb_Init(
self,
pressureField,
velocityGradientsField,
- strainRateField,
materialPointsSwarm,
Stg_ComponentFactory_GetDouble( cf, self->name, "cohesion", 0.0 ),
Stg_ComponentFactory_GetDouble( cf, self->name, "cohesionAfterSoftening", 0.0 ),
@@ -210,15 +201,8 @@
void _MohrCoulomb_Initialise( void* rheology, void* data ) {
MohrCoulomb* self = (MohrCoulomb*) rheology;
- double* ptr;
Particle_Index lParticle_I;
Particle_Index particleLocalCount;
- double normalLength2;
- double invNormalLength;
- double initialDamageFraction;
- MaterialPoint* materialPoint;
- Index dof_I;
- Dimension_Index dim = self->materialPointsSwarm->dim;
AbstractContext* context = (AbstractContext*)data;
_YieldRheology_Initialise( self, data );
@@ -248,9 +232,6 @@
Coord xi )
{
MohrCoulomb* self = (MohrCoulomb*) rheology;
- double postFailureWeakening;
- double yieldCriterion;
- double yieldIndicator; /* A materialPoint will yield if yieldCriterion < yieldIndicator */
/* Don't want to yield on the first ever solve */
if ( constitutiveMatrix->previousSolutionExists == False ) {
@@ -272,48 +253,33 @@
Coord xi )
{
MohrCoulomb* self = (MohrCoulomb*) rheology;
- double* stress = self->currentStress;
Eigenvector* eigenvectorList = self->currentEigenvectorList;
double sigma_ns;
Dimension_Index dim = constitutiveMatrix->dim;
double stressMin;
double stressMax;
- double theta;
- /* Failure criterion from stress field - we can calculate
- * whether the material has failed at the current stress from
- * the principle stresses */
-
- theta = 0.5 * atan( 1.0/ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint ) );
+ /* The yield indicator is simply the difference in principal
+ stresses. In 2D, this is the 2nd invariant of the stress
+ tensor. This is not the actual shear stress. Rather, it is
+ the terms with deviatoric stress, because that scales
+ linearly with viscosity. So if we scale the viscosity, we
+ only affect the yield indicator. */
stressMin = eigenvectorList[0].eigenvalue;
stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
-
sigma_ns = 0.5 * ( stressMax - stressMin );
-
-/* printf("theta %lf\n",theta); */
- return fabs(sigma_ns);
+ if(dim==3)
+ {
+ double alpha=_MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint );
+ double theta = 0.5 * atan( 1.0/ alpha);
+ double factor=1/(sin(2*theta) - alpha*cos(2*theta));
- {
-/* SymmetricTensor stress; */
-/* SymmetricTensor strainRate; */
-/* double stressInv; */
+ sigma_ns-=(stressMax+stressMin)*0.5*alpha*factor;
+ }
-/* /\* Get Strain Rate *\/ */
-/* FeVariable_InterpolateWithinElement( self->strainRateField, lElement_I, xi, strainRate ); */
-
-/* /\* Get Stress *\/ */
-/* ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, stress ); */
-/* stressInv = SymmetricTensor_2ndInvariant( stress, constitutiveMatrix->dim ); */
-
-/* printf("indicator %lf %lf\n",fabs(sigma_ns),stressInv); */
-
-/* return stressInv; */
- }
-
-/* return fabs(sigma_ns); */
-
+ return fabs(sigma_ns);
}
double _MohrCoulomb_GetYieldCriterion(
@@ -329,7 +295,6 @@
double effectiveCohesion;
double effectiveFrictionCoefficient;
double frictionalStrength;
- double sigma_nn;
FeVariable* pressureField = self->pressureField;
double pressure;
@@ -339,35 +304,37 @@
FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
- /* Calculate frictional strength */
+ /* Calculate frictional strength. We modify the friction and
+ cohesion because we have grouped terms from the normal
+ stresses and moved it to the yield indicator. */
effectiveFrictionCoefficient = _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint );
effectiveCohesion = _MohrCoulomb_EffectiveCohesion( self, materialPoint );
-/* sigma_nn = _MohrCoulomb_Sigma_nn( self, materialPoint, constitutiveMatrix->dim ); */
-/* frictionalStrength = effectiveFrictionCoefficient * sigma_nn + effectiveCohesion; */
-
-
theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
factor=1/(sin(2*theta) - effectiveFrictionCoefficient*cos(2*theta));
-
frictionalStrength = effectiveFrictionCoefficient * pressure*factor + effectiveCohesion*factor;
-/* printf("criterion %lf %f\n",sigma_nn,pressure); */
-
- /* Make sure frictionalStrength is above the minimum */
+ /* If the minimumYieldStress is not set, then use the
+ effective cohesion. Maybe it should be the modified
+ effective cohesion, though that probably should not matter
+ much. */
minimumYieldStress = self->minimumYieldStress;
if(minimumYieldStress==0.0)
minimumYieldStress=effectiveCohesion;
+ /* Make sure frictionalStrength is above the minimum */
if ( frictionalStrength < minimumYieldStress)
frictionalStrength = minimumYieldStress;
return frictionalStrength;
}
-/* This sets the viscosity to yield_stress/(2*second stress invariant) if
- the material has failed. */
+/* This scales the viscosity so that the stress is the yield stress.
+ We can simply scale the viscosity by yieldCriterion/yieldIndicator
+ because yieldCriterion does not depend on the deviatoric stress
+ (and thus the viscosity), while the yieldIndicator depends linearly
+ on the deviatoric stresses. */
void _MohrCoulomb_HasYielded(
void* rheology,
@@ -378,46 +345,11 @@
double yieldCriterion,
double yieldIndicator )
{
- MohrCoulomb* self = (MohrCoulomb*)rheology;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- Dimension_Index dim = constitutiveMatrix->dim;
double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
-
-/* double frictionCoefficient = _MohrCoulomb_EffectiveFrictionCoefficient(self,materialPoint); */
-/* double cohesion = _MohrCoulomb_EffectiveCohesion(self,materialPoint); */
-/* double theta = 0.5 * atan( 1.0/ frictionCoefficient); */
-/* double beta = 1.0 - (yieldCriterion - frictionCoefficient*yieldIndicator/tan(2*theta)) */
-/* /(yieldIndicator*(1-frictionCoefficient/tan(2*theta))); */
-
-/* double stressMin = eigenvectorList[0].eigenvalue; */
-/* double stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue); */
-/* double pressure = self->currentPressure; */
-
-/* if(yieldCriterion==self->minimumYieldStress */
-/* || (yieldCriterion>=cohesion*0.99999 && yieldCriterion<=cohesion*1.000001)) */
-/* beta=1.0-yieldCriterion/yieldIndicator; */
-
double beta=1.0-yieldCriterion/(yieldIndicator);
-
-/* printf("hasyielded %lf %lf %lf %lf %lf %lf %lf %lf\n",viscosity,viscosity*(1-beta), */
-/* frictionCoefficient*pressure+cohesion, */
-/* (stressMin-stressMax)/2, yieldIndicator, yieldCriterion, */
-/* beta, 1-(frictionCoefficient*pressure+cohesion)/((stressMin-stressMax)/2)); */
-
- ConstitutiveMatrix_IsotropicCorrection( constitutiveMatrix, -viscosity * beta );
-
-/* printf("yielded %lf %lf %lf %lf %lf\n", */
-/* viscosity, yieldCriterion, */
-/* yieldIndicator, viscosity*beta, */
-/* ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix )); */
-
-/* viscosity = viscosity * yieldCriterion/yieldIndicator; */
-
-/* viscosity = yieldCriterion/(2*self->strainRateSecondInvariant); */
-
-/* ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity ); */
+ ConstitutiveMatrix_IsotropicCorrection( constitutiveMatrix, -viscosity * beta );
}
double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) {
@@ -444,34 +376,6 @@
return effectiveFrictionCoefficient;
}
-double _MohrCoulomb_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double pressure = self->currentPressure;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double tau_nn;
- double sigma_nn;
-
- double stressMin = eigenvectorList[0].eigenvalue;
- double stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
- double theta = 0.5 * atan( 1.0/ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint ) );
-
- tau_nn =cos( 2.0 * theta ) * ( stressMax - stressMin )*0.5;
-/* tau_nn = 0.5 * (( stressMax + stressMin ) + cos( 2.0 * theta ) * ( stressMax - stressMin )); */
-
-/* sigma_nn = -tau_nn; */
-/* sigma_nn = pressure - tau_nn; */
-/* sigma_nn = pressure; */
- sigma_nn = pressure + tau_nn;
-
-/* printf("sigma_nn %lf %lf %lf %lf %lf\n",pressure,tau_nn,stressMax,stressMin,cos(2*theta)); */
-
-/* printf("ratios %lf %lf %lf %lf %lf %lf %lf\n",(stressMax-stressMin)/2,pressure,sigma_nn,sin(2*theta)*(stressMax-stressMin)/2, */
-/* cos(2*theta)*(stressMax-stressMin)/2, */
-/* stressMax,stressMin); */
-
- return sigma_nn;
-}
-
void _MohrCoulomb_StoreCurrentParameters(
void* rheology,
ConstitutiveMatrix* constitutiveMatrix,
@@ -506,16 +410,4 @@
SymmetricTensor_CalcAllEigenvectors( self->currentStress, dim, self->currentEigenvectorList );
SymmetricTensor_CalcAllEigenvectors( strainRate, dim, evectors);
- if(dim == 2)
- {
- self->strainRateSecondInvariant=0.5*fabs(evectors[0].eigenvalue-evectors[1].eigenvalue);
- }
- else
- {
- e0=evectors[0].eigenvalue;
- e1=evectors[1].eigenvalue;
- e2=evectors[2].eigenvalue;
- self->strainRateSecondInvariant=
- sqrt(((e0-e1)*(e0-e1) + (e1-e2)*(e1-e2) + (e2-e0)*(e2-e0))/6.0);
- }
}
More information about the cig-commits
mailing list