Fortran CFD 代码中出现非数字

Question

我正在尝试求解通过双收敛-发散喷嘴的流体流动的欧拉一维方程。

这是我的完整代码：

module constants

double precision ::GAM, GAMM1 ,rkconst(4),conA,conB,conC,conK,R, T0
double precision, allocatable :: xc(:),xf(:),Sc(:),Sf(:), Scp(:)

end module constants

Program Euler1d
  
  use constants

  Implicit None

  double precision, allocatable :: U(:,:),U0(:,:),cell_flux(:,:),P(:)

  double precision :: XMIN, XMAX, DX, DT,TIME

  double precision :: UCL(3),UCR(3),UPL(3),UPR(3), NORM,Flux(3)

  integer          :: i, j, k, NX,NXC,NTIME, NTIMETOT

R = 287.05;

T0 = 276;

GAM=1.4; GAMM1=GAM - 1

conA= 0.408

conB= 0.93

conC= 1.07

conK= 0.17

!Define runge kutta constants

  rkconst(1) = 0.1084d0

  rkconst(2) = 0.2602d0

  rkconst(3) = 0.5052d0

  rkconst(4) = 1.d0

!Input Parameters

  NTIMETOT = 5000            ! Total timesteps

  XMIN = 0.0                 ! XMIN

  XMAX = 2.0                 ! XMAX

  DX = 1.0E-2                ! DX

  DT = 1.0E-6                ! DT

!Finite Volume Scheme-Nodes represent faces

!Number of points/faces

  NX = int((XMAX-XMIN)/DX)

!Calculate new DX for accurate division by NX

  DX = (XMAX-XMIN)/(NX-1)

  !write (*,*) 'New DX =', DX

!Number of cell centers

  NXC= NX -1

! 0: represents boundary cell at start and NXC+1 boundary cell at end

! U is the solution vector (Dens, Dens *U, Dens *E), U0 the vector after the 1st Runge-Kutta step

! cell_flux is the sum of the fluxes for each cell

  allocate (U(3,0:NXC+1),U0(3,0:NXC+1))

  allocate (cell_flux(1:3,0:NXC+1))
  
  allocate (P(1:NXC))
 
allocate(xf(1:NX))
allocate(xc(1:NXC+1))
allocate(Sf(1:NX))
allocate(Sc(0:NXC+1))
allocate(Scp(1:NXC+1))

do i=1,NX
xf(i)= XMIN + (i-1)*DX
xc(i)= xf(i) + DX/2
Sf(i)= conK+conA*((xf(i)-conC)**2)*((xf(i)-conC)**2-conB**2)
Sc(i)= conK+conA*((xc(i)-conC)**2)*((xc(i)-conC)**2-conB**2)
Scp(i)=4*conA*(xc(i)-conC)**3-(2*conA*conB**2)*(xc(i)-conC)

  open (2,file='areas.dat')

  !write(2,'(a)') '#index, x of face, x of center, face, center, prime'
  
  write(2,'(6(e28.17,1x))') i, xf(i), xc(i), Sf(i), Sc(i), Scp(i)

!write (*,*) i, Sc(i)

enddo
Sc(0)= Sc(1)

write (2,*)

TIME=0

!Set initial conditions

  call ic
  
   !write (*,*) 'U initial is', U(1:3,0:NXC+1)

  do NTIME = 1, NTIMETOT

     k=1

!k 1:4 four step RK

     do while (k.le.4) 

!Set boundary conditions

        !write (*,*) NTIME, U(1,1)
!=========================== SOMETHING happens as soon as NTIME = 56, and from k = 1 to k = 2. molis k = 2 xtypaei NaN, sto k = 1 den xtypaei ==============================================================================================================================================
!if (NTIME .eq. 56 .and. k .eq. 2) then
   !write (*,*) NTIME, U(1:3,0:NXC+1)
!endif
        call bc

   !write (*,*) 'U is', U(1:3,0:NXC+1)

        cell_flux=0.d0

