[cig-commits] commit: Bring in Gale's version of NonNewtonian

[Mercurial]

changeset:   690:be936d65eb8c
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Tue Dec 22 16:37:51 2009 -0800
files:       Rheology/src/NonNewtonian.c Rheology/src/NonNewtonian.h
description:
Bring in Gale's version of NonNewtonian


diff -r d058abdcb8f2 -r be936d65eb8c Rheology/src/NonNewtonian.c
--- a/Rheology/src/NonNewtonian.c	Mon Dec 21 13:27:10 2009 -0800
+++ b/Rheology/src/NonNewtonian.c	Tue Dec 22 16:37:51 2009 -0800
@@ -95,10 +95,24 @@ NonNewtonian* _NonNewtonian_New(
 	return self;
 }
 
-void _NonNewtonian_Init( NonNewtonian* self, FeVariable* strainRateInvField, double stressExponent ) {
+void _NonNewtonian_Init( NonNewtonian* self, FeVariable* strainRateInvField,
+                         FeVariable* temperatureField,
+                         double n,
+                         double T_0,
+                         double A,
+                         double refStrainRate,
+                         double minViscosity,
+                         double maxViscosity) {
+
 	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;
 
 	Rheology_SetToNonLinear( self );
 }
@@ -122,7 +136,10 @@ void* _NonNewtonian_DefaultNew( Name nam
 
 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 );
 
 	/* Construct Parent */
 	_Rheology_Construct( self, cf, data );
@@ -136,12 +153,30 @@ void _NonNewtonian_Construct( void* rheo
 		FeVariable, 
 		True,
 		data );
-	/*strainRateInvField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
-				"StrainRateInvariantField", FeVariable, True);*/
+	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");
+
 	_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 ) );
 }
 
 void _NonNewtonian_ModifyConstitutiveMatrix( 
@@ -154,21 +189,49 @@ void _NonNewtonian_ModifyConstitutiveMat
 {
 	NonNewtonian*	              self              = (NonNewtonian*) rheology;
 	double                        strainRateInv;
+	double                        minVis;
+	double                        maxVis;
 	double                        viscosity;
-	double                        n;
+	double                        n, T_0, A, T;
+	Stream* errorStream = Journal_Register( Error_Type,
+                                                NonNewtonian_Type );
 
-	/* Don't want to yield on the first ever solve */
-	if ( !constitutiveMatrix->previousSolutionExists )
-		return;
+	n = self->n;
+        T_0=self->T_0;
+        A=self->A;
+	minVis = self->minViscosity;
+	maxVis = self->maxViscosity;
 
-	/* Calculate Parameters */
-	FeVariable_InterpolateWithinElement( self->strainRateInvField, lElement_I, xi, &strainRateInv );
-	n = self->stressExponent;
-	if ( fabs( n - 1.0 ) < 1.0e-10 )
-		return;
+        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);
 
-	/* Calculate New Viscosity */
-	viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
-	viscosity = pow(2.0 * strainRateInv, 1.0/n - 1.0) * pow(viscosity,1.0/n);
-	ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
+        viscosity=exp(T_0/(n*T))/(2*A);
+
+	if ( fabs( n - 1.0 ) > 1.0e-10 )
+          {
+            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);
+          }
+        if(maxVis!=0 && viscosity>maxVis)
+          {
+            viscosity=maxVis;
+          }
+        if(minVis!=0 && viscosity<minVis)
+          {
+            viscosity=minVis;
+          }
+	ConstitutiveMatrix_SetIsotropicViscosity(constitutiveMatrix,viscosity);
 }
diff -r d058abdcb8f2 -r be936d65eb8c Rheology/src/NonNewtonian.h
--- a/Rheology/src/NonNewtonian.h	Mon Dec 21 13:27:10 2009 -0800
+++ b/Rheology/src/NonNewtonian.h	Tue Dec 22 16:37:51 2009 -0800
@@ -57,7 +57,15 @@
 		/* 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;  \
+
+
 
 	struct NonNewtonian { __NonNewtonian };