!******************************************************************************** !* Educationally-Designed Unstructured 2D (EDU2D) Code !* !* --- EDU2D_MMS_TWODNS !* !* !* This program shows how to compute the source terms arising from the method !* of manufactured solutions for the 2D compressible Navier-Stokes equations. !* !* !* Method of manufactured solutions: !* !* (0) Consider a target equation, df(u)/dx + dg(u)/dy = 0. !* (1) Choose a function, v(x,y). !* (2) Substitute v into the tartget equation: s = df(v)/dx + dg(v)/dy. !* (3) Then, v(x,y) is the exact solution for the following eq.: !* !* df(u)/dx + dg(u)/dy = s. !* !* (4) Numerically solve this equation for u, then you can !* measure the error: error = |u - v|. !* !* This program shows how to evaluate the source term 's' numerically !* for the 2D Navier-Stokes equation, for a chosen function v(x,y), !* which is given here by !* !* v(x,y) = c0 + cs*sin(cx*x+cy*y), !* !* where the coefficients, c0, cs, cx, cy, can be chosen differently !* for different variables. You can change this function as you like. !* !* !* Note: The main program calls the subroutine: !* !* 'compute_manufactured_sol_and_source' !* !* for a specified (x,y), and prints on screen the solutions and !* source terms. !* !* Note: You can use the subroutine 'compute_manufactured_sol_and_source' !* to compute the manufactured solution and the source term in !* your code. !* !* Note: The source (or forcing) term, 's', is a function of space only, not !* of the solution. Therefore, it is computed only once and stored. !* !* Note: It is quite straightforward to extend the subroutine to generate !* source terms for the 3D Navier-Stokes equations. !* !* Note: The viscosity is computed by Sutherland's law with a suitable !* scaling for the Navier-Stokes equations nondimensionalized with !* the freestream speed of sound: (u,v)=(u/a_inf,v/a_inf), !* p=p/(rho_inf*a_inf^2). Modify the factor if other (or no) !* non-dimensionalization is used. See "I do like CFD, VOL.1" for !* details, which is available in PDF at http://www.cfdbooks.com. !* !* !* This is Version 2 (December 16, 2017). !* !* This F90 code is written and made available for an educational purpose. !* This file may be updated in future. !* !* Katate Masatsuka, http://www.cfdbooks.com !******************************************************************************** program edu2d_mms_twodns implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2) :: T_inf, M_inf, Re_inf, gamma, Prandtl real(p2) :: x, y ! 1 2 3 4 5 6 7 8 9 10 11 ! Solution: w(1:11) = [rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort, T] real(p2), dimension(11) :: w ! Gradients: gradw(1:11,1) = d[rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort,T]/dx ! gradw(1:11,2) = d[rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort,T]/dy real(p2), dimension(11,2) :: gradw ! Second derivatives: wxx = d^2[rho,u,v,p]/dx^2 ! wyy = d^2[rho,u,v,p]/dy^2 ! wyy = d^2[rho,u,v,p]/dx/dy real(p2), dimension( 4) :: wxx, wyy, wxy ! Source terms: Source terms for [continuity, x/y-momentum, energy eqns] real(p2), dimension( 4) :: source !------------------------------------------------------------------------- ! Exmaple parameter values: T_inf = 300.0_p2 ! Freestream temperature M_inf = 0.3_p2 ! Freestream Mach number Re_inf = 1000.0_p2 ! Freestream Reynolds number gamma = 1.4_p2 ! Ratio of specific heats Prandtl = 0.72_p2 ! Prandtl number !------------------------------------------------------------------------- ! The position where the manufactured (exact) soluton ! and the source terms are sought. Below is just an example. x = 0.42113973644362_p2 y = 0.81842548218445_p2 !------------------------------------------------------------------------- ! Compute the manufactured solution and source terms at (x,y): call compute_manufactured_sol_and_source(x,y, T_inf,M_inf,Re_inf,gamma, & Prandtl, w,gradw,wxx,wyy,wxy,source ) write(*,*) write(*,*) " ****************************************************" write(*,*) " Position where the solution and source is sought:" write(*,*) " ****************************************************" write(*,*) write(*,*) " x = ", x write(*,*) " y = ", y write(*,*) write(*,*) " ****************************************************" write(*,*) " Manufactured solution (exact solution):" write(*,*) " ****************************************************" write(*,*) write(*,*) " rho = ", w( 1) write(*,*) " u = ", w( 2) write(*,*) " v = ", w( 3) write(*,*) " p = ", w( 4) write(*,*) " tauxx = ", w( 5) write(*,*) " tauyy = ", w( 6) write(*,*) " tauxy = ", w( 7) write(*,*) " qx = ", w( 8) write(*,*) " qy = ", w( 9) write(*,*) " vorticity = ", w(10) write(*,*) " temperature = ", w(11) write(*,*) write(*,*) " --------------------------------------" write(*,*) " First derivatives:" write(*,*) write(*,*) " rhox = ", gradw( 1,1) write(*,*) " ux = ", gradw( 2,1) write(*,*) " vx = ", gradw( 3,1) write(*,*) " px = ", gradw( 4,1) write(*,*) " (tauxx)x = ", gradw( 5,1) write(*,*) " (tauyy)x = ", gradw( 6,1) write(*,*) " (tauxy)x = ", gradw( 7,1) write(*,*) " qxx = ", gradw( 8,1) write(*,*) " qyx = ", gradw( 9,1) write(*,*) " (vorticity)x = ", gradw(10,1) write(*,*) " (temperature)x = ", gradw(11,1) write(*,*) write(*,*) " rhoy = ", gradw( 1,2) write(*,*) " uy = ", gradw( 2,2) write(*,*) " vy = ", gradw( 3,2) write(*,*) " py = ", gradw( 4,2) write(*,*) " (tauxx)y = ", gradw( 5,2) write(*,*) " (tauyy)y = ", gradw( 6,2) write(*,*) " (tauxy)y = ", gradw( 7,2) write(*,*) " qxy = ", gradw( 8,2) write(*,*) " qyy = ", gradw( 9,2) write(*,*) " (vorticity)y = ", gradw(10,2) write(*,*) " (temperature)y = ", gradw(11,2) write(*,*) write(*,*) " --------------------------------------" write(*,*) " Second derivatives of [rho,u,v,p]:" write(*,*) write(*,*) " rhoxx = ", wxx(1) write(*,*) " uxx = ", wxx(2) write(*,*) " vxx = ", wxx(3) write(*,*) " pxx = ", wxx(4) write(*,*) write(*,*) " rhoyy = ", wyy(1) write(*,*) " uyy = ", wyy(2) write(*,*) " vyy = ", wyy(3) write(*,*) " pyy = ", wyy(4) write(*,*) write(*,*) " rhoxy = ", wxy(1) write(*,*) " uxy = ", wxy(2) write(*,*) " vxy = ", wxy(3) write(*,*) " pxy = ", wxy(4) write(*,*) write(*,*) " ****************************************************" write(*,*) " Source terms:" write(*,*) " ****************************************************" write(*,*) write(*,*) " Continuity: s(1) = ", source(1) write(*,*) " X momemtum: s(2) = ", source(2) write(*,*) " Y momemtum: s(3) = ", source(3) write(*,*) " Energy : s(4) = ", source(4) write(*,*) write(*,*) " Add these to the Navier-Stokes equations: fx+gy=0," write(*,*) " on the right hand side, and the manufactured solutions" write(*,*) " are the exact solution to the resulting system, fx+gy=s." write(*,*) stop contains !******************************************************************************** !* !* Educationally-Designed Unstructured 2D (EDU2D) Code !* !* --- EDU2D_MMS !* !* !* This subroutine computes, for a given position (x,y), the manufactured !* solution values and the source terms. The manufactured solution takes !* the following form: !* !* c0 + cs*sin(cx*x+cy*y) !* !* with the coefficients, c0, cs, cx, cy chosen differently for different !* variable. !* !* ------------------------ !* Input: !* !* (x,y) = Position at which the solution and source terms are sought. !* !* ------------------------ !* Output: !* 1 2 3 4 5 6 7 8 9 10 11 !* w(1:11) = [rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort, T] !* !* gradw(1:11,1) = d[rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort, T]/dx !* !* gradw(1:11,2) = d[rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort, T]/dy !* !* source(1:5) = Source terms for continuity, x/y-momentum, and energy eqns !* !* ------------------------ !* !* Note: Second derivatives are computed only for [rho,u,v,p] !* !* uxx, uyy, uxy, vxx, vyy, vxy, rhoxx, rhoyy, rhoxy, pxx, pyy, pxy. !* !* just because these are required anyway to compute the first-derivatives !* of [tauxx,tauyy,tauxy,qx,qy,vort]. !* !* Note: You can extend this subroutine to compute higher-order derivatives !* of the solution and the source terms as well. Personally, I have !* an extended version that computes up to second derivatives of the !* source terms. !* !* Note: You can extend this also to the unsteady NS equations by including !* time derivatives when generating the source term: s=dv/dt+df(v)/dx+dg(v)/dy. !* !* Note: For high-order schemes, you can use this subroutine to compute the !* source terms at quadrature points. !* !* !* This is Version 2 (December 16, 2017). !* !* This F90 code is written and made available for an educational purpose. !* This file may be updated in future. !* !* Katate Masatsuka, http://www.cfdbooks.com !******************************************************************************** subroutine compute_manufactured_sol_and_source(x,y,T_inf,M_inf,Re_inf,gamma, & Prandtl, w,gradw,wxx,wyy,wxy,source ) implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2), parameter :: pi = 3.141592653589793238_p2 real(p2), parameter :: half = 0.5_p2 real(p2), parameter :: one = 1.0_p2 real(p2), parameter :: two = 2.0_p2 real(p2), parameter :: two_third = 2.0_p2/3.0_p2 real(p2), parameter :: four_third = 4.0_p2/3.0_p2 !Input real(p2), intent( in) :: x, y real(p2), intent( in) :: T_inf, M_inf, Re_inf, gamma, Prandtl !Output real(p2), dimension(11 ), intent(out) :: w real(p2), dimension(11,2), intent(out) :: gradw real(p2), dimension( 4 ), intent(out) :: wxx, wyy, wxy real(p2), dimension( 4 ), intent(out) :: source !Local variables real(p2) :: rho , u , v , p real(p2) :: rhox, ux, vx, px real(p2) :: rhoy, uy, vy, py real(p2) :: mu, mux, muy real(p2) :: t, tx, ty real(p2) :: a, ax, ay real(p2) :: b, bx, by real(p2) :: k, kx, ky real(p2) :: u2 , u2x , u2y real(p2) :: v2 , v2x , v2y real(p2) :: au , aux , auy real(p2) :: av , avx , avy real(p2) :: bu , bux , buy real(p2) :: bv , bvx , bvy real(p2) :: rhoH , rhoHx , rhoHy real(p2) :: rhouH, rhouHx, rhouHy real(p2) :: rhovH, rhovHx, rhovHy real(p2) :: rhoxx, rhoxy, rhoyy real(p2) :: uxx, uxy, uyy real(p2) :: vxx, vxy, vyy real(p2) :: pxx, pxy, pyy real(p2) :: axx, axy, ayy real(p2) :: bxx, bxy, byy real(p2) :: txx, txy, tyy integer :: ix, iy integer :: irho, iu, iv, ip integer :: itauxx, itauyy, itauxy integer :: iqx, iqy integer :: ivort, iTemp integer :: iconti, ixmom, iymom, ienrgy real(p2) :: cr0, crs, crx, cry real(p2) :: cu0, cus, cux, cuy real(p2) :: cv0, cvs, cvx, cvy real(p2) :: cp0, cps, cpx, cpy !----------------------------------------------------------- ! Define some indices ix = 1 ! x component iy = 2 ! y component irho = 1 ! density iu = 2 ! x-velocity iv = 3 ! y-velocity ip = 4 ! pressure itauxx = 5 ! viscous stress, xx itauyy = 6 ! viscous stress, yy itauxy = 7 ! viscous stress, xy iqx = 8 ! heat flux in x iqy = 9 ! heat flux in y ivort = 10 ! vorticity iTemp = 11 ! temperature iconti = 1 ! continuity equation ixmom = 2 ! x-momentum equation iymom = 3 ! y-momentum equation ienrgy = 4 ! energy equation !----------------------------------------------------------- ! Constants for the exact solution: c0 + cs*sin(cx*x+cy*y). ! ! Note: Make sure the density and pressure are positive. ! Note: These values are passed to the subroutine, ! manufactured_sol(c0,cs,cx,cy, nx,ny,x,y), ! whcih returns the solution value or derivative. !----------------------------------------- ! Density = cr0 + crs*sin(crx*x+cry*y) cr0 = 1.12_p2 crs = 0.15_p2 crx = 3.12_p2*pi cry = 2.92_p2*pi !----------------------------------------- ! X-velocity = cu0 + cus*sin(cux*x+cuy*y) cu0 = 1.32_p2 cus = 0.06_p2 cux = 2.09_p2*pi cuy = 3.12_p2*pi !----------------------------------------- ! Y-velocity = cv0 + cvs*sin(cvx*x+cvy*y) cv0 = 1.18_p2 cvs = 0.03_p2 cvx = 2.15_p2*pi cvy = 3.32_p2*pi !----------------------------------------- ! Pressure = cp0 + cps*sin(cpx*x+cpy*y) cp0 = 1.62_p2 cps = 0.31_p2 cpx = 3.79_p2*pi cpy = 2.98_p2*pi !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Part I: Compute w = [rho,u,v,p,tauxx,tauyy,tauxy,qx,qy,vort,T] and grad(w). !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !------------------------------------------------------------------------ ! rho: Density and its derivatives w(irho) = manufactured_sol(cr0,crs,crx,cry, 0,0,x,y) gradw(irho,ix) = manufactured_sol(cr0,crs,crx,cry, 1,0,x,y) gradw(irho,iy) = manufactured_sol(cr0,crs,crx,cry, 0,1,x,y) rhoxx = manufactured_sol(cr0,crs,crx,cry, 2,0,x,y) rhoyy = manufactured_sol(cr0,crs,crx,cry, 0,2,x,y) rhoxy = manufactured_sol(cr0,crs,crx,cry, 1,1,x,y) wxx(irho) = rhoxx wyy(irho) = rhoyy wxy(irho) = rhoxy !------------------------------------------------------------------------ ! u: x-velocity and its derivatives w(iu) = manufactured_sol(cu0,cus,cux,cuy, 0,0,x,y) gradw(iu,ix) = manufactured_sol(cu0,cus,cux,cuy, 1,0,x,y) gradw(iu,iy) = manufactured_sol(cu0,cus,cux,cuy, 0,1,x,y) uxx = manufactured_sol(cu0,cus,cux,cuy, 2,0,x,y) uyy = manufactured_sol(cu0,cus,cux,cuy, 0,2,x,y) uxy = manufactured_sol(cu0,cus,cux,cuy, 1,1,x,y) wxx(iu) = uxx wyy(iu) = uyy wxy(iu) = uxy !------------------------------------------------------------------------ ! v: y-velocity and its derivatives w(iv) = manufactured_sol(cv0,cvs,cvx,cvy, 0,0,x,y) gradw(iv,ix) = manufactured_sol(cv0,cvs,cvx,cvy, 1,0,x,y) gradw(iv,iy) = manufactured_sol(cv0,cvs,cvx,cvy, 0,1,x,y) vxx = manufactured_sol(cv0,cvs,cvx,cvy, 2,0,x,y) vyy = manufactured_sol(cv0,cvs,cvx,cvy, 0,2,x,y) vxy = manufactured_sol(cv0,cvs,cvx,cvy, 1,1,x,y) wxx(iv) = vxx wyy(iv) = vyy wxy(iv) = vxy !------------------------------------------------------------------------ ! p: pressure and its derivatives w(ip) = manufactured_sol(cp0,cps,cpx,cpy, 0,0,x,y) gradw(ip,ix) = manufactured_sol(cp0,cps,cpx,cpy, 1,0,x,y) gradw(ip,iy) = manufactured_sol(cp0,cps,cpx,cpy, 0,1,x,y) pxx = manufactured_sol(cp0,cps,cpx,cpy, 2,0,x,y) pyy = manufactured_sol(cp0,cps,cpx,cpy, 0,2,x,y) pxy = manufactured_sol(cp0,cps,cpx,cpy, 1,1,x,y) wxx(ip) = pxx wyy(ip) = pyy wxy(ip) = pxy !----------------------------------------------------------------------------- ! Store the exact solution and derivatives in local variables. rho = w(irho) u = w(iu) v = w(iv) p = w(ip) rhox = gradw(irho,ix) ux = gradw(iu ,ix) vx = gradw(iv ,ix) px = gradw(ip ,ix) rhoy = gradw(irho,iy) uy = gradw(iu ,iy) vy = gradw(iv ,iy) py = gradw(ip ,iy) !------------------------------------------------------------------------ ! w: Vorticity and its derivatives. w(ivort) = vx - uy gradw(ivort,ix) = vxx - uxy gradw(ivort,iy) = vxy - uyy !------------------------------------------------------------------------ ! T: Temperature and its derivatives. T=(gamma*p)*(1/rho). a = gamma*p ax = gamma*px ay = gamma*py axx = gamma*pxx axy = gamma*pxy ayy = gamma*pyy b = rho**(-1) bx = -rho**(-2)*rhox by = -rho**(-2)*rhoy bxx = two*rho**(-3)*rhox*rhox - rho**(-2)*rhoxx bxy = two*rho**(-3)*rhox*rhoy - rho**(-2)*rhoxy byy = two*rho**(-3)*rhoy*rhoy - rho**(-2)*rhoyy !Note: This subroutine returns t=a*b, tx=ax*b+a*bx, and ty=ay*b+a*by, ! and also the second derivatives. call derivatives_ab2(a,ax,ay,axx,ayy,axy, b,bx,by,bxx,byy,bxy, & t,tx,ty,txx,tyy,txy ) w(iTemp) = t gradw(iTemp,ix) = tx gradw(iTemp,iy) = ty !------------------------------------ ! mu: Viscosity and its derivatives. mu = viscosity( t,T_inf,M_inf,Re_inf) mux = dviscosity_dT(t,T_inf,M_inf,Re_inf) * tx ! = mux = dmu/dT*dT/dx muy = dviscosity_dT(t,T_inf,M_inf,Re_inf) * ty ! = muy = dmu/dT*dT/dy !------------------------------------ !------------------------------------------------------------------------ ! tauxx: viscous stress tauxx and its derivatives b = four_third*ux - two_third*vy bx = four_third*uxx - two_third*vxy by = four_third*uxy - two_third*vyy !Note: This subroutine returns k=mu*b, kx=mux*b+mu*bx, and ky=muy*b+mu*by. call derivatives_ab(mu,mux,muy, b,bx,by, k,kx,ky) w(itauxx) = mu*( four_third*ux - two_third*vy ) ! = k gradw(itauxx,ix) = kx gradw(itauxx,iy) = ky !------------------------------------------------------------------------ ! tauyy: viscous stress tauyy and its derivatives b = four_third*vy - two_third*ux bx = four_third*vxy - two_third*uxx by = four_third*vyy - two_third*uxy call derivatives_ab(mu,mux,muy, b,bx,by, k,kx,ky) w(itauyy) = mu*( four_third*vy - two_third*ux ) ! = k gradw(itauyy,ix) = kx gradw(itauyy,iy) = ky !------------------------------------------------------------------------ ! tauxy: viscous stress tauxy and its derivatives b = uy + vx bx = uxy + vxx by = uyy + vxy call derivatives_ab(mu,mux,muy, b,bx,by, k,kx,ky) w(itauxy) = mu*( uy + vx ) ! = k gradw(itauxy,ix) = kx gradw(itauxy,iy) = ky !------------------------------------------------------------------------ ! qx: Heat flux qx and its derivatives b = -tx /(Prandtl*(gamma-one)) bx = -txx /(Prandtl*(gamma-one)) by = -txy /(Prandtl*(gamma-one)) call derivatives_ab(mu,mux,muy, b,bx,by, k,kx,ky) w(iqx) = -mu*tx/(Prandtl*(gamma-one)) ! = k gradw(iqx,ix) = kx gradw(iqx,iy) = ky !------------------------------------------------------------------------ ! qy: Heat flux qy and its derivatives b = -ty /(Prandtl*(gamma-one)) bx = -txy /(Prandtl*(gamma-one)) by = -tyy /(Prandtl*(gamma-one)) call derivatives_ab(mu,mux,muy, b,bx,by, k,kx,ky) w(iqy) = -mu*ty/(Prandtl*(gamma-one)) ! = k gradw(iqy,ix) = kx gradw(iqy,iy) = ky !----------------------------------------------------------------------------- ! ! Exact (manufactured) solutons and derivatives have been computed. ! We move on the source terms, which are the governing equations we wish ! to solve (the compressible NS) evaluated by the exact solution (i.e., ! the function we wish to make exact with the source terms). ! ! Navier-Stokes: dF(w)/dx + dG(w)/dy = 0. ! ! Our manufactured solutions, wm, are the exact soltuions to the following: ! ! dF(w)/dx + dG(w)/dy = S, ! ! where S = dF(wm)/dx + dG(wm)/dy. ! !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Part II: Compute the source terms. !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Inviscid terms !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Derivatives of u^2 call derivatives_ab(u,ux,uy, u,ux,uy, u2,u2x,u2y) ! a=u^2 ! Derivatives of v^2 call derivatives_ab(v,vx,vy, v,vx,vy, v2,v2x,v2y) ! b=v^2 ! Derivatives of k=(u^2+v^2)/2 k = half*(u*u + v*v) kx = half*(u2x + v2x) ky = half*(u2y + v2y) ! Derivatives of rho*k = rho*(u^2+v^2)/2 call derivatives_ab(rho,rhox,rhoy, k,kx,ky, a,ax,ay) !a=rho*(u^2+v^2)/2 ! Derivatives of rho*H = gamma/(gamma-1)*p + rho*k rhoH = gamma/(gamma-one)*p + a rhoHx = gamma/(gamma-one)*px + ax rhoHy = gamma/(gamma-one)*py + ay !----------------------------------------------------------------------------- ! Compute derivatives of (rho*u) call derivatives_ab(rho,rhox,rhoy, u,ux,uy, a,ax,ay) !a=(rho*u) ! Compute derivatives of (rho*v) call derivatives_ab(rho,rhox,rhoy, v,vx,vy, b,bx,by) !ab=(rho*v) !----------------------------------------------------------------------------- ! Compute derivatives of (a*u)=(rho*u*u) call derivatives_ab(a,ax,ay, u,ux,uy, au,aux,auy) ! Compute derivatives of (a*v)=(rho*u*v) call derivatives_ab(a,ax,ay, v,vx,vy, av,avx,avy) ! Compute derivatives of (b*u)=(rho*v*u) call derivatives_ab(b,bx,by, u,ux,uy, bu,bux,buy) ! Compute derivatives of (b*v)=(rho*v*v) call derivatives_ab(b,bx, by, v,vx,vy, bv,bvx,bvy) !----------------------------------------------------------------------------- ! Compute derivatives of (u*rH) call derivatives_ab( u,ux,uy, rhoH,rhoHx,rhoHy, rhouH,rhouHx,rhouHy) ! Compute derivatives of (v*rH) call derivatives_ab( v,vx,vy, rhoH,rhoHx,rhoHy, rhovH,rhovHx,rhovHy) !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! Store the inviscid terms in the source term array, source(:). !--------------------------------------------------------------------- !------------------------------------------------------ ! Continuity: (rho*u)_x + (rho*v)_y source(iconti) = (rhox*u + rho*ux) + (rhoy*v + rho*vy) !------------------------------------------------------ ! Momentum: (rho*u*u)_x + (rho*u*v)_y + px source(ixmom) = aux + buy + px !------------------------------------------------------ ! Momentum: (rho*u*v)_x + (rho*v*v)_y + px source(iymom) = avx + bvy + py !------------------------------------------------------ ! Energy: (rho*u*H)_x + (rho*v*H) source(ienrgy) = rhouHx + rhovHy !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Viscous terms !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------- ! Derivatives of (tauxx*u) a = w(itauxx) ax = gradw(itauxx,ix) ay = gradw(itauxx,iy) call derivatives_ab(a,ax,ay, u,ux,uy, au,aux,auy) !----------------------------------------------------------- ! Derivatives of (tauyy*v) a = w(itauyy) ax = gradw(itauyy,ix) ay = gradw(itauyy,iy) call derivatives_ab(a,ax,ay, v,vx,vy, bv,bvx,bvy) !----------------------------------------------------------- ! Derivatives of (tauxy*v) and (tauxy*u) a = w(itauxy) ax = gradw(itauxy,ix) ay = gradw(itauxy,iy) call derivatives_ab(a,ax,ay, v,vx,vy, bu,bux,buy) call derivatives_ab(a,ax,ay, u,ux,uy, av,avx,avy ) !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! Add viscous terms to the inviscid terms. !--------------------------------------------------------------------- !--------------------------------------------------------------------- !------------------------------------------------------ ! X-momentum: (u*u)_x + (u*v)_y + px - [ (tauxx)_x + (tauxy)_y ] source(ixmom ) = source(ixmom ) - ( gradw(itauxx,ix) + gradw(itauxy,iy)) !------------------------------------------------------ ! Y-momentum: (u*v)_x + (v*v)_y + px - [ (tauxy)_x + (tauyy)_y ] source(iymom ) = source(iymom ) - ( gradw(itauxy,ix) + gradw(itauyy,iy)) !------------------------------------------------------ ! Energy: (u*v)_x + (v*v)_y + px - (tauxx*u+tauxy*v)_x - (tauxy*u+tauyy*v)_y + (qx)_x + (qy)_y source(ienrgy) = source(ienrgy) - (aux + bux) - (avy + bvy) + gradw(iqx,ix) + gradw(iqy,iy) end subroutine compute_manufactured_sol_and_source !******************************************************************************** !******************************************************************************** !******************************************************************************** ! ! This subroutine computes 1st and 2nd derivatives of a quadratic term ! ! Input: a, ax, ay !Function value a, and its derivatives, (ax,ay,axx,ayy,axy). ! b, bx, by !Function value b, and its derivatives, (bx,by,bxx,byy,bxy). ! ! Output: ab = a*b, abx = d(a*b)/dx, aby = d(a*b)/dy, and second derivatives. ! !******************************************************************************** subroutine derivatives_ab2(a,ax,ay,axx,ayy,axy, b,bx,by,bxx,byy,bxy, & ab,abx,aby,abxx,abyy,abxy ) implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2), parameter :: two = 2.0_p2 !Input real(p2), intent( in) :: a, ax, ay, axx, ayy, axy real(p2), intent( in) :: b, bx, by, bxx, byy, bxy !Output real(p2), intent(out) :: ab, abx, aby, abxx, abyy, abxy ab = a*b abx = ax*b + a*bx aby = ay*b + a*by abxx = axx*b + two*ax*bx + a*bxx abyy = ayy*b + two*ay*by + a*byy abxy = axy*b + ax*by + ay*bx + a*bxy end subroutine derivatives_ab2 !******************************************************************************** ! ! This subroutine computes first derivatives of a quadratic term ! ! Input: a, ax, ay !Function value a, and its derivatives, (ax,ay). ! b, bx, by !Function value b, and its derivatives, (bx,by). ! ! Output: ab = a*b, abx = d(a*b)/dx, aby = d(a*b)/dy. ! !******************************************************************************** subroutine derivatives_ab(a,ax,ay, b,bx,by, ab,abx,aby) implicit none integer , parameter :: p2 = selected_real_kind(p=15) !Input real(p2), intent( in) :: a, ax, ay real(p2), intent( in) :: b, bx, by !Output real(p2), intent(out) :: ab, abx, aby ab = a*b abx = ax*b + a*bx aby = ay*b + a*by end subroutine derivatives_ab !******************************************************************************** !* !* This function computes the sine function: !* !* f = a0 + as*sin(ax*x+ay*y) !* !* and its derivatives: !* !* df/dx^nx/dy^ny = d^{nx+ny}(a0+as*sin(ax*x+ay*y))/(dx^nx*dy^ny) !* !* depending on the input parameters: !* !* !* Input: !* !* a0,as,ax,ay = coefficients in the function: f = a0 + as*sin(ax*x+ay*y). !* x = x-coordinate at which the function/derivative is evaluated. !* y = y-coordinate at which the function/derivative is evaluated. !* nx = nx-th derivative with respect to x (nx >= 0). !* ny = ny-th derivative with respect to y (ny >= 0). !* !* Output: The function value. !* !* !* Below are some examples: !* !* f = a0 + as*sin(ax*x+ay*y) !<- (nx,ny)=(0,0) !* !* fx = ax * as*cos(ax*x+ay*y) !<- (nx,ny)=(1,0) !* fy = ay * as*cos(ax*x+ay*y) !<- (nx,ny)=(0,1) !* !* fxx = -ax**2 * as*sin(ax*x+ay*y) !<- (nx,ny)=(2,0) !* fxy = -ax*ay * as*sin(ax*x+ay*y) !<- (nx,ny)=(1,1) !* fyy = -ay**2 * as*sin(ax*x+ay*y) !<- (nx,ny)=(0,2) !* !* fxxx = -ax**3 * as*cos(ax*x+ay*y) !<- (nx,ny)=(3,0) !* fxxy = -ax**2 *ay * as*cos(ax*x+ay*y) !<- (nx,ny)=(2,1) !* fxyy = -ax *ay**2 * as*cos(ax*x+ay*y) !<- (nx,ny)=(1,2) !* fyyy = - ay**3 * as*cos(ax*x+ay*y) !<- (nx,ny)=(0,3) !* !* fxxxx = ax**4 * as*sin(ax*x+ay*y) !<- (nx,ny)=(4,0) !* fxxxy = ax**3 *ay * as*sin(ax*x+ay*y) !<- (nx,ny)=(3,1) !* fxxyy = ax**2 *ay**2 * as*sin(ax*x+ay*y) !<- (nx,ny)=(2,2) !* fxyyy = ax *ay**3 * as*sin(ax*x+ay*y) !<- (nx,ny)=(1,3) !* fyyyy = ay**4 * as*sin(ax*x+ay*y) !<- (nx,ny)=(0,4) !* !* and so on. !* !* !******************************************************************************** function manufactured_sol(a0,as,ax,ay, nx,ny,x,y) result(fval) implicit none integer , parameter :: p2 = selected_real_kind(p=15) !Input real(p2), intent(in) :: a0, as, ax, ay, x, y integer , intent(in) :: nx, ny !Output real(p2) :: fval if (nx < 0 .or. ny < 0) then write(*,*) " Invalid input: nx and ny must be greater or equal to zero... Try again." stop endif if ( nx+ny == 0 ) then fval = a0 + as*sin(ax*x + ay*y) elseif ( mod(nx+ny,2) == 0 ) then fval = - (ax**nx * ay**ny)*as*sin(ax*x + ay*y) if ( mod(nx+ny,4) == 0 ) fval = -fval else fval = (ax**nx * ay**ny)*as*cos(ax*x + ay*y) if ( mod(nx+ny+1,4) == 0 ) fval = -fval endif end function manufactured_sol !******************************************************************************** !******************************************************************************** !* !* This function computes viscosity by Sutherland's law !* !* ------------------------------------------------------------------------------ !* Input: T = temperature !* !* Output: mu = viscosity computed by Sutherland's law !* ------------------------------------------------------------------------------ !* !* Note: The viscosity is scaled by the ratio of the freestream Mach number !* to the freestream Reynolds number: M_inf/Re_inf. This corresponds to !* the non-dimensionalized Navier-Stokes equations with the freestream !* speed of sound: (u,v)=(u/a_inf,v/a_inf), p=p/(rho_inf*a_inf^2). !* Modify the factor if other (or no) non-dimensionalization is used. !* See "I do like CFD, VOL.1" for details, which is available in PDF !* at http://www.cfdbooks.com. !* !******************************************************************************** function viscosity(T,T_inf,M_inf,Re_inf) result(mu) implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2), parameter :: one = 1.0_p2 real(p2), parameter :: three = 3.0_p2 real(p2), parameter :: half = 0.5_p2 !Input real(p2), intent(in) :: T, T_inf, M_inf, Re_inf !Output real(p2) :: mu !Parameter for Sutherland's law: real(p2), parameter :: C = 110.5_p2 ! Sutherland's law (nondimensional form) mu = (one+C/T_inf)/(T+C/T_inf)*T**(three*half) ! Scaling for the non-dimensionalized NS equations. mu = mu * M_inf/Re_inf end function viscosity !******************************************************************************** !******************************************************************************** !* !* This function computes the derivative of Sutherland's law w.r.t. temperature. !* !* ------------------------------------------------------------------------------ !* Input: T = temperature !* !* Output: dmu_dT = derivative of Sutherland's law w.r.t. temperature !* ------------------------------------------------------------------------------ !* !* Note: The viscosity is scaled by the ratio of the freestream Mach number !* to the freestream Reynolds number: M_inf/Re_inf. This corresponds to !* the non-dimensionalized Navier-Stokes equations with the freestream !* speed of sound: (u,v)=(u/a_inf,v/a_inf), p=p/(rho_inf*a_inf^2). !* Modify the factor if other (or no) non-dimensionalization is used. !* See "I do like CFD, VOL.1" for details, which is available in PDF !* at http://www.cfdbooks.com. !* !******************************************************************************** function dviscosity_dT(T,T_inf,M_inf,Re_inf) result(dmu_dT) implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2), parameter :: one = 1.0_p2 real(p2), parameter :: three = 3.0_p2 real(p2), parameter :: half = 0.5_p2 !Input real(p2), intent(in) :: T, T_inf, M_inf, Re_inf !Output real(p2) :: dmu_dT !Parameter for Sutherland's law: real(p2), parameter :: C = 110.5_p2 real(p2) :: mu ! Sutherland's law scaled for the non-dimensionalized NS equations by M_inf/Re_inf. mu = (one+C/T_inf)/(T+C/T_inf)*T**(three*half) * M_inf/Re_inf ! Derivative with respect to the temperature. dmu_dT = half*mu/(T+C/T_inf)*(one + three*(C/T_inf)/T) end function dviscosity_dT !******************************************************************************** end program edu2d_mms_twodns