!For all Faces

        do i=1,NX
        
  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) i, U(1:3,i)
  !endif
  
           !Left cell

           UCL = U(1:3,i-1)  
                 
  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) i, UCL
  !endif
  
           !Right cell

           UCR = U(1:3,i)        

           !Get primitive variables
           
  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) 'UCL(1) is', UCL(1)
  !endif

           call co_to_pri(UCL,UPL,Sf(i))

           call co_to_pri(UCR,UPR,Sf(i))

           NORM=1

           !calculate Roe-flux

           call calc_flux(UPL, UPR, UCL, UCR,NORM,Flux)

           !Normal points from left to right

           !For bc's it's obsolete

           cell_flux(1:3,i-1) = cell_flux(1:3,i-1)  + Flux  

           cell_flux(1:3,i)   = cell_flux(1:3,i)    - Flux

        enddo

!    Time integration

       if(k.eq.1) U0(1:3,1:NXC)=U(1:3,1:NXC)
       do i=1,NXC
        P(i) = GAMM1 * (U(3,i) - 0.5d0 / U(1,i) * U(2,i)**2/Sf(i))
  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) 'P vector is ', P(i)
  !endif
  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) 'U matrix is ', U(1:3,i)
  !endif
  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) i, U(1,i)
  !endif
        U(1:3,i) = U0(1:3,i) - DT/DX *rkconst(k)*(cell_flux(1:3,i) - cell_flux(1:3,i-1))
        U(2,i) = U(2,i) + (DT*P(i))*Scp(i)
       enddo
        k=k+1

     enddo
!if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.540E-4) then
   !write (*,*) 'Umatrix is ', U(1:3,1:NX)
!endif
     TIME = TIME + DT

     !write(*,*) 'NTIME= ',NTIME,'TIME= ', TIME

!write results 

     if (mod(NTIME,100).eq.0) call writeres(TIME)

  enddo

contains

!Set initial conditions

Subroutine ic

 Implicit None

 double precision :: UC(3),UP(3)

!Sod's shock tube initial conditions

! set pressure, velocity, density and then transform to conservative
       
 UP(3)= 3.3E5; UP(2) = 0.d0; UP(1) = UP(3)/(R*T0);

 do i=1,NXC

       call pri_to_co(UC,UP,Sf(i))

       U(1:3,i) = UC(1:3)

   !write (*,*) 'initial matrix is ', U(1:3,i)
   !write (*,*) NTIME, U(1,1)

 enddo

End Subroutine ic

!Set boundary conditions at ends

Subroutine bc

double precision :: UP(3),UC(3)

!give bc in a primitive form 

!start

 UP(1) =   4E5/(R*T0) !Density

 UP(2) =  U(2,1)/U(1,1) !Velocity

 UP(3) =   4E5   !Pressure
 
 !write (*,*) 'UP is', UP(1:3)
 
  !if (UP(2) .gt. 374) then
   !write (*,*) NTIME, UP(1:3)
  !endif

 call pri_to_co(UC,UP,Sc(0))

 U(1,0) =  UC(1) !Density

 U(2,0) =  UC(2) !Velocity

 U(3,0) =  UC(3) !Pressure

!end

 UP(1) =  U(1,NXC)/(conK+conA*((2-conC)**2)*((2-conC)**2-conB**2)) !Density

 UP(2) =  U(2,NXC)/U(1,NXC) !Velocity

 UP(3) =  3.2e5     !Pressure

 call pri_to_co(UC,UP,Sc(NXC+1))

 U(1,NXC+1) =  UC(1) !Density

 U(2,NXC+1) =  UC(2) !Velocity

 U(3,NXC+1) =  UC(3) !Pressure

end Subroutine bc

!write results in primitive variable form

Subroutine writeres(TIME)

double precision :: TIME,UP(3),UC(3)

  open (1,file='results.dat')

  write(1,'(a)') '#Time, X, Density, Velocity, Pressure'

  do i=1,NXC

     UC(1:3) = U(1:3,i)

     call co_to_pri(UC,UP,Sc(i))

     write(1,'(5(e28.17,1x))') TIME,XMIN + (i-1) * DX + 0.5 * DX, UP(1:3)

  enddo

  write(1,*) 

End Subroutine writeres

! Flux:Calculated Roe flux

Subroutine calc_flux(UPL, UPR, UCL, UCR,NORM,Flux)
use constants

  Implicit None

  double precision, intent(in)  :: UPL(3), UPR(3),UCL(3),UCR(3),NORM

  double precision, intent(out) :: Flux(3)

  double precision              :: Flux_left(3), Flux_right(3), Aroe(3,3), DF1(3)

