1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine con2phys(vold,vcon,nk,ntmat_,shcon,nshcon,rhcon,nrhcon, 20 & physcon,ithermal,compressible,iturbulent,inomat,mi, 21 & ierr,ifreesurface,dgravity,depth) 22! 23! calculates the physical variable from the conservative 24! variables: 25! for gases: temperature,velocity and pressure 26! for thermal liquids: temperature, velocity and density 27! for athermal liquids: velocity and density 28! 29 implicit none 30! 31 integer compressible,ifreesurface,nrhcon(*),ntmat_,iturbulent, 32 & mi(*),nshcon(*),nk,ithermal(*),i,j,k,imat,ierr,inomat(*) 33! 34 real*8 vold(0:mi(2),*),vcon(nk,0:mi(2)),rhcon(0:1,ntmat_,*),rho, 35 & c1,cp,r,temp,temp0,c2,shcon(0:3,ntmat_,*),physcon(*), 36 & dgravity,depth(*) 37! 38 if(ithermal(1).gt.1) then 39! 40 do i=1,nk 41 imat=inomat(i) 42 temp=vold(0,i) 43! 44 if(compressible.eq.1) then 45! 46! gas : rho*epsilon_t, rho and rho*v known 47! 48! calculation of the static temperature 49! 50 rho=vcon(i,4) 51 r=shcon(3,1,imat) 52 c1=(vcon(i,0)-(vcon(i,1)**2+vcon(i,2)**2+ 53 & vcon(i,3)**2)/(2.d0*rho))/rho 54! 55! the energy density per volume unit cannot be 56! negative 57! 58 if(c1.lt.1.d-10) c1=1.d-10 59! 60 temp0=temp 61 j=0 62 do 63 call materialdata_cp_sec(imat,ntmat_,temp,shcon, 64 & nshcon,cp,physcon) 65 temp=c1/(cp-r)+physcon(1) 66 j=j+1 67 if((dabs(temp-temp0).lt.1.d-4*dabs(temp)).or. 68 & (dabs(temp-temp0).lt.1.d-10)) then 69 vold(0,i)=temp 70 exit 71 endif 72 if(j.gt.100) then 73 write(*,*) 74 & '*ERROR in con2phys: too many iterations' 75 write(*,*) ' for node',i 76 write(*,*) ' increment is recalculated' 77 write(*,*) ' with a higher shock smoothing' 78 ierr=1 79 return 80 endif 81 temp0=temp 82 enddo 83! 84 if(ifreesurface.eq.0) then 85! 86! determining the pressure (gas equation) 87! 88 vold(4,i)=rho*r*(temp-physcon(1)) 89 else 90! 91! determining the pressure (shallow water) 92! 93 vold(4,i)=dgravity*(vcon(i,4)**2-depth(i)**2)/2.d0 94 endif 95! 96! calculating the velocity from the conservative variables 97! 98 do k=1,3 99 vold(k,i)=vcon(i,k)/rho 100 enddo 101! 102! calculating the turbulent quantities from the conservative 103! variables 104! 105 if(iturbulent.ne.0) then 106 do k=5,6 107 vold(k,i)=vcon(i,k)/rho 108 enddo 109 endif 110! 111 cycle 112 else 113! 114! thermal liquid 115! 116 c1=vcon(i,0) 117 c2=(vcon(i,1)**2+vcon(i,2)**2+vcon(i,3)**2)/2.d0 118 temp0=temp 119 j=0 120! 121! iterating to find the temperature 122! 123 do 124 call materialdata_cp_sec(imat,ntmat_,temp,shcon, 125 & nshcon,cp,physcon) 126 call materialdata_rho(rhcon,nrhcon,imat,rho, 127 & temp,ntmat_,ithermal) 128 temp=(c1-c2/rho)/(rho*cp)+physcon(1) 129 j=j+1 130 if((dabs(temp-temp0).lt.1.d-4*dabs(temp)).or. 131 & (dabs(temp-temp0).lt.1.d-10)) then 132 vold(0,i)=temp 133 exit 134 endif 135 if(j.gt.100) then 136 write(*,*) 137 & '*ERROR in con2phys: too many iterations' 138 write(*,*) ' for node',i 139 write(*,*) ' actual temperature ',temp,' K' 140 stop 141 endif 142 temp0=temp 143 enddo 144! 145! storing the density 146! 147 vcon(i,4)=rho 148! 149! calculating the velocity 150! 151 do k=1,3 152 vold(k,i)=vcon(i,k)/rho 153 enddo 154! 155! calculating the turbulent quantities 156! 157 if(iturbulent.ne.0) then 158 do k=5,6 159 vold(k,i)=vcon(i,k)/rho 160 enddo 161 endif 162! 163 endif 164 enddo 165 else 166! 167! athermal calculations 168! 169 do i=1,nk 170! 171! calculating the fictitious pressure for shallow water 172! calculations 173! 174 if(ifreesurface.eq.1) then 175 vold(4,i)=dgravity*(vcon(i,4)**2-depth(i)**2)/2.d0 176 endif 177! 178! calculating the velocity 179! 180 rho=vcon(i,4) 181 do k=1,3 182 vold(k,i)=vcon(i,k)/rho 183 enddo 184! 185! calculating the turbulent quantities 186! 187 if(iturbulent.ne.0) then 188 do k=5,6 189 vold(k,i)=vcon(i,k)/rho 190 enddo 191 endif 192! 193 enddo 194 endif 195! 196 return 197 end 198