CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

2D AUSM program (FORTRAN): problems with the pressure and vector implementations

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   January 7, 2016, 03:19
Default 2D AUSM program (FORTRAN): problems with the pressure and vector implementations
  #1
Member
 
garry
Join Date: Oct 2013
Posts: 41
Rep Power: 5
yzeyue is on a distinguished road
Hi everyone,

I program a 2D AUSM code. But I do not know how to deal with the pressure part and the use of vector.

There is pressure part on the energy equation in some codes I see on the website, but some there is not. I see the paper about AUSM scheme, and I find that the pressure part is only in one motion equation. Should the pressure part in all the motion equation?

Appreciate it very much if anyone can help me with the problem.

My AUSM code is as follows.

Code:
! i direction
            nx = XLNI(I,J) / SLNI(I,J)
        ny = YLNI(I,J) / SLNI(I,J)

!  Left state
          rhoL = DENSL
            vL = VELXL  ! velocity of x direction
           vyL = VELYL  ! velocity of y direction
           qnL = vL*nx + vyL*ny
            pL = PRESL
            aL = sqrt(gamma*pL/rhoL)
            mL = qnL / aL
            eL = pL/(gamma-one)+half*(vL*vL+vyL+vyL)
            hL = eL + pL/rhoL
!  Left state
          rhoR = DENSR
            vR = VELXR  ! velocity of x direction
           vyR = VELYR  ! velocity of y direction
           qnR = vR*nx + vyR*ny
            pR = PRESR
            aR = sqrt(gamma*pR/rhoR)
            mR = qnR / aR
            eR = pR/(gamma-one)+half*(vR*vR+vyR+vyR)
            hR = eR + pR/rhoR


           if(ABS(mL).le.one)then
                mLP = quarter*(mL + one)**two
                pLP = quarter*pL*(mL+one)**two*(two-mL)
           else
                mLP = half*(mL + ABS(mL))
                pLP = half*pL*(mL+ABS(mL))/mL
           end if

           if(ABS(mR).le.one)then
                mRN = -quarter*(mR - one)**two
                pRN = quarter*pR*(mR-one)**two*(two+mR)
           else
                mRN = half*(mR - ABS(mR))
                pRN = half*pR*(mR-ABS(mR))/mR
           end if

         m_Mid = mLP+mRN
         

   FLUXI(1,I,J) = half*m_Mid*(rhoL*aL + rhoR*aR)-half*ABS(m_Mid)*(rhoR*aR - rhoL*aL)

   FLUXI(2,I,J) = half*m_Mid*(rhoL*aL*vL + rhoR*aR*vR) - half*ABS(m_Mid)*(rhoR*aR*vR - rhoL*aL*vL)+ (pLP + pRN)*nx
                  
   FLUXI(3,I,J) = half*m_Mid*(rhoL*aL*vyL + rhoR*aR*vyR) - half*ABS(m_Mid)*(rhoR*aR*vyR - rhoL*aL*vyL)+ (pLP + pRN)*ny

   FLUXI(4,I,J) = half*m_Mid*(rhoL*aL*hL + rhoR*aR*hR) - half*ABS(m_Mid)*(rhoR*aR*hR - rhoL*aL*hL)

    

! j direction
            nx = XLNJ(I,J) / SLNJ(I,J)
            ny = YLNJ(I,J) / SLNJ(I,J)

!  Left state
          rhoL = DENSL
            vL = VELXL  ! velocity of x direction
           vyL = VELYL  ! velocity of y direction
           qnL = vL*nx + vyL*ny
            pL = PRESL
            aL = sqrt(gamma*pL/rhoL)
            mL = qnL / aL
            eL = pL/(gamma-one)+half*(vL*vL+vyL+vyL)
            hL = eL + pL/rhoL
!  Left state
          rhoR = DENSR
            vR = VELXR  ! velocity of x direction
           vyR = VELYR  ! velocity of y direction
           qnR = vR*nx + vyR*ny
            pR = PRESR
            aR = sqrt(gamma*pR/rhoR)
            mR = qnR / aR
            eR = pR/(gamma-one)+half*(vR*vR+vyR+vyR)
            hR = eR + pR/rhoR


           if(ABS(mL).le.one)then
                mLP = quarter*(mL + one)**two
                pLP = quarter*pL*(mL+one)**two*(two-mL)
           else
                mLP = half*(mL + ABS(mL))
                pLP = half*pL*(mL+ABS(mL))/mL
           end if

           if(ABS(mR).le.one)then
                mRN = -quarter*(mR - one)**two
                pRN = quarter*pR*(mR-one)**two*(two+mR)
           else
                mRN = half*(mR - ABS(mR))
                pRN = half*pR*(mR-ABS(mR))/mR
           end if

         m_Mid = mLP+mRN
         

   FLUXJ(1,I,J) = half*m_Mid*(rhoL*aL + rhoR*aR)-half*ABS(m_Mid)&
                  *(rhoR*aR - rhoL*aL)

   FLUXJ(2,I,J) = half*m_Mid*(rhoL*aL*vL + rhoR*aR*vR) - half*ABS(m_Mid)*(rhoR*aR*vR - rhoL*aL*vL)+ (pLP + pRN)*nx

   FLUXJ(3,I,J) = half*m_Mid*(rhoL*aL*vyL + rhoR*aR*vyR) - half*ABS(m_Mid)*(rhoR*aR*vyR - rhoL*aL*vyL)+ (pLP + pRN)*ny
                  
   FLUXJ(4,I,J) = half*m_Mid*(rhoL*aL*hL + rhoR*aR*hR) - half*ABS(m_Mid)*(rhoR*aR*hR - rhoL*aL*hL)

    
12        CONTINUE
      DO 31 J=1,JN-1
      DO 31 I=1,IN-1
        FLUX(1,I,J)=FLUXI(1,I+1,J)*SLNI(I+1,J)-FLUXI(1,I,J)*SLNI(I,J)&
               +FLUXJ(1,I,J+1)*SLNJ(I,J+1)-FLUXJ(1,I,J)*SLNJ(I,J)
        FLUX(2,I,J)=FLUXI(2,I+1,J)*SLNI(I+1,J)-FLUXI(2,I,J)*SLNI(I,J)&
               +FLUXJ(2,I,J+1)*SLNJ(I,J+1)-FLUXJ(2,I,J)*SLNJ(I,J)
        FLUX(3,I,J)=FLUXI(3,I+1,J)*SLNI(I+1,J)-FLUXI(3,I,J)*SLNI(I,J)&
               +FLUXJ(3,I,J+1)*SLNJ(I,J+1)-FLUXJ(3,I,J)*SLNJ(I,J)
        FLUX(4,I,J)=FLUXI(4,I+1,J)*SLNI(I+1,J)-FLUXI(4,I,J)*SLNI(I,J)&
               +FLUXJ(4,I,J+1)*SLNJ(I,J+1)-FLUXJ(4,I,J)*SLNJ(I,J)

31    CONTINUE

Last edited by wyldckat; January 9, 2016 at 11:05. Reason: Added [CODE][/CODE] markers
yzeyue is offline   Reply With Quote

Old   January 7, 2016, 03:45
Default
  #2
Member
 
garry
Join Date: Oct 2013
Posts: 41
Rep Power: 5
yzeyue is on a distinguished road
Hi, how to deal with this problem?
yzeyue is offline   Reply With Quote

Reply

Tags
ausm, cfd, code, discretion

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 03:41.