!1D fluxes

  call flux1d(UPL,UCL,NORM,Flux_left,Sf(i))

  call flux1d(UPR,UCR,NORM,Flux_right,Sf(i))

!Roe matrix

  !if (TIME .gt. 0.539E-4 .and. TIME .lt. 0.54E-4) then
   !write (*,*) 'UPL(1) matrix is ', UPL(1)
  !endif

  call roematrix(UPL,UCL,UPR,UCR,NORM,Aroe,Sf(i))

  DF1(1:3) = UCR(1:3)-UCL(1:3)

!Final flux 

 Flux(1 : 3) = 0.5d0 * (Flux_left(1:3) + Flux_right(1:3) - matmul(Aroe,DF1))

End subroutine calc_flux

End Program Euler1d

!==============================================================================

!Subroutine to  convert conservative (UC) to primitive (UP)

Subroutine co_to_pri(UC,UP,S)

  use ieee_arithmetic, only: ieee_is_nan

  use constants

  Implicit None

  double precision, intent(in)  :: UC(3),S

  double precision, intent(out) :: UP(3)

  UP(1) = UC(1)/S

  UP(2) = UC(2)/UC(1)

  UP(3) = GAMM1 * (UC(3) - 0.5d0 / UC(1) * UC(2)**2)/S
  
    !write(*,*) 'error is',UC(1)

End Subroutine co_to_pri

!Subroutine to  convert primitive(UP) to conservative (UC)

Subroutine pri_to_co(UC,UP,S)

  use constants

  Implicit None

  double precision, intent(out)  :: UC(3)

  double precision, intent(in)   :: UP(3),S

  UC(1) = UP(1)*S

  UC(2) = UP(2)*UP(1)*S

  UC(3) = (UP(3) / GAMM1 + 0.5d0 * UP(1) * UP(2)**2)*S

End Subroutine pri_to_co

!Subroutine to calculate Roe flux

! UPL :Primitive Variable left 

! UPR :Primitive Variable right

! UCL :Conservative Variable left

! UCR :Conservative Variable Right

! NORM:Normal

!1D euler equaations flux functions

! UP :Primitive Variable

! UC :Conservative Variable 

! NORM:Normal

Subroutine flux1d(UP,UC,NORM,Flux,S)

  Implicit None

  double precision, intent(in)  :: UP(3), UC(3),NORM, S

  double precision, intent(out) :: Flux(3)

  double precision :: VN

  VN =  UP(2) * NORM

  Flux(1) = UP(1) * VN * S

  Flux(2) = (UP(1) * UP(2)* VN + UP(3) * NORM) * S

  Flux(3) = (UC(3) + UP(3)) * VN * S

End subroutine flux1d

!Subroutine to calculate Roe matrix

! UPL :Primitive Variable left 

! UPR :Primitive Variable right

! UCL :Conservative Variable left

! UCR :Conservative Variable Right

! NORM:Normal

! Aroe: ROe matrix in conservative variables

Subroutine roematrix(UPL,UCL,UPR,UCR,NORM,Aroe,S)

  use ieee_arithmetic, only: ieee_is_nan
  
  use constants

  Implicit None

  double precision, intent(in)  :: UPL(3),UCL(3),UPR(3),UCR(3),NORM, S

  double precision, intent(out) :: Aroe(3,3)

  double precision              :: DensRoe, Uroe, Hl,Hr, Hroe,croe,L1,L2,L3,delta

  double precision              :: left(3,3),right(3,3),eigen(3,3)

!Calculate Roe averaged variables (Density,Velocity,Enthalpy)

  DensRoe = sqrt(UPL(1))*sqrt(UPR(1))
  
  Uroe    = (sqrt(UPL(1)) * UPL(2) + sqrt(UPR(1)) * UPR(2)) / ( sqrt(UPL(1)) + sqrt(UPR(1)) )

  Hl     = (UCL(3)/S + UPL(3))/UPL(1)

  Hr     = (UCR(3)/S + UPR(3))/UPR(1)

  Hroe    = (sqrt(UPL(1)) * Hl + sqrt(UPR(1)) * Hr) / ( sqrt(UPL(1)) + sqrt(UPR(1)) )

  croe    =  sqrt( GAMM1 * (Hroe - 0.5d0 * Uroe**2) )
  
  !if (ieee_is_nan(croe)) then
  
     !write (*,*) UPL(1), UPL(2)
     
  !endif

