1 subroutine getthm (pres, T, Sclr, rk, rho, 2 & ei, h, s, cv, cp, 3 & alfap, betaT, gamb, c ) 4c 5c----------------------------------------------------------------------- 6c 7c This subroutine calculates the thermodynamic properties. 8c 9c The different possibilities are: 10c ipress = 0 : calorifically perfect gas 11c ipress = 1 : thermally perfect gas 12c ipress = 2 : mixture of thermally perfect gases in 13c thermo-chemical equilibrium 14c 15c The options available are: 16c 17c ithm = 2 : given rho and T, compute pres and engBC 18c ithm = 3 : given pres and T 19c ithm = 4 : given pres and T, engBC 20c ithm = 6 : given pres and T, compute rho, ei 21c (also given sclr for levelset) 22c ithm = 7 : given pres and T, compute rho, ei, h, cv, cp, 23c alfap, betaT, gamb, c 24c 25c Variables: 26c 27c pres (npro) : pressure 28c T (npro) : temperature 29c rk (npro) : specific kinetic energy 30c rho (npro) : density 31c ei (npro) : internal energy 32c h (npro) : enthalpy 33c s (npro) : entropy 34c cv (npro) : specific heat at constant volume 35c cp (npro) : specific heat at constant pressure 36c alfap (npro) : expansivity 37c betaT (npro) : isothermal compressibility 38c gamb (npro) : gamma-bar (defined in paper by Chalot et al.) 39c c (npro) : speed of sound 40c 41c 42c Zdenek Johan, Spring 1990. 43c Frederic Chalot, Summer 1990. 44c Zdenek Johan, Winter 1991. (Fortran 90) 45c----------------------------------------------------------------------- 46c 47 include "common.h" 48c 49 dimension pres(npro), Sclr(npro), 50 & T(npro), rk(npro), 51 & rho(npro), ei(npro), 52 & h(npro), s(npro), 53 & cv(npro), cp(npro), 54 & alfap(npro), betaT(npro), 55 & gamb(npro), c(npro), 56 & rsrhol(npro), rsrhog(npro), 57 & tmpg(npro), tmpl(npro) 58c 59 dimension Texp1(npro), Texp2(npro) 60 real*8 prop_blend(npro),test_it(npro) 61 62c ttim(27) = ttim(27) - secs(0.0) 63c 64c.... get the property type flag 65c 66 ipress = matflg(1,1) 67c 68c.... ***********************> IPRESS = 0 <*************************** 69c 70 if (ipress .eq. 0) then 71c 72c.... ---------------------> ithm = 1 or 2 <-------------------------- 73c 74c 75 if (ithm .eq. 2) then 76 pres = Rgas * rho * T 77c 78 79c.... compute engBC (internal energy in this case) 80c 81c engBC = T * (Rgas / gamma1) 82c 83! flops = flops + npro 84c 85 endif 86c 87c.... ---------------------> ithm = 3 or 4 <-------------------------- 88c 89 if ((ithm .eq. 3) .or. (ithm .eq. 4)) then 90c 91 endif 92c 93c if (ithm .eq. 4) then 94c 95c.... compute engBC (enthalpy in this case) 96c 97c engBC = T * (Rgas * gamma / gamma1) 98c 99! flops = flops + npro 100c 101c endif 102c 103c.... --------------------> ithm = 5, 6 or 7 <------------------------ 104c 105c 106 if (ithm .ge. 6) then 107c 108c.... compute density and internal energy 109c 110 if (iLSet .eq. 0)then 111 rho = pres / ( Rgas * T ) 112c 113 else ! two fluid properties used in this model 114 115! Smooth the tranistion of properties for a "distance" of epsilon_ls 116! around the interface. Here "distance" is define as the value of the 117! levelset function. If the levelset function is properly defined, 118! this is the true distance normal from the front. Of course, the 119! distance is in a driection normal to the front. 120 121 do i= 1, npro 122 if (sclr(i) .lt. - epsilon_ls)then 123 prop_blend(i) = zero 124 elseif (abs(sclr(i)) .le. epsilon_ls)then 125 prop_blend(i) = 0.5*(one + Sclr(i)/epsilon_ls + 126 & (sin(pi*Sclr(i)/epsilon_ls))/pi ) 127 elseif (sclr(i) .gt. epsilon_ls) then 128 prop_blend(i) = one 129 endif 130 enddo 131 fact = datmat(1,1,2)/datmat(1,1,1) 132c call eqs(pres,T,rsrhol) 133c rsrhol(:) = pres(:) / ( Rgas * T(:) ) 134c rsrhog(:) = fact* pres(:) / ( Rgas * T(:)) 135 rsrhog(:) = pres(:) / ( Rgas * T(:)) 136 rsrhol(:) = datmat(1,1,1)* 137 & (1+0.000000000517992*(pres-18.02)) 138 rho(:) = rsrhol(:)*prop_blend(:)+rsrhog(:) 139 & *(1-prop_blend(:)) 140c ... for the VOF case..just in case if we want to run VOF 141c$$$ prop_blend(:) = min((max(sclr(:),0.0)),1.0) 142c$$$ rho(:)=rsrhol(:) * prop_blend(:) + rsrhog(:) * (1-prop_blend(:)) 143c 144c 145 endif 146c Calculate Internal Energy 147 if (ilset .eq. 0) then 148 ei = T * (Rgas / gamma1) 149 else 150 tmpg = T * (Rgas / gamma1) !for gas phase 151 tmpl = 3.264*1000*T !cv*T for liquid phase 152c tmpl = T* 8.314* 18.0/ ((3.598/3.264) -1.0) 153 ei(:) = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 154 endif 155! flops = flops + 4*npro 156 endif ! end if(ithm==6) 157c 158 if (ithm .ge. 7) then 159c 160c.... compute enthalpy, cv, cp, alfap, betaT, gamb and c 161c 162 if (ilset .eq. 0) then 163 h = T * (Rgas * gamma / gamma1) 164c 165 cv = Rgas / gamma1 166 cp = Rgas * gamma / gamma1 167c 168 alfap = one / T 169 betaT = one / pres 170c 171 gamb = gamma1 172 c = sqrt( (gamma * Rgas) * T ) 173c 174 else 175c 176 do i= 1, npro 177 if (sclr(i) .lt. - epsilon_ls)then 178 prop_blend(i) = zero 179 elseif (abs(sclr(i)) .le. epsilon_ls)then 180 prop_blend(i) = 0.5*(one + Sclr(i)/epsilon_ls + 181 & (sin(pi*Sclr(i)/epsilon_ls))/pi ) 182 elseif (sclr(i) .gt. epsilon_ls) then 183 prop_blend(i) = one 184 endif 185 enddo 186 rsrhol(:) = datmat(1,1,1)* 187 & (1+0.000000000517992*(pres-18.02)) 188 tmpg = T * (Rgas * gamma / gamma1) ! enthalpy of gas phase 189 tmpl = T*3.264*1000 + pres/rsrhol ! enthalpy of liquid phase 190c tmpl = T*8.314*18.0*3.598/(3.598-3.264) 191 h = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 192c 193 tmpg = Rgas / gamma1 ! cv for gas phase 194 tmpl = 3.264*1000.0 ! cv=cp for liquid phase 195 cv = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 196c 197 tmpg = Rgas * gamma / gamma1 ! cp for gas phase 198 tmpl = 3.264*1000.0 ! cp=cv for liquid phase 199 cp = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 200c 201 tmpg = one / T ! alfap for gas phase 202 tmpl = zero+epsM ! alfap for nearly incompressible liquid 203 alfap = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 204c 205 tmpg = one / pres ! betaT for gas phase 206 tmpl = 0.000000000517992 ! betaT for nearly incompressible liquid 207 betaT = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 208c 209 tmpg = gamma1 ! gamb for gas phase 210 tmpl = 0.0 ! gamb for liquid phase 211 gamb = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 212c 213 tmpg = sqrt( (gamma * Rgas) * T ) ! c for gas phase 214c tmpl = sqrt( (3.598/3.264 * 8.314*18) * T ) ! for liquid 215 tmpl = sqrt(pres/rsrhol) 216 c = tmpl(:)*prop_blend(:)+tmpg*(1-prop_blend(:)) 217c 218 219 endif 220! flops = flops + 12*npro 221c 222 endif 223c 224c 225c.... end of ipress = 0 226c 227 endif 228c 229c.... ***********************> IPRESS = 1 <*************************** 230c 231 if (ipress .eq. 1) then 232c 233c.... ---------------------> ithm = 1 or 2 <-------------------------- 234c 235c 236 if (ithm .eq. 2) then 237c 238c.... compute engBC (internal energy in this case) 239c 240c engBC = Rgas * T / gamma1 241c 242 endif 243c 244c.... ---------------------> ithm = 3 or 4 <-------------------------- 245c 246c 247c if (ithm .eq. 4) then 248c 249c.... compute engBC (enthalpy in this case) 250c 251c engBC = Rgas * T * gamma / gamma1 252c 253c endif 254c 255c.... --------------------> ithm = 5, 6 or 7 <------------------------ 256c 257c 258 if (ithm .ge. 6) then 259c 260c.... compute density and internal energy 261c 262 Texp1 = exp ( - Tvib(1)/T ) 263 Texp2 = exp ( - Tvib(2)/T ) 264c 265 if (iLSet .eq. 0)then 266 rho = pres / ( Rgas * T ) 267c 268 else ! two fluid properties used in this model 269 270! Smooth the tranistion of properties for a "distance" of epsilon_ls 271! around the interface. Here "distance" is define as the value of the 272! levelset function. If the levelset function is properly defined, 273! this is the true distance normal from the front. Of course, the 274! distance is in a driection normal to the front. 275 do i= 1, npro 276 if (sclr(i) .lt. - epsilon_ls)then 277 prop_blend(i) = zero 278 elseif (abs(sclr(i)) .le. epsilon_ls)then 279 prop_blend(i) = 0.5*(one + Sclr(i)/epsilon_ls + 280 & (sin(pi*Sclr(i)/epsilon_ls))/pi ) 281 elseif (sclr(i) .gt. epsilon_ls) then 282 prop_blend(i) = one 283 endif 284 enddo 285 fact = datmat(1,1,2)/datmat(1,1,1) 286c rsrhol(:) = pres(:) / ( Rgas * T(:) ) 287c call eqs(pres,T,rsrhol) 288c rsrhog(:) = fact*pres(:) / ( Rgas * T(:)) 289 rsrhog(:) = pres(:) / ( Rgas * T(:)) 290 rsrhol(:) = datmat(1,1,1)* 291 & (1+0.000000000517992*(pres-18.02)) 292 rho(:)=rsrhol(:)*prop_blend(:) 293 & +rsrhog(:)*(1-prop_blend(:)) 294c ..... for the VOF case .. in case if we want to run VOF 295c$$$ prop_blend(:) = min((max(sclr(:),0.0)),1.0) 296c$$$ rho(:)=rsrhol(:) * prop_blend(:) + rsrhog(:) * (1-prop_blend(:)) 297c 298 endif 299c 300 ei = yN2 * ( cvs(1) * T 301 & + Rs(1) * Tvib(1) * Texp1 / ( one - Texp1 ) ) 302 & + yO2 * ( cvs(2) * T 303 & + Rs(2) * Tvib(2) * Texp2 / ( one - Texp2 ) ) 304c 305 endif 306c 307 if (ithm .ge. 7) then 308c 309c.... compute enthalpy, cp, alfap, betaT, cv, gamb and c 310c 311 h = ei + Rgas * T 312c 313 alfap = one / T 314c 315 betaT = one / pres 316c 317 cp = yN2 * ( cps(1) + Rs(1) * Tvib(1)**2 * Texp1 318 & / ( ( one - Texp1 ) * T )**2 ) 319 & + yO2 * ( cps(2) + Rs(2) * Tvib(2)**2 * Texp2 320 & / ( ( one - Texp2 ) * T )**2 ) 321c 322 cv = cp - Rgas 323c 324 gamb = Rgas / cv 325c 326 c = sqrt( cp * gamb * T ) 327c 328 endif 329c 330c 331c.... end of ipress = 1 332c 333 endif 334 335c ttim(27) = ttim(27) + secs(0.0) 336c 337c.... end 338c 339 return 340 end 341