0001 function viscousheating=ViscousHeating(pentaelem,grids,materials,inputs)
0002
0003
0004
0005
0006
0007
0008
0009
0010 viscousheating=0;
0011 volume=0;
0012
0013
0014 NDOF1=1;
0015 numgrids=6;
0016
0017
0018 matice=materials(pentaelem.matid).material;
0019 B=matice.B;
0020
0021
0022 xyz_list=getgriddata(pentaelem,grids);
0023
0024
0025 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0026 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0027
0028
0029
0030 if ~velocity_is_present,
0031 error('ViscousHeating error message: input velocity not present!');
0032 end
0033
0034
0035 vxvyvz_list=zeros(numgrids,3);
0036
0037
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
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
0055
0056
0057
0058
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
0064 for igarea=1:num_area_gauss,
0065 for igvert=1:num_vert_gauss,
0066
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
0071 epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
0072 viscosity=GetViscosity3d(matice,epsilon);
0073
0074
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
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
0090 end
0091
0092
0093 viscousheating=viscousheating/volume;
0094
0095
0096 end