!Get left and right eigenvectors in conservative form

  call leftrightcons(DensRoe,croe,Uroe,left,right)

  eigen     = 0.d0

   L1=abs(Uroe)

   L2=abs(Uroe+croe)

   L3=abs(Uroe-croe)

!Harten's entropy fix  and fill eigenvalue matrix

   delta = 0.05*croe

   if (L1.le.delta) L1 = (L1**2 + delta**2) / (2 * delta)

   if (L2.le.delta) L2 = (L2**2 + delta**2) / (2 * delta)

   if (L3.le.delta) L3 = (L3**2 + delta**2) / (2 * delta)

   eigen=0

   eigen(1,1)= L1

   eigen(2,2)= L2

   eigen(3,3)= L3

!get Roe matrix

   Aroe = matmul(right,eigen)

   Aroe = matmul(Aroe,left)

End Subroutine roematrix

!Subroutine to calculate left and right eigenvectors in conservative form

Subroutine leftrightcons(Dens,C,U,left,right)

   use constants

   Implicit None

   double precision, intent(in) :: Dens,U,C

   double precision, intent(out):: left(3,3),right(3,3)

   left(1,1) =  1.d0 - 0.5d0 * GAMM1 * U**2/C**2

   left(1,2) =  GAMM1 * U /C**2

   left(1,3) = -GAMM1/C**2

   left(2,1) = (0.5d0*GAMM1*U**2 - U*C) * (1.d0/(Dens*C))

   left(2,2) = (C - GAMM1*U) *(1.d0/(Dens*C))

   left(2,3) = GAMM1/(Dens*C)

   left(3,1) = -(0.5d0*GAMM1*U**2 + U*C) * (1.d0/(Dens*C))

   left(3,2) = (C + GAMM1*U) *(1.d0/(Dens*C))

   left(3,3) = -GAMM1/(Dens*C)

   right(1,1) =  1.d0

   right(1,2) =  0.5d0 * Dens/C

   right(1,3) = -0.5d0 * Dens/C
   
   right(2,1) = U

   right(2,2) = 0.5d0 * (U+C) * Dens/C

   right(2,3) =-0.5d0 * (U-C) * Dens/C

   right(3,1) = 0.5d0 * U**2

   right(3,2) = (0.5d0 * U**2 + U*C + C**2/GAMM1) * 0.5d0 * Dens/C

   right(3,3) =-(0.5d0 * U**2 - U*C + C**2/GAMM1) * 0.5d0 * Dens/C

End Subroutine leftrightcons

!==============================================================================

这段代码应该可以正常工作，并用密度、速度和压力的预期结果写入文件 results.dat。但是，在某些时候，某处出现 NaN 值，它会影响所有后续计算，并且 results.dat 文件最终会崩溃，充满 NaN 值。通过一些搜索，NaN 的第一次出现可能是在 NTIME = 53, k = 2 迭代或 NTIME = 54, k = 2 附近。另外，其他人告诉我，在某些时候平方根的参数在croe的计算中（第486行）是负数。我对 Fortran 和调试的了解非常有限，而且我现在没有时间学习，因为我很快就会报告这一点。如果有人有任何想法，我们将不胜感激。

Answer 1

您的代码有一些重大错误。

最大的问题是在时间积分中减去通量

cell_flux(1:3,i) - cell_flux(1:3,i-1)

。数组 cell_flux 已保存该单元格的通量散度：您不需要再次减去。

P(i) 和 Flux(3) 的表达式也有错误。

恐怕我必须重构你的代码才能阅读它。这是一个可以运行的版本。请注意，您必须离开并将 NTIMETOT 增加至少 10 倍才能接近稳定解 - 这将大大增加运行时间，但会显示出口附近形成冲击。

我没有改变你的时间步长例程（因为我不认识它），也没有改变你的 Roe 通量（因为我处理黎曼求解器已经太久了，而且我也忘记了 Phi Roe 的论文）。

module constants

   double precision, parameter :: GAM = 1.4, GAMM1 = GAM - 1         ! ratio of specific heats
   double precision, parameter :: R = 287.05                         ! ideal gas constant for air
   double precision, parameter :: T0 = 276.0                         ! static temperature (K)
   double precision, parameter :: rkconst(4) = [ 0.1084d0, 0.2602d0, 0.5052d0, 1.0d0 ]

   double precision, allocatable :: xc(:), Sc(:)
   double precision, allocatable :: xf(:), Sf(:)

