[cig-commits] r14750 - in long/3D/Gale/trunk: . src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Sat Apr 18 05:05:54 PDT 2009
Author: walter
Date: 2009-04-18 05:05:53 -0700 (Sat, 18 Apr 2009)
New Revision: 14750
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.h
Log:
r2662 at dante: boo | 2009-04-18 05:05:17 -0700
Make NonNewtonian model the power law directly, rather than relying upon Arrhenius
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2661
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2662
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.c 2009-04-18 12:05:43 UTC (rev 14749)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.c 2009-04-18 12:05:53 UTC (rev 14750)
@@ -94,10 +94,21 @@
return self;
}
-void _NonNewtonian_Init( NonNewtonian* self, FeVariable* strainRateInvField, double stressExponent, double refStrainRate, double maxViscosity, double minViscosity) {
+void _NonNewtonian_Init( NonNewtonian* self, FeVariable* strainRateInvField,
+ FeVariable* temperatureField,
+ double n,
+ double T_0,
+ double A,
+ double refStrainRate,
+ double maxViscosity,
+ double minViscosity) {
+
self->strainRateInvField = strainRateInvField;
+ self->temperatureField = temperatureField;
- self->stressExponent = stressExponent;
+ self->n = n;
+ self->T_0 = T_0;
+ self->A = A;
self->refStrainRate = refStrainRate;
self->minViscosity = minViscosity;
self->maxViscosity = maxViscosity;
@@ -124,7 +135,7 @@
void _NonNewtonian_Construct( void* rheology, Stg_ComponentFactory* cf, void* data ){
NonNewtonian* self = (NonNewtonian*)rheology;
- FeVariable* strainRateInvField;
+ FeVariable *strainRateInvField, *temperatureField;
double refStrainRate;
Stream* errorStream = Journal_Register( Error_Type,
NonNewtonian_Type );
@@ -133,25 +144,34 @@
_Rheology_Construct( self, cf, data );
/* TODO: 'Keyfallback' soon to be deprecated/updated */
- strainRateInvField = Stg_ComponentFactory_ConstructByNameWithKeyFallback(
+ strainRateInvField = Stg_ComponentFactory_ConstructByKey(
cf,
self->name,
"StrainRateInvariantField",
- "StrainRateInvariantField",
FeVariable,
True,
data );
+ temperatureField = Stg_ComponentFactory_ConstructByKey(
+ cf,
+ self->name,
+ "TemperatureField",
+ FeVariable,
+ True,
+ data );
refStrainRate=Stg_ComponentFactory_GetDouble( cf, self->name,
"refStrainRate", 0.0 );
if(refStrainRate==0.0)
Journal_Firewall(0,errorStream,
- "refStrainRate must be set to a non-zero number in a NonNewtonian Rheology");
+ "refStrainRate must be set to a non-zero number in a NonNewtonian rheology");
_NonNewtonian_Init(
self,
strainRateInvField,
- Stg_ComponentFactory_GetDouble( cf, self->name, "stressExponent", 1.0 ),
+ temperatureField,
+ Stg_ComponentFactory_GetDouble( cf, self->name, "n", 1.0 ),
+ Stg_ComponentFactory_GetDouble( cf, self->name, "T_0", 1.0 ),
+ Stg_ComponentFactory_GetDouble( cf, self->name, "A", 1.0 ),
refStrainRate,
Stg_ComponentFactory_GetDouble( cf, self->name, "minViscosity", 0.0 ) ,
Stg_ComponentFactory_GetDouble( cf, self->name, "maxViscosity", 0.0 ) );
@@ -170,28 +190,39 @@
double minVis;
double maxVis;
double viscosity;
- double n;
+ double n, T_0, A, T;
+ Stream* errorStream = Journal_Register( Error_Type,
+ NonNewtonian_Type );
- n = self->stressExponent;
+ n = self->n;
+ T_0=self->T_0;
+ A=self->A;
minVis = self->minViscosity;
maxVis = self->maxViscosity;
- if ( fabs( n - 1.0 ) < 1.0e-10 )
- return;
- if ( !constitutiveMatrix->previousSolutionExists )
- /* on the first ever solve, use a ref strainrate */
+
+ FeVariable_InterpolateWithinElement(self->temperatureField,
+ lElement_I,xi,&T);
+ Journal_Firewall(n!=0 && T!=0 && A!=0,errorStream,
+ "Error in NonNewtonian: T, n, and A must all be non-zero:\n\tT=%g\n\tn=%g\n\tA=%g\n",T,n,A);
+
+ viscosity=exp(T_0/(2*n*T))/(2*A);
+
+ if ( fabs( n - 1.0 ) > 1.0e-10 )
{
- strainRateInv = self->refStrainRate;
+ if ( !constitutiveMatrix->previousSolutionExists )
+ /* on the first ever solve, use a ref strainrate */
+ {
+ strainRateInv = self->refStrainRate;
+ }
+ else
+ {
+ /* Calculate Parameters */
+ FeVariable_InterpolateWithinElement(self->strainRateInvField,
+ lElement_I,xi,&strainRateInv);
+ }
+ /* Calculate New Viscosity */
+ viscosity*= pow(strainRateInv,1.0/n - 1.0);
}
- else
- {
- /* Calculate Parameters */
- FeVariable_InterpolateWithinElement(self->strainRateInvField,
- lElement_I,xi,&strainRateInv);
- }
- /* Calculate New Viscosity */
- viscosity =ConstitutiveMatrix_GetIsotropicViscosity(constitutiveMatrix);
- viscosity = pow(2.0*strainRateInv,1.0/n - 1.0) * pow(viscosity,1.0/n);
-
if(maxVis!=0 && viscosity>maxVis)
{
viscosity=maxVis;
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.h 2009-04-18 12:05:43 UTC (rev 14749)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/NonNewtonian.h 2009-04-18 12:05:53 UTC (rev 14750)
@@ -57,7 +57,10 @@
/* Virtual functions go here */ \
/* Material Parameters */\
FeVariable* strainRateInvField; \
- double stressExponent; \
+ FeVariable* temperatureField; \
+ double n; \
+ double T_0; \
+ double A; \
double refStrainRate; \
double minViscosity; \
double maxViscosity; \
More information about the CIG-COMMITS
mailing list