ViscousHeating

PURPOSE ^

VISCOUSHEATING - compute the viscous heating on each pentaelem

SYNOPSIS ^

function viscousheating=ViscousHeating(pentaelem,grids,materials,inputs)

DESCRIPTION ^

VISCOUSHEATING - compute the viscous heating on each pentaelem

   Usage:
      viscousheating=ViscousHeating(pentaelem,grids,materials,inputs)

   See also GETSTRAINRATE, GETVISCOSITY3D

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function viscousheating=ViscousHeating(pentaelem,grids,materials,inputs)
0002 %VISCOUSHEATING - compute the viscous heating on each pentaelem
0003 %
0004 %   Usage:
0005 %      viscousheating=ViscousHeating(pentaelem,grids,materials,inputs)
0006 %
0007 %   See also GETSTRAINRATE, GETVISCOSITY3D
0008 
0009     %initialize
0010     viscousheating=0;
0011     volume=0;
0012 
0013     %some variables
0014     NDOF1=1;
0015     numgrids=6;
0016 
0017     %recover material parameters
0018     matice=materials(pentaelem.matid).material;
0019     B=matice.B;
0020 
0021     %Get all element grid data:
0022     xyz_list=getgriddata(pentaelem,grids);
0023 
0024     %recover extra inputs
0025     [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0026     [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0027 
0028     %we need velocities to compute thermal profiles (even if it is a zero
0029     %vector).
0030     if ~velocity_is_present,
0031         error('ViscousHeating error message: input velocity not present!');
0032     end
0033 
0034     %initialize vxvyvz_list
0035     vxvyvz_list=zeros(numgrids,3);
0036 
0037     %Build row indices for elementary vector.
0038     for i=1:numgrids,
0039         doflist=grids(pentaelem.g(i)).grid.doflist;
0040         for j=1:3,
0041             dof=doflist(j);
0042             vxvyvz_list(i,j)=velocity_param(dof);
0043         end    
0044         dof=doflist(1);
0045         if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
0046     end
0047     
0048     %Update material parameter that deals with ice rigidity:
0049     if flow_law_is_present,
0050         B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
0051         matice.B=B_param; clear B_param.
0052     end
0053 
0054     % Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
0055     %get tria gaussian points as well as segment gaussian points. For tria gaussian
0056     %points, order of integration is 2, because we need to integrate the product tB*D*B'
0057     %which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
0058     %points, same deal, which yields 3 gaussian points.
0059     area_order=2;
0060     num_vert_gauss=3;
0061     [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);
0062 
0063     %Start  looping on the number of gaussian points:
0064     for igarea=1:num_area_gauss,
0065         for igvert=1:num_vert_gauss,
0066             %Pick up the gaussian point and its weight:
0067             gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
0068             gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
0069             
0070             %Build Deformational heating
0071             epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0072             viscosity=GetViscosity3d(matice,epsilon); 
0073 
0074             %Build the strain rate matrix
0075             epsilon_matrix=[epsilon(1) epsilon(4) epsilon(5)
0076                     epsilon(4) epsilon(2) epsilon(6)
0077                     epsilon(5) epsilon(6) epsilon(3)];
0078 
0079             epsilon_e=effective_value(epsilon_matrix);
0080 
0081             phi=2*viscosity*epsilon_e^2;
0082 
0083             %Get Jacobian determinant:
0084             Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
0085 
0086             viscousheating=viscousheating+phi*Jdet*gauss_weight;
0087             volume=volume+Jdet*gauss_weight;
0088 
0089         end %for ig=1:vert_gauss,
0090     end %for ig=1:num_area_gauss,
0091 
0092 %Divide viscousheating, integrated over volume, by volume.
0093 viscousheating=viscousheating/volume;
0094 
0095 
0096 end %end function

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003