end module constants

!===============================================================================

program Euler1d
   use constants
   implicit none
   double precision, external :: area
   double precision, allocatable :: U(:,:), U0(:,:), cell_flux(:,:), P(:)
   double precision XMIN, XMAX, DX, DT,TIME
   double precision UCL(3), UCR(3), UPL(3), UPR(3), Flux(3)
   double precision NORM
   integer NX, NXC, NTIME, NTIMETOT
   integer i, k

   ! Input Parameters
   NTIMETOT = 5000            ! Total timesteps (needs to be raised to at least 10 times this)
   DT = 1.0E-6
   XMIN = 0.0;   XMAX = 2.0
   DX = 1.0E-2             

   ! Number of nodes and faces
   NX = int( ( XMAX - XMIN ) / DX )         ! Number of faces
   DX = ( XMAX - XMIN ) / ( NX - 1 )        ! Calculate new DX for accurate division by NX
   NXC = NX - 1                             ! Number of interior cell centers

   ! 0: represents ghost cell at inflow and NXC+1 ghost cell at end
   ! U is the solution vector [density * area, density * velocity * area, density * energy * area ],
   ! U0 is the vector after the 1st Runge-Kutta step
   ! Cell_flux is the NET FLUX OUT OF A CELL
   allocate( U(3,0:NXC+1), U0(3,0:NXC+1), P(0:NXC+1), xc(0:NXC+1), Sc(0:NXC+1) )
   allocate( cell_flux(1:3,0:NXC+1) )
   allocate( xf(1:NX), Sf(1:NX) )

   ! Set locations and areas
   do i = 1, NX
      xf(i) = XMIN + ( i - 1 ) * DX
      xc(i) = xf(i) + DX / 2
      Sf(i) = area( xf(i) )
      Sc(i) = area( xc(i) )
   end do
   xc(0   )  = xc(1 ) - DX
   xc(NX+1)  = xc(NX) + DX
   Sc(0)     = Sc(1)
   Sc(NXC+1) = Sc(NXC)

   ! Set initial conditions
   TIME = 0
   call ic

   ! Timestepping
   do NTIME = 1, NTIMETOT

      do k = 1, 4                      ! four-step RK
         ! Set boundary conditions
         call bc

         ! Calculate net fluxes from cells 1 to NXC (=NX-1)
         cell_flux = 0

         do i = 1, NX
            ! Left/right cells
            UCL = U(:,i-1)
            UCR = U(:,i)

            ! Get primitive variables
            call co_to_pri( UCL, UPL, Sc(i-1) )
            call co_to_pri( UCR, UPR, Sc(i  ) )

            ! Normal points from left to right
            NORM = 1

            ! Calculate flux on cell face
            call calc_flux( UPL, UPR, UCL, UCR, NORM, Flux )

            ! Compound the fluxes OUT of adjacent cells
            cell_flux(:,i-1) = cell_flux(1:3,i-1) + Flux
            cell_flux(:,i  ) = cell_flux(1:3,i  ) - Flux

         end do

         ! Time integration
         if( k == 1 ) U0 = U
         do i = 1, NXC
            P(i) = GAMM1 * ( U(3,i) - 0.5 * U(2,i) ** 2  / U(1,i) ) / Sc(i)
            U(:,i) = U0(:,i) - DT / DX * rkconst(k) * cell_flux(:,i)           ! d(U)/dt = -d(F)/dx
            U(2,i) = U(2,i) + DT * P(i) * ( Sf(i+1) - Sf(i) ) / DX             ! effect of side-wall expansion
         end do
      end do
      TIME = TIME + DT

      if ( mod( NTIME, 100 ) == 0 ) call writeres( TIME )

  end do

