


GETVISCOSITY2 - returns viscosity for control method
Return viscosity in 2d, for triaelem, used in the inverse control methods.
Careful, this is not the visocsity used by the diagnostic equations!
2*(1-n)/2n
mu2= -------------------------------------------------------------------
2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(3n-1)/2n]
where mu2 is the viscotiy, (u,v) the velocity
vector, and n the flow law exponent.
If epsilon is NULL, it means this is the first time Emg is being run, and we
return 10^14, initial viscosity.
usage:
mu2=getviscosity2(matice,epsilon)

0001 function mu2=GetViscosity2(matice,epsilon) 0002 %GETVISCOSITY2 - returns viscosity for control method 0003 % 0004 % Return viscosity in 2d, for triaelem, used in the inverse control methods. 0005 % Careful, this is not the visocsity used by the diagnostic equations! 0006 % 2*(1-n)/2n 0007 % mu2= ------------------------------------------------------------------- 0008 % 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(3n-1)/2n] 0009 % 0010 % where mu2 is the viscotiy, (u,v) the velocity 0011 % vector, and n the flow law exponent. 0012 % 0013 % If epsilon is NULL, it means this is the first time Emg is being run, and we 0014 % return 10^14, initial viscosity. 0015 % 0016 % usage: 0017 % mu2=getviscosity2(matice,epsilon) 0018 0019 eps0=10^-27; 0020 if matice.n==1, 0021 error('control method not implemented yet for viscous case'); 0022 else 0023 %Non linear creep 0024 if norm(epsilon,2)~=0, 0025 0026 n=matice.n; 0027 0028 exx=epsilon(1,:); 0029 eyy=epsilon(2,:); 0030 exy=epsilon(3,:); 0031 0032 %Build viscosity: mu2=2*B/(2*A^e): 0033 A=exx.^2+eyy.^2+exy.^2+exx.*eyy+eps0^2; 0034 e=(3*n-1)/2/n; 0035 mu2=(2*(1-n)/2/n*ones(1,size(epsilon,2)))./(2*A.^e); 0036 else 0037 mu2=10^14*size(epsilon,2); 0038 end 0039 end 0040