!***************************************************************************** !* A collection of two-dimensional Euler numerical fluxes, Version 2 (2010), !* !* written by Dr. Katate Masatsuka (info[at]cfdbooks.com), !* !* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com). !* !*------------------------ !* List of Flux Functions: ! ! Roe ! Rotated-RHLL ! !*------------------------ !* !* These F90 routines were written and made available for download !* for an educational purpose. Detailed descripstion of each numerical !* flux can be found in the original paper, or in popular textbooks, or !* in the (will-be-available) second volume of "I do like CFD". !* !* Note that all routines are not efficietly implemented for clarity; you can !* improve the efficiecy and also you can covert it to double precision !* version if you wish. !* !* NOTES: There were bugs in the Rotated-RHLL solver in version 1. !* It has been fixed and tested for a shock-diffraction problem. !* !* This file may be updated in future. (to add more flux functions.) !***************************************************************************** !***************************************************************************** !* -- Roe's Flux Function --- !* !* P. L. Roe, Approximate Riemann Solvers, Parameter Vectors and Difference !* Schemes, Journal of Computational Physics, 43, pp. 357-372. !* !* Katate Masatsuka, February 2009. http://www.cfdbooks.com !***************************************************************************** function Roe(uL, uR, nx, ny) real :: uL(4), uR(4) ! Input: conservative variables rho*[1, u, v, E] real :: nx, ny ! Input: face normal vector, [nx, ny] (Left-to-Right) real :: Roe(4) ! Output: Roe flux function (upwind) !Local constants real :: gamma ! Ratio of specific heat. real :: zero, fifth, half, one, two ! Numbers !Local variables real :: tx, ty ! Tangent vector (perpendicular to the face normal) real :: vxL, vxR, vyL, vyR ! Velocity components. real :: rhoL, rhoR, pL, pR ! Primitive variables. real :: vnL, vnR, vtL, vtR ! Normal and tangent velocities real :: aL, aR, HL, HR ! Speeds of sound. real :: RT,rho,vx,vy,H,a,vn, vt ! Roe-averages real :: drho,dvx,dvy,dvn,dvt,dp,dV(4) ! Wave strenghs real :: ws(4),dws(4), Rv(4,4) ! Wave speeds and right-eigevectors real :: fL(4), fR(4), diss(4) ! Fluxes ad dissipation term integer :: i, j !Constants. gamma = 1.4 zero = 0.0 fifth = 0.2 half = 0.5 one = 1.0 two = 2.0 !Tangent vector (Do you like it? Actually, Roe flux can be implemented ! without any tangent vector. See "I do like CFD, VOL.1" for details.) tx = -ny ty = nx !Primitive and other variables. ! Left state rhoL = uL(1) vxL = uL(2)/uL(1) vyL = uL(3)/uL(1) vnL = vxL*nx+vyL*ny vtL = vxL*tx+vyL*ty pL = (gamma-one)*( uL(4) - half*rhoL*(vxL*vxL+vyL*vyL) ) aL = sqrt(gamma*pL/rhoL) HL = ( uL(4) + pL ) / rhoL ! Right state rhoR = uR(1) vxR = uR(2)/uR(1) vyR = uR(3)/uR(1) vnR = vxR*nx+vyR*ny vtR = vxR*tx+vyR*ty pR = (gamma-one)*( uR(4) - half*rhoR*(vxR*vxR+vyR*vyR) ) aR = sqrt(gamma*pR/rhoR) HR = ( uR(4) + pR ) / rhoR !First compute the Roe Averages RT = sqrt(rhoR/rhoL) rho = RT*rhoL vx = (vxL+RT*vxR)/(one+RT) vy = (vyL+RT*vyR)/(one+RT) H = ( HL+RT* HR)/(one+RT) a = sqrt( (gamma-one)*(H-half*(vx*vx+vy*vy)) ) vn = vx*nx+vy*ny vt = vx*tx+vy*ty !Wave Strengths drho = rhoR - rhoL dp = pR - pL dvn = vnR - vnL dvt = vtR - vtL dV(1) = (dp - rho*a*dvn )/(two*a*a) dV(2) = rho*dvt/a dV(3) = drho - dp/(a*a) dV(4) = (dp + rho*a*dvn )/(two*a*a) !Wave Speed ws(1) = abs(vn-a) ws(2) = abs(vn) ws(3) = abs(vn) ws(4) = abs(vn+a) !Harten's Entropy Fix JCP(1983), 49, pp357-393: ! only for the nonlinear fields. dws(1) = fifth if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) ) dws(4) = fifth if ( ws(4) < dws(4) ) ws(4) = half * ( ws(4)*ws(4)/dws(4)+dws(4) ) !Right Eigenvectors Rv(1,1) = one Rv(2,1) = vx - a*nx Rv(3,1) = vy - a*ny Rv(4,1) = H - vn*a Rv(1,2) = zero Rv(2,2) = a*tx Rv(3,2) = a*ty Rv(4,2) = vt*a Rv(1,3) = one Rv(2,3) = vx Rv(3,3) = vy Rv(4,3) = half*(vx*vx+vy*vy) Rv(1,4) = one Rv(2,4) = vx + a*nx Rv(3,4) = vy + a*ny Rv(4,4) = H + vn*a !Dissipation Term diss = zero do i=1,4 do j=1,4 diss(i) = diss(i) + ws(j)*dV(j)*Rv(i,j) end do end do !Compute the flux. fL(1) = rhoL*vnL fL(2) = rhoL*vnL * vxL + pL*nx fL(3) = rhoL*vnL * vyL + pL*ny fL(4) = rhoL*vnL * HL fR(1) = rhoR*vnR fR(2) = rhoR*vnR * vxR + pR*nx fR(3) = rhoR*vnR * vyR + pR*ny fR(4) = rhoR*vnR * HR Roe = half * (fL + fR - diss) end function Roe !***************************************************************************** !* -- Rotated-Roe-HLL Flux Function --- !* !* H. Nishikawa and K. Kitamura, Very Simple, Carbuncle-Free, Boundary-Layer !* Resolving, Rotated-Hybrid Riemann Solvers, !* Journal of Computational Physics, 227, pp. 2560-2581, 2008. !* !* The most robust Riemann solver known to the author (in terms of nonlinear !* instability such as carbuncle). !* !* NB: At a boundary face, need to switch to a geometric normal vector: !* (nx2,ny2)=(nx, ny), (nx1,ny1)=(-ny,nx). !* This is not implemented here. It requires information on whether !* the geometric normal, (nx,ny), is on a boundary face or not. !* It shouldn't be difficult for you to implement it. !* !* Katate Masatsuka, February 2010. http://www.cfdbooks.com !***************************************************************************** function Rotated_RHLL(uL, uR, nx, ny) real :: uL(4), uR(4) ! Input: conservative variables rho*[1, u, v, E] real :: nx, ny ! Input: face normal vector, [nx, ny] (Left-to-Right) real :: Rotated_RHLL(4) ! Output: Rotated_RHLL flux function. !Local constants real :: gamma ! Ratio of specific heat. real :: zero, fifth, half, one, two ! Numbers real :: eps ! !Local variables real :: nx1, ny1, nx2, ny2 ! Rotated normals, n1 and n2 real :: tx, ty ! Tangent vector (taken as n1) real :: alpha1, alpha2 ! Projections of the new normals real :: vxL, vxR, vyL, vyR ! Velocity components. real :: rhoL, rhoR, pL, pR ! Primitive variables. real :: vnL, vnR, vtL, vtR ! Normal and tagent velocities real :: aL, aR, HL, HR ! Speeds of sound and total enthalpy real :: RT,rho,vx,vy,H,a ! Roe-averages real :: vn, vt ! Normal and tagent velocities(Roe-average) real :: drho,dvx,dvy,dvn,dvt,dp,dV(4) ! Wave strenghs real :: abs_dq ! Magnitude of the velocity difference real :: abs_ws(4),ws(4),dws(4), Rv(4,4)! Wave speeds and right-eigevectors real :: SRp,SLm ! Wave speeds for the HLL part real :: fL(4), fR(4), diss(4) ! Fluxes ad dissipation term real :: temp integer :: i, j !Constants. gamma = 1.4 zero = 0.0 fifth = 0.2 half = 0.5 one = 1.0 two = 2.0 eps = 1.0e-5 ! 1.0e-12 in the original paper (double precision) !Primitive and other variables. ! Left state rhoL = uL(1) vxL = uL(2)/uL(1) vyL = uL(3)/uL(1) pL = (gamma-one)*( uL(4) - half*rhoL*(vxL*vxL+vyL*vyL) ) aL = sqrt(gamma*pL/rhoL) HL = ( uL(4) + pL ) / rhoL ! Right state rhoR = uR(1) vxR = uR(2)/uR(1) vyR = uR(3)/uR(1) pR = (gamma-one)*( uR(4) - half*rhoR*(vxR*vxR+vyR*vyR) ) aR = sqrt(gamma*pR/rhoR) HR = ( uR(4) + pR ) / rhoR vnL = vxL*nx + vyL*ny vnR = vxR*nx + vyR*ny !Compute the flux. fL(1) = rhoL*vnL fL(2) = rhoL*vnL * vxL + pL*nx fL(3) = rhoL*vnL * vyL + pL*ny fL(4) = rhoL*vnL * HL fR(1) = rhoR*vnR fR(2) = rhoR*vnR * vxR + pR*nx fR(3) = rhoR*vnR * vyR + pR*ny fR(4) = rhoR*vnR * HR !Define n1 and n2, and compute alpha1 and alpha2: (4.2) in the original paper. !(NB: n1 and n2 may need to be frozen at some point during ! a steady calculation to fully make it converge. For time-accurate ! calculation, this is fine.) ! NB: For a boundary face, set (nx2,ny2)=(nx,ny), (nx1,ny1)=(-ny,nx). abs_dq = sqrt( (vxR-vxL)**2+(vyR-vyL)**2 ) if ( abs_dq > eps) then nx1 = (vxR-vxL)/abs_dq ny1 = (vyR-vyL)/abs_dq else nx1 = -ny ny1 = nx endif alpha1 = nx * nx1 + ny * ny1 ! To make alpha1 always positive. temp = sign(one,alpha1) nx1 = temp * nx1 ny1 = temp * ny1 alpha1 = temp * alpha1 ! Take n2 as perpendicular to n1. nx2 = -ny1 ny2 = nx1 alpha2 = nx * nx2 + ny * ny2 ! To make alpha2 always positive. temp = sign(one,alpha2) nx2 = temp * nx2 ny2 = temp * ny2 alpha2 = temp * alpha2 !Now we are going to compute the Roe flux with n2 as the normal !and n1 as the tagent vector, with modified wave speeds (5.12) !Compute the Roe Averages RT = sqrt(rhoR/rhoL) rho = RT*rhoL vx = (vxL+RT*vxR)/(one+RT) vy = (vyL+RT*vyR)/(one+RT) H = ( HL+RT* HR)/(one+RT) a = sqrt( (gamma-one)*(H-half*(vx*vx+vy*vy)) ) vn = vx*nx2+vy*ny2 vt = vx*nx1+vy*ny1 !Wave Strengths (remember that n2 is the normal and n1 is the tangent.) vnL = vxL*nx2 + vyL*ny2 vnR = vxR*nx2 + vyR*ny2 vtL = vxL*nx1 + vyL*ny1 vtR = vxR*nx1 + vyR*ny1 drho = rhoR - rhoL dp = pR - pL dvn = vnR - vnL dvt = vtR - vtL dV(1) = (dp - rho*a*dvn )/(two*a*a) dV(2) = rho*dvt/a dV(3) = drho - dp/(a*a) dV(4) = (dp + rho*a*dvn )/(two*a*a) !Wave Speeds for Roe flux part. ws(1) = vn-a ws(2) = vn ws(3) = vn ws(4) = vn+a abs_ws = abs(ws) !Harten's Entropy Fix JCP(1983), 49, pp357-393: !only for the nonlinear fields. dws(1) = fifth if (abs_ws(1)