[cig-commits] r11843 - in long/3D/Gale/trunk: . input/benchmarks src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Tue Apr 22 11:23:00 PDT 2008
Author: walter
Date: 2008-04-22 11:23:00 -0700 (Tue, 22 Apr 2008)
New Revision: 11843
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/input/benchmarks/extension.xml
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
Log:
r2103 at earth: boo | 2008-04-22 10:22:24 -0700
Rough implementation of corrected MohrCoulomb. Not working in 3D
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2102
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2103
Modified: long/3D/Gale/trunk/input/benchmarks/extension.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/extension.xml 2008-04-22 18:22:57 UTC (rev 11842)
+++ long/3D/Gale/trunk/input/benchmarks/extension.xml 2008-04-22 18:23:00 UTC (rev 11843)
@@ -331,7 +331,7 @@
</struct>
<struct name="crustViscosity">
<param name="Type">MaterialViscosity</param>
- <param name="eta0">5.0e7</param>
+ <param name="eta0">1.0e12</param>
</struct>
<struct name="strainWeakening">
<param name="Type">StrainWeakening</param>
@@ -345,9 +345,11 @@
<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>
@@ -422,25 +424,25 @@
<struct name="conditionFunctions">
<param name="Type">StgFEM_StandardConditionFunctions</param>
</struct>
- <struct name="kineticFriction">
- <param name="Type">StressBC</param>
- <param name="ForceVector">mom_force</param>
- <param name="Swarm">picIntegrationPoints</param>
- <param name="wall">bottom</param>
- <param name="StressField">StressField</param>
- <param name="PressureField">PressureField</param>
- <param name="VelocityField">VelocityField</param>
- <param name="friction_type">double</param>
- <param name="friction_value">0.30</param>
- <param name="vx_type">func</param>
- <param name="vx_value">StepFunction</param>
- </struct>
+<!-- <struct name="kineticFriction"> -->
+<!-- <param name="Type">StressBC</param> -->
+<!-- <param name="ForceVector">mom_force</param> -->
+<!-- <param name="Swarm">picIntegrationPoints</param> -->
+<!-- <param name="wall">bottom</param> -->
+<!-- <param name="StressField">StressField</param> -->
+<!-- <param name="PressureField">PressureField</param> -->
+<!-- <param name="VelocityField">VelocityField</param> -->
+<!-- <param name="friction_type">double</param> -->
+<!-- <param name="friction_value">0.30</param> -->
+<!-- <param name="vx_type">func</param> -->
+<!-- <param name="vx_value">StepFunction</param> -->
+<!-- </struct> -->
</struct>
<list name="plugins">
<param>Underworld_EulerDeform</param>
- <param>Underworld_DumpSwarm</param>
+<!-- <param>Underworld_DumpSwarm</param> -->
<param>Underworld_VTKOutput</param>
- <param>Gale_SurfaceProcess</param>
+<!-- <param>Gale_SurfaceProcess</param> -->
</list>
<param name="maxTimeSteps">500</param>
<param name="dumpEvery">1</param>
@@ -554,17 +556,6 @@
<param name="type">double</param>
<param name="value">0</param>
</struct>
- </list>
- </struct>
-
- <struct>
- <param name="type">StaticFrictionVC</param>
- <param name="wall">bottom</param>
- <param name="StaticFriction">0.33</param>
- <param name="includeUpperX">False</param>
- <param name="StressField">StressField</param>
- <param name="PressureField">PressureField</param>
- <list name="variables">
<struct>
<param name="name">vx</param>
<param name="type">func</param>
@@ -572,6 +563,22 @@
</struct>
</list>
</struct>
+
+<!-- <struct> -->
+<!-- <param name="type">StaticFrictionVC</param> -->
+<!-- <param name="wall">bottom</param> -->
+<!-- <param name="StaticFriction">0.33</param> -->
+<!-- <param name="includeUpperX">False</param> -->
+<!-- <param name="StressField">StressField</param> -->
+<!-- <param name="PressureField">PressureField</param> -->
+<!-- <list name="variables"> -->
+<!-- <struct> -->
+<!-- <param name="name">vx</param> -->
+<!-- <param name="type">func</param> -->
+<!-- <param name="value">StepFunction</param> -->
+<!-- </struct> -->
+<!-- </list> -->
+<!-- </struct> -->
</list>
</struct>
@@ -596,12 +603,12 @@
<param name="StepFunctionValue">6.94444444444e-6</param>
<param name="StepFunctionDim">0</param>
<param name="StepFunctionLessThan">False</param>
- <param name="checkpointEvery">1</param>
+<!-- <param name="checkpointEvery">1</param> -->
<param name="mgLevels">2</param>
<param name="dtFactor">1.0</param>
-<!-- <param name="journal.info">True</param> -->
-<!-- <param name="journal.debug">True</param> -->
-<!-- <param name="journal-level.info">2</param> -->
-<!-- <param name="journal-level.debug">2</param> -->
+ <param name="journal.info">True</param>
+ <param name="journal.debug">True</param>
+ <param name="journal-level.info">2</param>
+ <param name="journal-level.debug">2</param>
</StGermainData>
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:22:57 UTC (rev 11842)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-04-22 18:23:00 UTC (rev 11843)
@@ -272,40 +272,48 @@
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; */
+ 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 *\/ */
+ /* Failure criterion from stress field - we can calculate
+ * whether the material has failed at the current stress from
+ * the principle stresses */
-/* stressMin = eigenvectorList[0].eigenvalue; */
-/* stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue); */
+ theta = 0.5 * atan( 1.0/ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint ) );
+
+ stressMin = eigenvectorList[0].eigenvalue;
+ stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
-/* sigma_ns = 0.5 * ( stressMax - stressMin ); */
+ sigma_ns = 0.5 * ( stressMax - stressMin );
-/* return fabs(sigma_ns); */
- SymmetricTensor stress;
- SymmetricTensor strainRate;
- double stressInv;
+/* printf("theta %lf\n",theta); */
+ return fabs(sigma_ns);
+
+ {
+/* SymmetricTensor stress; */
+/* SymmetricTensor strainRate; */
+/* double stressInv; */
- /* Get Strain Rate */
- FeVariable_InterpolateWithinElement( self->strainRateField, lElement_I, xi, strainRate );
+/* /\* Get Strain Rate *\/ */
+/* FeVariable_InterpolateWithinElement( self->strainRateField, lElement_I, xi, strainRate ); */
- /* Get Stress */
- ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, stress );
- stressInv = SymmetricTensor_2ndInvariant( stress, constitutiveMatrix->dim );
+/* /\* Get Stress *\/ */
+/* ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, stress ); */
+/* stressInv = SymmetricTensor_2ndInvariant( stress, constitutiveMatrix->dim ); */
-/* printf("stress %lf %lf %lf\n",stressInv,xi[0],xi[1]); */
+/* printf("indicator %lf %lf\n",fabs(sigma_ns),stressInv); */
- return stressInv;
+/* return stressInv; */
+ }
+/* return fabs(sigma_ns); */
+
}
double _MohrCoulomb_GetYieldCriterion(
@@ -322,14 +330,32 @@
double effectiveFrictionCoefficient;
double frictionalStrength;
double sigma_nn;
+
+ FeVariable* pressureField = self->pressureField;
+ double pressure;
+ double theta;
+ double factor;
+
+ FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
+
/* Calculate frictional strength */
effectiveFrictionCoefficient = _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint );
effectiveCohesion = _MohrCoulomb_EffectiveCohesion( self, materialPoint );
- sigma_nn = _MohrCoulomb_Sigma_nn( self, materialPoint, constitutiveMatrix->dim );
+/* sigma_nn = _MohrCoulomb_Sigma_nn( self, materialPoint, constitutiveMatrix->dim ); */
- frictionalStrength = effectiveFrictionCoefficient * sigma_nn + effectiveCohesion;
+/* 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 */
minimumYieldStress = self->minimumYieldStress;
if(minimumYieldStress==0.0)
@@ -357,17 +383,41 @@
Dimension_Index dim = constitutiveMatrix->dim;
double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
-/* printf("yielded %lf %lf %lf %lf %lf %lf\n", */
+/* 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, self->strainRateSecondInvariant, */
-/* viscosity * yieldCriterion/yieldIndicator, */
-/* yieldCriterion/(2*self->strainRateSecondInvariant)); */
+/* yieldIndicator, viscosity*beta, */
+/* ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix )); */
/* viscosity = viscosity * yieldCriterion/yieldIndicator; */
- viscosity = yieldCriterion/(2*self->strainRateSecondInvariant);
+/* viscosity = yieldCriterion/(2*self->strainRateSecondInvariant); */
- ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
+/* ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity ); */
}
double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) {
@@ -405,10 +455,20 @@
double stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
double theta = 0.5 * atan( 1.0/ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint ) );
- tau_nn = 0.5 * (( stressMax + stressMin ) + cos( 2.0 * theta ) * ( stressMax - stressMin ));
+ tau_nn =cos( 2.0 * theta ) * ( stressMax - stressMin )*0.5;
+/* tau_nn = 0.5 * (( stressMax + stressMin ) + cos( 2.0 * theta ) * ( stressMax - stressMin )); */
- sigma_nn = pressure - tau_nn;
+/* 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;
}
@@ -425,11 +485,22 @@
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 );
+
+ SymmetricTensor_GetTrace(strainRate, dim, &trace);
+
+ /* Subtract the trace (which should be zero anyway). 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;
+
ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, self->currentStress );
SymmetricTensor_CalcAllEigenvectors( self->currentStress, dim, self->currentEigenvectorList );
More information about the cig-commits
mailing list