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

AUSM+, MUSCL coding problem

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

Reply
 
LinkBack Thread Tools Display Modes
Old   September 5, 2013, 04:27
Default AUSM+, MUSCL coding problem
  #1
New Member
 
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 3
Onim is on a distinguished road
Thanks for visiting this thread.

I am working on coding numerical algorithm for 3D Euler equation with general coordinate. Since the flow is compressible, I need to utilize upwinding.

After I implement AUSM+ and MUSCL for upwinding/flux limiting, I tried to solve Sod shock tube problem, and got an errorneous result. Rarefaction wave and shock wave are well refined, but contact discontinuity is oscillating.

If I use 1st order upwind difference for numerical flux, I can get a proper solution. I tried to figure out any bug in MUSCL code, but it was unsuccessful.

In this case, which part of the code should be reviewed?
Attached Images
File Type: jpg sod.jpg (26.4 KB, 7 views)
Onim is offline   Reply With Quote

Old   September 5, 2013, 05:33
Default
  #2
Senior Member
 
duri
Join Date: May 2010
Posts: 130
Rep Power: 7
duri is on a distinguished road
You said it is 1st order scheme but using MUSCL. This is quite confusing. 1st order scheme is monotone preserving no need to use techniques like MUSCL to preserve monotonicity.
Basic first order AUSM+ is sufficient enough to produce 1st order solution for sod problem.
Could you provide the details of face flux estimation from cell averaged data that would help in fixing this issue.
duri is offline   Reply With Quote

Old   September 5, 2013, 05:34
Default
  #3
New Member
 
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 3
Onim is on a distinguished road
I attach some MUSCL code for your information.
Code:
#define indcn i,j,k
#define indup i-ish,j-jsh,k-ksh
#define inddn i+ish,j+jsh,k+ksh

C     For left flux
      DO k=0,nz; DO j=0,ny; DO i=0,nx

C     Obtain gradient ratio r
c
c     Definition of r:
c           r_f= (phi_dnstream - phi_center)
c               /(phi_center - phi_upstream)
      rrho = slopeFrac(rho(indup),rho(indcn),rho(inddn))
      rpx  = slopeFrac(cpx(indup),cpx(indcn),cpx(inddn))

C     Obtain flux limiter function
      psi_rho =fluxlimit(rrho)
      psi_px  =fluxlimit(rpx)
C     Apply flux limiter and obtain final vector
      phi_left(indcn,1)   =rho(indcn)
     +                 + 0.5d0*psi_rho*(rho(indcn)-rho(indup))
      p_left(indcn,2)     =cpx(indcn)
     +                 + 0.5d0*psi_px*(cpx(indcn)-cpx(indup))
c     index 2 means that pressure in xi direction affects
c        quantity rho*u
      ENDDO;ENDDO;ENDDO

#undef indup
#undef indcn
#undef inddn


#define indup i+ish,j+jsh,k+ksh
#define indcn i,j,k
#define inddn i-ish,j-jsh,k-ksh

C     For right flux
      DO k=1,nz+1; DO j=1,ny+1; DO i=1,nx+1

      rrho = slopeFrac(rho(indup),rho(indcn),rho(inddn))


C     Obtain flux limiter function
      psirho =fluxlimit(rrho)

C     Apply flux limiter and obtain final vector
      phi_right(inddn,1)   =rho(indcn)
     +                 + 0.5d0*psirho*(rho(indcn)-rho(indup))

      ENDDO;ENDDO;ENDDO

      FUNCTION slopeFrac(upwind,center,downwind) RESULT(frac)

      nom   = downwind-center
      denom = center-upwind
c     make sure division by 0 does not happen
      if(abs(nom).lt.1d-14)then
        nom = 0.0d0
        denom = 1.0d0
      elseif(nom.gt.1d-14.and.abs(denom).lt.1d-14)then
        nom = 1d14
        denom = 1d0
      elseif(nom.lt.-1d-14.and.abs(denom).lt.1d-14)then
        nom = -1d14
        denom = 1d0
      endif
      frac = nom/denom

      FUNCTION fluxlimit(var) RESULT(fl)

      fl=MAX(0.0d0,MIN(2.0d0*var,0.5d0*(var+1.0d0),2.0d0))

Last edited by Onim; September 5, 2013 at 06:09. Reason: comment was not up to date
Onim is offline   Reply With Quote

Old   September 5, 2013, 05:38
Default
  #4
New Member
 
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 3
Onim is on a distinguished road
Thank your for your reply, duri.

Well, I want to make things clear, I used both 1st order and 2nd order differencing(simple upwind scheme and Fromm's scheme) and 1st order case was successfule while 2nd order case was unsuccessful.
Onim is offline   Reply With Quote

Old   September 5, 2013, 05:54
Default
  #5
Senior Member
 
duri
Join Date: May 2010
Posts: 130
Rep Power: 7
duri is on a distinguished road
It is not quite clear how you are applying limiter. It seems you are limiting the slope but comments says flux limiter. I hope you are using slope limiter. Stencil size for slope limiter will be 5, but your code says it is only 3. on left face slope need to be limited from i-2,i-1 and i-1,i on the right face it should be from i,i+1 and i+1,i+2.
duri is offline   Reply With Quote

Old   September 5, 2013, 06:04
Default
  #6
New Member
 
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 3
Onim is on a distinguished road
You`re right. Because I omitted definition of indices, it is still confusing. I`ll add their definition.
By the way, I don`t know what is slope limiter. I thought flux limiter gets slope rate as its argument and the result limits numerical flux. Are slope limiter and flux limiter different?
Onim is offline   Reply With Quote

Old   September 5, 2013, 06:11
Default
  #7
New Member
 
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 3
Onim is on a distinguished road
I added some definitions to the code.

Last edited by Onim; September 5, 2013 at 06:18. Reason: previous message was my misunderstanding.
Onim is offline   Reply With Quote

Old   September 10, 2013, 06:00
Default
  #8
New Member
 
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 3
Onim is on a distinguished road
I have made my program work, and I would like to post the reason for oscillation.

At the first time, when I "extrapolate" physical quantities from the cell center(?) to cell interface, I just calculated for each of them independently. That made these physical quantities be inconsistent to each other. Especially for the case of contact discontinuity, where some quantities are continuous while others are not, so on the interface some quantities are extrapolated at high order accuracy but others are at first order . That was resource of unphysical phenomenom. After I reconstruct all physical quantities on the interface from selected primitive variables which are extrapolated with flux limiter, the problem has gone away.

Additionally, I understood the difference between slope limiter and flux limiter, but the difference does not matter that much. We may find some discussion here. Slope limiter vs. Flux limiter
Onim is offline   Reply With Quote

Reply

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


Similar Threads
Thread Thread Starter Forum Replies Last Post
UDF compiling problem Wouter Fluent UDF and Scheme Programming 6 June 6, 2012 04:43
area does not match neighbour by ... % -- possible face ordering problem St.Pacholak OpenFOAM 9 November 22, 2011 11:02
Problem in implementing cht tilek CFX 3 May 8, 2011 08:39
natural convection problem for a CHT problem Se-Hee CFX 2 June 10, 2007 06:29
Adiabatic and Rotating wall (Convection problem) ParodDav CFX 5 April 29, 2007 19:13


All times are GMT -4. The time now is 14:28.