Changeset 27164

Timestamp:

07/28/22 06:53:10 (3 years ago)

Author:

vverjans

Message:

CHG frictioncoulomb uses effective presure instead of thickness above floatation

Location:

issm/trunk-jpl/src

Files:

• c/classes/Loads/Friction.cpp (modified) (3 diffs)
• m/classes/frictioncoulomb.m (modified) (1 diff)
• m/classes/frictioncoulomb.py (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -321 +321 @@
 
         /*This routine calculates the basal friction coefficient
-          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
+          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p
+          alpha2= min(drag^2 * Neff ^r * | vel | ^(s-1), drag_coulomb^2 * Neff*/
 
         /*diverse: */
         IssmDouble  r,s;
         IssmDouble  drag_p, drag_q;
-        IssmDouble  thickness,base,floatation_thickness,sealevel;
         IssmDouble  drag_coefficient,drag_coefficient_coulomb;
         IssmDouble  alpha2,alpha2_coulomb;
@@ -333 +333 @@
         element->GetInputValue(&drag_p,gauss,FrictionPEnum);
         element->GetInputValue(&drag_q,gauss,FrictionQEnum);
-        element->GetInputValue(&thickness, gauss,ThicknessEnum);
-        element->GetInputValue(&base, gauss,BaseEnum);
-        element->GetInputValue(&sealevel, gauss,SealevelEnum);
         element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
         element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
@@ -360 +357 @@
         }
 
-        floatation_thickness=0;
-        if(base<0) floatation_thickness=-(rho_water/rho_ice)*base;
         if(vmag==0.){
                 alpha2_coulomb=0.;
         }
         else{
-                alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
+                //alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
+                alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*Neff/vmag;
         }
 

issm/trunk-jpl/src/m/classes/frictioncoulomb.m

@@ -64 +64 @@
                 end % }}}
                 function disp(self) % {{{
-                        disp(sprintf('Basal shear stress parameters: Sigma_b = min( coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b\n, coefficientcoulomb^2 * rho_i * g * (h-h_f)) (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p, floatation thickness h_f=max(0,-rho_sw / rho_i * bed))'));
+                        disp(sprintf('Basal shear stress parameters: Sigma_b = min( coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b\n, coefficientcoulomb^2 * Neff) (effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p, floatation thickness h_f=max(0,-rho_sw / rho_i * bed))'));
                         fielddisplay(self,'coefficient','power law (Weertman) friction coefficient [SI]');
                         fielddisplay(self,'coefficientcoulomb','Coulomb friction coefficient [SI]');

issm/trunk-jpl/src/m/classes/frictioncoulomb.py

@@ -28 +28 @@
     def __repr__(self):  # {{{
         s = 'Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n'
-        s += 'coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n'
+        s += 'coefficientcoulomb^2 * Neff), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n'
         s += '{}\n'.format(fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
         s += '{}\n'.format(fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))