


GETVISCOSITY3D - returns viscosity from a tensor, using Glen's flow law
From a string tensor and a material object, return viscosity, using Glen's flow law.
2*B
mu= -------------------------------------------------------------------
2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
where mu is the viscotiy, B the flow law parameter , (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.
This routine is appropriate for 3d elements.
Usage:
mu=GetViscosity3d(matice,epsilon)

0001 function mu=GetViscosity3d(matice,epsilon) 0002 %GETVISCOSITY3D - returns viscosity from a tensor, using Glen's flow law 0003 % 0004 % From a string tensor and a material object, return viscosity, using Glen's flow law. 0005 % 2*B 0006 % mu= ------------------------------------------------------------------- 0007 % 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] 0008 % 0009 % where mu is the viscotiy, B the flow law parameter , (u,v) the velocity 0010 % vector, and n the flow law exponent. 0011 % 0012 % If epsilon is NULL, it means this is the first time Emg is being run, and we 0013 % return 10^14, initial viscosity. 0014 % 0015 % This routine is appropriate for 3d elements. 0016 % 0017 % Usage: 0018 % mu=GetViscosity3d(matice,epsilon) 0019 0020 eps0=10^-27; 0021 if matice.n==1, 0022 %Viscous behaviour. mu=B; 0023 mu=matice.B; 0024 else 0025 %Non linear creep 0026 if norm(epsilon,2)~=0 0027 0028 if size(epsilon,1)==5, 0029 B=matice.B; 0030 n=matice.n; 0031 0032 exx=epsilon(1); 0033 eyy=epsilon(2); 0034 exy=epsilon(3); 0035 exz=epsilon(4); 0036 eyz=epsilon(5); 0037 0038 %Build viscosity: mu=2*B/(2*A^e): 0039 A=exx^2+eyy^2+exy^2+exx*eyy+exz^2+eyz^2+eps0^2; 0040 e=(n-1)/2/n; 0041 mu=2*B/(2*A^e); 0042 0043 elseif size(epsilon,1)==6, 0044 B=matice.B; 0045 n=matice.n; 0046 0047 exx=epsilon(1); 0048 eyy=epsilon(2); 0049 ezz=epsilon(3); 0050 exy=epsilon(4); 0051 exz=epsilon(5); 0052 eyz=epsilon(6); 0053 0054 %Build viscosity: mu=2*B/(2*A^e): 0055 A=exx^2+eyy^2+exy^2+exx*eyy+exz^2+eyz^2+eps0^2; 0056 e=(n-1)/2/n; 0057 mu=2*B/(2*A^e); 0058 0059 else 0060 error('GetViscosity3d error message : wrong size of epsilon'); 0061 end 0062 0063 else 0064 %If epsilon==0, mu=minimum viscosity 0065 mu=10^14; 0066 end 0067 end 0068