ViscousHeating

PURPOSE ^

VISCOUSHEATING - compute the viscous heating on each triaelem

SYNOPSIS ^

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

DESCRIPTION ^

VISCOUSHEATING - compute the viscous heating on each triaelem

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

   See also GETSTRAINRATE, GETVISCOSITY3D

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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