Specify the density in the ghost cell based on the temperature on the wall. Isothermal or adiabatic
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in), | dimension(1:6) | :: | temperature | ||
| character(len=*), | intent(in) | :: | face | |||
| type(boundarytype), | intent(in) | :: | bc | |||
| type(extent), | intent(in) | :: | dims |
subroutine temp_based_density(temperature, face, bc, dims)
!< Specify the density in the ghost cell based on the
!< temperature on the wall. Isothermal or adiabatic
implicit none
type(extent), intent(in) :: dims
type(boundarytype), intent(in) :: bc
real(wp), dimension(1:6) , intent(in) :: temperature
character(len=*) , intent(in) :: face
real(wp) :: stag_temp
integer :: i,j,k
select case(face)
case("imin")
if(temperature(1)<0.0)then
do k = 1,kmx-1
do j = 1,jmx-1
do i = 1,1
stag_temp = (pressure(i,j,k)/(R_gas*density(i,j,k)))*(1 + (0.5*(gm-1.)*gm*pressure(i,j,k)/density(i,j,k)))
density(i-1,j,k) = pressure(i-1,j,k)/(R_gas*stag_temp)
density(i-2,j,k) = pressure(i-2,j,k)/(R_gas*stag_temp)
density(i-3,j,k) = pressure(i-3,j,k)/(R_gas*stag_temp)
end do
end do
end do
elseif(temperature(1)>1.0)then
do k = 1,kmx-1
do j = 1,jmx-1
do i = 1,1
density(i-1,j,k) = pressure(i-1,j,k)/(R_gas*(2*temperature(1)-(pressure(i+0,j,k)/(R_gas*density(i+0,j,k)))))
density(i-2,j,k) = pressure(i-2,j,k)/(R_gas*(2*temperature(1)-(pressure(i+1,j,k)/(R_gas*density(i+1,j,k)))))
density(i-3,j,k) = pressure(i-3,j,k)/(R_gas*(2*temperature(1)-(pressure(i+2,j,k)/(R_gas*density(i+2,j,k)))))
end do
end do
end do
else
call copy3(density , "symm", face, bc, dims)
end if
case("imax")
if(temperature(2)<0.0)then
do k = 1,kmx-1
do j = 1,jmx-1
do i = imx-1,imx-1
stag_temp = (pressure(i,j,k)/(R_gas*density(i,j,k)))*(1 + (0.5*(gm-1.)*gm*pressure(i,j,k)/density(i,j,k)))
density(i+1,j,k) = pressure(i+1,j,k)/(R_gas*stag_temp)
density(i+2,j,k) = pressure(i+2,j,k)/(R_gas*stag_temp)
density(i+3,j,k) = pressure(i+3,j,k)/(R_gas*stag_temp)
end do
end do
end do
elseif(temperature(2)>1.0)then
do k = 1,kmx-1
do j = 1,jmx-1
do i = imx-1,imx-1
density(i+1,j,k) = pressure(i+1,j,k)/(R_gas*(2*temperature(2)-(pressure(i-0,j,k)/(R_gas*density(i-0,j,k)))))
density(i+2,j,k) = pressure(i+2,j,k)/(R_gas*(2*temperature(2)-(pressure(i-1,j,k)/(R_gas*density(i-1,j,k)))))
density(i+3,j,k) = pressure(i+3,j,k)/(R_gas*(2*temperature(2)-(pressure(i-2,j,k)/(R_gas*density(i-2,j,k)))))
end do
end do
end do
else
call copy3(density , "symm", face, bc, dims)
end if
case("jmin")
if(temperature(3)<0.0)then
do k = 1,kmx-1
do j = 1,1
do i = 1,imx-1
stag_temp = (pressure(i,j,k)/(R_gas*density(i,j,k)))*(1 + (0.5*(gm-1.)*gm*pressure(i,j,k)/density(i,j,k)))
density(i,j-1,k) = pressure(i,j-1,k)/(R_gas*stag_temp)
density(i,j-2,k) = pressure(i,j-2,k)/(R_gas*stag_temp)
density(i,j-3,k) = pressure(i,j-3,k)/(R_gas*stag_temp)
end do
end do
end do
elseif(temperature(3)>1.0)then
do k = 1,kmx-1
do j = 1,1
do i = 1,imx-1
density(i,j-1,k) = pressure(i,j-1,k)/(R_gas*(2*temperature(3)-(pressure(i,j+0,k)/(R_gas*density(i,j+0,k)))))
density(i,j-2,k) = pressure(i,j-2,k)/(R_gas*(2*temperature(3)-(pressure(i,j+1,k)/(R_gas*density(i,j+1,k)))))
density(i,j-3,k) = pressure(i,j-3,k)/(R_gas*(2*temperature(3)-(pressure(i,j+2,k)/(R_gas*density(i,j+2,k)))))
end do
end do
end do
else
call copy3(density , "symm", face, bc, dims)
end if
case("jmax")
if(temperature(4)<0.0)then
do k = 1,kmx-1
do j = jmx-1,jmx-1
do i = 1,imx-1
stag_temp = (pressure(i,j,k)/(R_gas*density(i,j,k)))*(1 + (0.5*(gm-1.)*gm*pressure(i,j,k)/density(i,j,k)))
density(i,j+1,k) = pressure(i,j+1,k)/(R_gas*stag_temp)
density(i,j+2,k) = pressure(i,j+2,k)/(R_gas*stag_temp)
density(i,j+3,k) = pressure(i,j+3,k)/(R_gas*stag_temp)
end do
end do
end do
elseif(temperature(4)>1.0)then
do k = 1,kmx-1
do j = jmx-1,jmx-1
do i = 1,imx-1
density(i,j+1,k) = pressure(i,j+1,k)/(R_gas*(2*temperature(4)-(pressure(i,j-0,k)/(R_gas*density(i,j-0,k)))))
density(i,j+2,k) = pressure(i,j+2,k)/(R_gas*(2*temperature(4)-(pressure(i,j-1,k)/(R_gas*density(i,j-1,k)))))
density(i,j+3,k) = pressure(i,j+3,k)/(R_gas*(2*temperature(4)-(pressure(i,j-2,k)/(R_gas*density(i,j-2,k)))))
end do
end do
end do
else
call copy3(density , "symm", face, bc, dims)
end if
case("kmin")
if(temperature(5)<0.0)then
do k = 1,1
do j = 1,jmx-1
do i = 1,imx-1
stag_temp = (pressure(i,j,k)/(R_gas*density(i,j,k)))*(1 + (0.5*(gm-1.)*gm*pressure(i,j,k)/density(i,j,k)))
density(i,j,k-1) = pressure(i,j,k-1)/(R_gas*stag_temp)
density(i,j,k-2) = pressure(i,j,k-2)/(R_gas*stag_temp)
density(i,j,k-3) = pressure(i,j,k-3)/(R_gas*stag_temp)
end do
end do
end do
elseif(temperature(5)>1.0)then
do k = 1,1
do j = 1,jmx-1
do i = 1,imx-1
density(i,j,k-1) = pressure(i,j,k-1)/(R_gas*(2*temperature(5)-(pressure(i,j,k+0)/(R_gas*density(i,j,k+0)))))
density(i,j,k-2) = pressure(i,j,k-2)/(R_gas*(2*temperature(5)-(pressure(i,j,k+1)/(R_gas*density(i,j,k+1)))))
density(i,j,k-3) = pressure(i,j,k-3)/(R_gas*(2*temperature(5)-(pressure(i,j,k+2)/(R_gas*density(i,j,k+2)))))
end do
end do
end do
else
call copy3(density , "symm", face, bc, dims)
end if
case("kmax")
if(temperature(6)<0.0)then
do k = kmx-1,kmx-1
do j = 1,jmx-1
do i = 1,imx-1
stag_temp = (pressure(i,j,k)/(R_gas*density(i,j,k)))*(1 + (0.5*(gm-1.)*gm*pressure(i,j,k)/density(i,j,k)))
density(i,j,k+1) = pressure(i,j,k+1)/(R_gas*stag_temp)
density(i,j,k+2) = pressure(i,j,k+2)/(R_gas*stag_temp)
density(i,j,k+3) = pressure(i,j,k+3)/(R_gas*stag_temp)
end do
end do
end do
elseif(temperature(6)>1.0)then
do k = kmx-1,kmx-1
do j = 1,jmx-1
do i = 1,imx-1
density(i,j,k+1) = pressure(i,j,k+1)/(R_gas*(2*temperature(6)-(pressure(i,j,k-0)/(R_gas*density(i,j,k-0)))))
density(i,j,k+2) = pressure(i,j,k+2)/(R_gas*(2*temperature(6)-(pressure(i,j,k-1)/(R_gas*density(i,j,k-1)))))
density(i,j,k+3) = pressure(i,j,k+3)/(R_gas*(2*temperature(6)-(pressure(i,j,k-2)/(R_gas*density(i,j,k-2)))))
end do
end do
end do
else
call copy3(density , "symm", face, bc, dims)
end if
case DEFAULT
!print*, "ERROR: wrong face for boundary condition"
Fatal_error
end select
end subroutine temp_based_density