contains

   ! Set initial conditions
   subroutine ic
      implicit None
      double precision :: UP(3)
      integer i

      ! set pressure, velocity, density and then transform to conservative
      UP(3) = 3.3E5;  UP(2) = 0;   UP(1) = UP(3) / ( R * T0 );

      do i = 0, NXC + 1
         call pri_to_co( U(:,i), UP, Sc(i) )
      enddo

   end subroutine ic

   !----------------------------------------------

   subroutine bc
      double precision UP(3)

      ! Inlet
      UP(3) = 4E5                           ! pressure
      UP(1) = UP(3) / ( R * T0 )            ! density (based on pressure and temperature)
      UP(2) = U(2,1) / U(1,1)               ! velocity (extrapolated from interior)
      call pri_to_co( U(:,0), UP, Sc(0) )
      call pri_to_co( U(:,0), UP, Sc(1) )

      ! Outlet
      UP(3) = 3.2e5                         ! pressure
      UP(1) = U(1,NXC) / Sc(NXC+1)          ! density (extrapolated from interior)
      UP(2) = U(2,NXC) / U(1,NXC)           ! velocity (extrapolated from interior)
      call pri_to_co( U(:,NXC+1), UP, Sc(NXC+1) )
      call pri_to_co( U(:,NXC+1), UP, Sc(NXC  ) )

   end subroutine bc

   !----------------------------------------------

   subroutine writeres( TIME )
      double precision, intent(in) :: TIME
      double precision UP(3)
      integer i

      open( 10, file="results.dat" )
!     write( 10, "(a)" ) "#Time, X, Density, Velocity, Pressure"
      do i = 0, NXC + 1
         call co_to_pri( U(:,i), UP, Sc(i) )
         write( 10, "( 5( e13.6, 1x ) )" ) TIME, xc(i), UP
      end do
      close( 10 )

   end subroutine writeres

   !----------------------------------------------

   ! Calculated flux
   subroutine calc_flux( UPL, UPR, UCL, UCR, NORM, Flux )
      use constants
      implicit None
      double precision, intent(in)  :: UPL(3), UPR(3), UCL(3), UCR(3), NORM
      double precision, intent(out) :: Flux(3)
      double precision Flux_left(3), Flux_right(3), ARoe(3,3)

      ! 1D fluxes
      call flux1d( UPL, UCL, NORM, Flux_left , Sc(i) )
      call flux1d( UPR, UCR, NORM, Flux_right, Sc(i) )

      ! Roe matrix
      call Roematrix( UPL, UCL, UPR, UCR, ARoe, Sf(i) )
      Flux = 0.5 * ( Flux_left + Flux_right - matmul( ARoe, UCR - UCL ) )

   end subroutine calc_flux

end program Euler1d

!===============================================================================

! Function returning cross-sectional area
double precision function area( x )
   double precision, intent(in) :: x
   double precision, parameter :: conA = 0.408, conB= 0.93, conC= 1.07, conK= 0.17

   area = conK + conA * ( ( x - conC ) ** 2 ) * ( ( x - conC ) ** 2 - conB **2 )

end function area

!-------------------------------------------------------------------------------

! Convert conservative (UC) to primitive (UP) variables
subroutine co_to_pri( UC, UP, S )
  use constants
  implicit none
  double precision, intent(in ) :: UC(3), S
  double precision, intent(out) :: UP(3)

  UP(1) = UC(1) / S                                             ! density
  UP(2) = UC(2) / UC(1)                                         ! velocity
  UP(3) = GAMM1 * ( UC(3) - 0.5 * UC(2) ** 2 / UC(1) ) / S      ! pressure

end subroutine co_to_pri

!-------------------------------------------------------------------------------

! Convert primitive(UP) to conservative (UC)
subroutine pri_to_co( UC, UP, S )
  use constants
  implicit none
  double precision, intent(out) :: UC(3)
  double precision, intent(in ) :: UP(3), S

  UC(1) = UP(1) * S                                             ! density * area
  UC(2) = UP(2) * UP(1) * S                                     ! density * velocity * area
  UC(3) = ( UP(3) / GAMM1 + 0.5 * UP(1) * UP(2) ** 2 ) * S      ! volumetric internal energy * area

end subroutine pri_to_co

!-------------------------------------------------------------------------------

! Calculate flux vector
subroutine flux1d( UP, UC, NORM, Flux, S )
  implicit none
  double precision, intent(in ) :: UP(3), UC(3), NORM, S
  double precision, intent(out) :: Flux(3)
  double precision VN

  VN =  UP(2) * NORM
  Flux(1) = UP(1) * VN * S                    
  Flux(2) = ( UP(1) * UP(2) * VN + UP(3) ) * S
  Flux(3) = ( UC(3) + UP(3) * S ) * VN

