0001 function viscousheating=ViscousHeating(triaelem,grids,materials,inputs)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 viscousheating=0;
0012 volume=0;
0013
0014
0015 NDOF1=1;
0016 numgrids=3;
0017
0018
0019 matice=materials(triaelem.matid).material;
0020 B=matice.B;
0021
0022
0023 xyz_list=getgriddata(triaelem,grids);
0024
0025
0026 [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
0027 [flow_law_param flow_law_is_present]=recover_input(inputs,'B');
0028
0029
0030
0031 if ~velocity_is_present,
0032 error('ViscousHeating error message: input velocity not present!');
0033 end
0034
0035
0036 vxvy_list=zeros(numgrids,2);
0037
0038
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
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
0056
0057
0058
0059
0060 [num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights]=GaussTria(2);
0061
0062
0063 for igarea=1:num_area_gauss,
0064
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
0069 epsilon=GetStrainRate(triaelem,vxvy_list,xyz_list,gauss_coord);
0070 viscosity=GetViscosity2d(matice,epsilon);
0071
0072
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
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
0088
0089
0090 viscousheating=viscousheating/volume;
0091
0092
0093 end