[cig-commits] r13670 - in long/3D/Gale/trunk: . src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 12 13:46:35 PST 2008
Author: walter
Date: 2008-12-12 13:46:35 -0800 (Fri, 12 Dec 2008)
New Revision: 13670
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
Log:
r2421 at dante: boo | 2008-12-12 13:44:35 -0800
3D fixes and improved computation of the new viscosity for MohrCoulomb
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2420
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2421
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-12-12 21:46:32 UTC (rev 13669)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-12-12 21:46:35 UTC (rev 13670)
@@ -259,25 +259,17 @@
double stressMin;
double stressMax;
- /* 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. */
+ /* For the yield indicator, we use the difference between the
+ maximum and minimum principal stresses. This is not the
+ actual shear stress. The actual shear stress has another
+ geometric term (sin(2*theta)). Using just the difference
+ simplifies the term for the yield criterion and viscosity
+ correction. */
stressMin = eigenvectorList[0].eigenvalue;
- stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
+ stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue
+ : eigenvectorList[2].eigenvalue);
sigma_ns = 0.5 * ( stressMax - stressMin );
-
- 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));
-
- sigma_ns-=(stressMax+stressMin)*0.5*alpha*factor;
- }
return fabs(sigma_ns);
}
@@ -328,58 +320,73 @@
}
-double _MohrCoulomb_GetYieldCriterion(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi )
+double _MohrCoulomb_GetYieldCriterion(void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi)
{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double minimumYieldStress;
- double effectiveCohesion;
- double effectiveFrictionCoefficient;
- double frictionalStrength;
+ MohrCoulomb* self = (MohrCoulomb*) rheology;
+ double minimumYieldStress;
+ double effectiveCohesion;
+ double effectiveFrictionCoefficient;
+ double frictionalStrength;
+ Dimension_Index dim = constitutiveMatrix->dim;
+
+ FeVariable* pressureField = self->pressureField;
+ double pressure;
+ double theta;
+ double factor;
+ Cell_Index cell_I;
+ Coord coord;
+
+ FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
+
+ cell_I=CellLayout_MapElementIdToCellId(materialPointsSwarm->cellLayout,
+ lElement_I );
+ FeMesh_CoordLocalToGlobal(pressureField->feMesh, cell_I, xi, coord);
+ pressure+=PressureFunction(coord);
+
+ /* 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 );
+
+ theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
+ factor=sin(2*theta);
+
+ if(dim==2)
+ {
+ frictionalStrength = effectiveFrictionCoefficient*pressure*factor
+ + effectiveCohesion*factor;
+ }
+ else
+ {
+ Eigenvector* eigenvectorList = self->currentEigenvectorList;
+ double stressMin = eigenvectorList[0].eigenvalue;
+ double stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue
+ : eigenvectorList[2].eigenvalue);
- FeVariable* pressureField = self->pressureField;
- double pressure;
- double theta;
- double factor;
- Cell_Index cell_I;
- Coord coord;
-
- FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
-
- cell_I=CellLayout_MapElementIdToCellId(materialPointsSwarm->cellLayout,
- lElement_I );
- FeMesh_CoordLocalToGlobal(pressureField->feMesh, cell_I, xi, coord);
- pressure+=PressureFunction(coord);
-
- /* 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 );
-
- theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
-
- factor=sin(2*theta);
-
- frictionalStrength = effectiveFrictionCoefficient * pressure*factor + effectiveCohesion*factor;
-
- /* 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*factor)
- frictionalStrength = minimumYieldStress*factor;
- return frictionalStrength;
+ frictionalStrength =
+ effectiveFrictionCoefficient*factor*(pressure+(stressMin+stressMax)/2)
+ + effectiveCohesion*factor;
+ }
+
+ /* 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*factor)
+ frictionalStrength = minimumYieldStress*factor;
+ return frictionalStrength;
}
/* This scales the viscosity so that the stress is the yield stress.
@@ -397,11 +404,10 @@
double yieldCriterion,
double yieldIndicator )
{
- double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
+ MohrCoulomb* self = (MohrCoulomb*)rheology;
+ double viscosity = yieldCriterion/(2*self->strainRateSecondInvariant);
- double beta=1.0-yieldCriterion/(yieldIndicator);
-
- ConstitutiveMatrix_IsotropicCorrection( constitutiveMatrix, -viscosity * beta );
+ ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
}
double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) {
@@ -462,4 +468,17 @@
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