end subroutine flux1d

!-------------------------------------------------------------------------------

! Calculate Roe matrix
subroutine Roematrix( UPL, UCL, UPR, UCR, Aroe, S)
   use constants
   implicit none
   double precision, intent(in ) :: UPL(3), UCL(3), UPR(3), UCR(3), S
   double precision, intent(out) :: ARoe(3,3)
   double precision densRoe, uRoe, hl, hr, hRoe, cRoe, L1, L2, L3, delta
   double precision left(3,3), right(3,3), eigen(3,3)
   double precision fl, fr

   ! Calculate Roe-averaged variables
   fl = sqrt( UPL(1) );   fr = sqrt( UPR(1) )
   densRoe = fl * fr
   uRoe    = ( fl * UPL(2) + fr * UPR(2) ) / ( fl + fr )
   hl      = ( UCL(3) / S + UPL(3) ) / UPL(1)
   hr      = ( UCR(3) / S + UPR(3) ) / UPR(1)
   hRoe    = ( fl * hl + fr * hr ) / ( fl + fr )
   cRoe    =  sqrt( GAMM1 * ( hRoe - 0.5 * uRoe **2 ) )

   ! Get left and right eigenvectors in conservative form
   call leftrightcons( densRoe, cRoe, uRoe, left, right )
   eigen = 0
   L1 = abs( uRoe        )
   L2 = abs( uRoe + cRoe )
   L3 = abs( Uroe - cRoe )
   ! Harten's entropy fix  and fill eigenvalue matrix
   delta = 0.05 * cRoe
   if ( L1 <= delta ) L1 = ( L1 ** 2 + delta ** 2 ) / ( 2 * delta )
   if ( L2 <= delta ) L2 = ( L2 ** 2 + delta ** 2 ) / ( 2 * delta )
   if ( L3 <= delta ) L3 = ( L3 ** 2 + delta ** 2 ) / ( 2 * delta )
   eigen = 0
   eigen(1,1) = L1
   eigen(2,2) = L2
   eigen(3,3) = L3
   ARoe = matmul( right, eigen )
   ARoe = matmul( ARoe , left  )

end subroutine Roematrix

!-------------------------------------------------------------------------------

! Calculate left and right eigenvectors in conservative form
subroutine leftrightcons( dens, c, u, left, right )
   use constants
   implicit none
   double precision, intent(in) :: dens, u, c
   double precision, intent(out):: left(3,3), right(3,3)

   left(1,1)  =  1 - 0.5 * GAMM1 * u ** 2 / c ** 2
   left(1,2)  =  GAMM1 * u / c ** 2
   left(1,3)  = -GAMM1 / c ** 2
   left(2,1)  = ( 0.5 * GAMM1 * u ** 2 - u * c ) / ( dens * c )
   left(2,2)  = ( c - GAMM1 * u ) / ( dens * c )
   left(2,3)  = GAMM1 / ( dens * c )
   left(3,1)  = -( 0.5 * GAMM1 * u ** 2 + u * c ) / ( dens * c )
   left(3,2)  = ( c + GAMM1 * u ) / ( dens * c )
   left(3,3)  = -GAMM1 / ( dens * c )
   right(1,1) =  1
   right(1,2) =  0.5 * dens / c
   right(1,3) = -0.5 * dens / c
   right(2,1) = u
   right(2,2) = 0.5 * ( u + c ) * dens / c
   right(2,3) =-0.5 * ( u - c ) * dens / c
   right(3,1) = 0.5 * u ** 2
   right(3,2) = ( 0.5 * u ** 2 + u * c + c ** 2 / GAMM1 ) * 0.5 * dens / c
   right(3,3) =-( 0.5 * u ** 2 - u * c + c ** 2 / GAMM1 ) * 0.5 * dens / c

end subroutine leftrightcons

!==============================================================================

当我将 NTIMETOT 增加到 500000（记得打开优化！）时，我得到以下压力图。出口前似乎发生了震动。这在哪里以及有多大取决于您的入口和出口压力。

Fortran CFD 代码中出现非数字

问题描述投票：0回答：1

1个回答

最新问题

Fortran CFD 代码中出现非数字

问题描述 投票：0回答：1

1个回答

最新问题

问题描述投票：0回答：1