[cig-commits] r14446 - in long/3D/Gale/trunk: . src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Wed Mar 25 06:39:25 PDT 2009
Author: walter
Date: 2009-03-25 06:39:25 -0700 (Wed, 25 Mar 2009)
New Revision: 14446
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
Log:
r2587 at dante: boo | 2009-03-25 06:39:00 -0700
Clean up Mohr Coulomb and make it handle 3D
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2586
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2587
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2009-03-25 13:39:18 UTC (rev 14445)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2009-03-25 13:39:25 UTC (rev 14446)
@@ -270,10 +270,15 @@
return;
}
- /* Calculate and store values of stress, pressure, velocity gradients, eigenvectors so they are only calculated once */
- _MohrCoulomb_StoreCurrentParameters( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint, xi );
+ /* Calculate and store values of stress, pressure, velocity
+ gradients, inavriants so they are only calculated once */
+ _MohrCoulomb_StoreCurrentParameters( self, constitutiveMatrix,
+ materialPointsSwarm, lElement_I,
+ materialPoint, xi );
- _YieldRheology_ModifyConstitutiveMatrix( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint, xi );
+ _YieldRheology_ModifyConstitutiveMatrix(self, constitutiveMatrix,
+ materialPointsSwarm,lElement_I,
+ materialPoint, xi );
}
double _MohrCoulomb_GetYieldIndicator(
@@ -285,25 +290,9 @@
Coord xi )
{
MohrCoulomb* self = (MohrCoulomb*) rheology;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double sigma_ns;
- Dimension_Index dim = constitutiveMatrix->dim;
- double stressMin;
- double stressMax;
- /* 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);
- sigma_ns = 0.5 * ( stressMax - stressMin );
-
- return fabs(sigma_ns);
+ return SymmetricTensor_2ndInvariant( self->stress,
+ constitutiveMatrix->dim );
}
double _MohrCoulomb_GetYieldCriterion(void* rheology,
@@ -329,8 +318,7 @@
unsigned inds[3];
Grid* elGrid;
Bool inside_boundary;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
+ FeVariable_InterpolateWithinElement(pressureField,lElement_I,xi,&pressure );
cell_I=CellLayout_MapElementIdToCellId(materialPointsSwarm->cellLayout,
lElement_I );
@@ -340,18 +328,12 @@
pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,coord);
}
- /* Add the average of the trace, since with linear pressure elements
- we are not guaranteed to have a stress field with zero stress */
+ /* Normally we add the average of the trace of the stress. With
+ compressible material, you have to do it. But with stabilized
+ linear pressure elements, the non-zero trace is a numerical
+ artifact. So we do not add it. */
- if(dim==2)
- {
- pressure+=(eigenvectorList[0].eigenvalue+eigenvectorList[1].eigenvalue)/2;
- }
- else
- {
- pressure+=(eigenvectorList[0].eigenvalue+eigenvectorList[1].eigenvalue
- +eigenvectorList[2].eigenvalue)/2;
- }
+/* pressure+=self->trace/dim; */
/* Calculate frictional strength. We modify the friction and
cohesion because we have grouped terms from the normal
@@ -379,25 +361,34 @@
effectiveCohesion = _MohrCoulomb_EffectiveCohesion(self,materialPoint,
inside_boundary);
- /* effectiveFrictionCoefficient=tan(phi).
- If factor=sin(atan(1/tan(phi)))
- => factor=cos(alpha)=1/sqrt(1+tan(phi)**2) */
- factor=1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
if(dim==2)
{
+ /* effectiveFrictionCoefficient=tan(phi). If
+ factor=sin(atan(1/tan(phi))) =>
+ factor=cos(phi)=1/sqrt(1+tan(phi)**2) */
+ factor=1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
frictionalStrength = effectiveFrictionCoefficient*pressure*factor
+ effectiveCohesion*factor;
}
else
{
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double stressMin = eigenvectorList[0].eigenvalue;
- double stressMax = eigenvectorList[2].eigenvalue;
+ double cos_phi, sin_phi;
+ /* cos(phi)=1/sqrt(1+tan(phi)**2) */
+ cos_phi=
+ 1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
+ sin_phi=effectiveFrictionCoefficient*cos_phi;
+ factor=2*cos_phi/(sqrt(3.0)*(3-sin_phi));
- frictionalStrength =
- effectiveFrictionCoefficient*factor*(pressure+(stressMin+stressMax)/2)
- + effectiveCohesion*factor;
+ /* The full expression is
+
+ sqrt(J2)=p*2*sin(phi)/(sqrt(3)*(3-sin(phi)))
+ + C*6*cos(phi)/(sqrt(3)*(3-sin(phi)))
+
+ Note the extra factor of 3 for cohesion */
+
+ frictionalStrength = effectiveFrictionCoefficient*factor*pressure
+ + effectiveCohesion*3*factor;
}
/* If the minimumYieldStress is not set, then use the
@@ -488,42 +479,31 @@
MaterialPoint* materialPoint,
Coord xi )
{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- SymmetricTensor strainRate;
- Dimension_Index dim = constitutiveMatrix->dim;
- Eigenvector evectors[3];
- double e0, e1, e2;
- double trace;
- int i;
-
- FeVariable_InterpolateWithinElement( self->pressureField, lElement_I, xi, &self->currentPressure );
- FeVariable_InterpolateWithinElement( self->velocityGradientsField, lElement_I, xi, self->currentVelocityGradient );
-
- TensorArray_GetSymmetricPart( self->currentVelocityGradient, dim, strainRate );
+ MohrCoulomb* self = (MohrCoulomb*) rheology;
+ SymmetricTensor strainRate;
+ Dimension_Index dim = constitutiveMatrix->dim;
+ double strainRateTrace;
+ int i;
+
+ FeVariable_InterpolateWithinElement( self->pressureField, lElement_I, xi,
+ &self->currentPressure );
+ FeVariable_InterpolateWithinElement(self->velocityGradientsField,lElement_I,
+ xi, self->currentVelocityGradient );
+ TensorArray_GetSymmetricPart(self->currentVelocityGradient,dim,strainRate);
- SymmetricTensor_GetTrace(strainRate, dim, &trace);
+ ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate,
+ self->stress );
+ SymmetricTensor_GetTrace(self->stress, dim, &(self->trace));
+ SymmetricTensor_GetTrace(strainRate, dim, &strainRateTrace);
+
+ /* Subtract the trace so that the 2nd invariant is not polluted by
+ the trace */
- /* Subtract the trace. We can use TensorMapST3D even for 2D,
- because it is the same for the xx and yy components */
- for(i=0;i<dim;++i)
- strainRate[TensorMapST3D[i][i]]-=trace/dim;
+ for(i=0;i<dim;++i)
+ {
+ strainRate[TensorMapST3D[i][i]]-=strainRateTrace/dim;
+ self->stress[TensorMapST3D[i][i]]-=self->trace/dim;
+ }
- ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, self->currentStress );
-
- 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);
- }
+ self->strainRateSecondInvariant=SymmetricTensor_2ndInvariant(strainRate,dim);
}
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2009-03-25 13:39:18 UTC (rev 14445)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2009-03-25 13:39:25 UTC (rev 14446)
@@ -78,9 +78,9 @@
Bool boundaryBack; \
double minimumYieldStress; \
/* Stored values that are calculated once for each particle */ \
- Eigenvector currentEigenvectorList[3]; \
+ double trace; \
TensorArray currentVelocityGradient; \
- SymmetricTensor currentStress; \
+ SymmetricTensor stress; \
double currentPressure; \
double strainRateSecondInvariant; \
FeVariable* strainRateField; \
More information about the CIG-COMMITS
mailing list