Changeset 24072

Timestamp:

07/09/19 08:21:00 (6 years ago)

Author:

Mathieu Morlighem

Message:

BUG: improved poorly coded law

File:

• issm/trunk-jpl/src/m/materials/nye.m (modified) (2 diffs)

issm/trunk-jpl/src/m/materials/nye.m ¶

@@ -1 +1 @@
 function rigidity = nye(temperature,ice_type)
-        %NYE - figure out the rigidity of ice (either CO2 or H2O) for a given temperature
-        %   rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law sigma=B*e(1/3) (Nye, p2000). 
-        %   temperature is in Kelvin degrees
-        %   Usage:
-        %      rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
+%NYE - Nye viscosity coefficient
+%
+%   Compute rigidity of ice (either CO2 or H2O) for a given temperature
+%   rigidity (in s^(1/3)Pa) is the flow law parameter in the flow law
+%   sigma=B*e(1/3) (Nye, p2000).  temperature is in Kelvin degrees
+%
+%   Usage:
+%      rigidity=nye(temperature,ice_type) % ice_type = 1: CO2 ice // ice_type = 2: H2O ice
 
         if (any(temperature < 130) || any(temperature > 273))
@@ -12 +15 @@
         Rg = 8.3144598; % J mol^-1 K^-1
 
-        if ice_type == 1 % CO2 ice
+        if(ice_type==1) %CO2 ice
                 A_const = 10^(10.8); % s^-1 MPa
-                Q = 63000; % J mol^-1
-                n = 7; % Glen's exponent
-                T = 250:-5:100; % K
-                A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa
-                B = A.^(-1/n)*1e6; % s^(1/n) Pa
+                Q       = 63000;     % J mol^-1
+                n       = 7;         % Glen's exponent
+
+        elseif(ice_type==2) % H2O ice
+                A_const = 9e4;   % s^-1 MPa
+                Q       = 60000; % J mol^-1
+                n       = 3;     % Glen's exponent
+
+        else
+                error('ice type not supported');
         end
 
-        if ice_type == 2 % H2O ice
-                A_const = 9e4; % s^-1 MPa
-                Q = 60000; %  J mol^-1
-                n = 3; % Glen's exponent
-                T = -100:5:0; % Celsius 
-                T = T + 273; % K
-                A = A_const*exp(-Q./(T*Rg)); % s^-1 MPa
-                B = A.^(-1/n)*1e6; % s^(1/n) Pa
-        end
-
-        % Now, do a cubic fit between Temp and B: 
-        fittedmodel = fit(T',B','cubicspline');
-        rigidity = fittedmodel(temperature);
+        %Arhenius law
+        A = A_const*exp(-Q./(temperature*Rg)); % s^-1 MPa
+        B = A.^(-1/n)*1e6; % s^(1/n) Pa
 
 end