
[Sponsors] 
September 5, 2013, 04:27 
AUSM+, MUSCL coding problem

#1 
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6 
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? 

September 5, 2013, 05:33 

#2 
Senior Member
duri
Join Date: May 2010
Posts: 160
Rep Power: 9 
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. 

September 5, 2013, 05:34 

#3 
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6 
I attach some MUSCL code for your information.
Code:
#define indcn i,j,k #define indup iish,jjsh,kksh #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 iish,jjsh,kksh 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 = downwindcenter denom = centerupwind c make sure division by 0 does not happen if(abs(nom).lt.1d14)then nom = 0.0d0 denom = 1.0d0 elseif(nom.gt.1d14.and.abs(denom).lt.1d14)then nom = 1d14 denom = 1d0 elseif(nom.lt.1d14.and.abs(denom).lt.1d14)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 

September 5, 2013, 05:38 

#4 
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6 
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. 

September 5, 2013, 05:54 

#5 
Senior Member
duri
Join Date: May 2010
Posts: 160
Rep Power: 9 
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 i2,i1 and i1,i on the right face it should be from i,i+1 and i+1,i+2.


September 5, 2013, 06:04 

#6 
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6 
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? 

September 5, 2013, 06:11 

#7 
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6 
I added some definitions to the code.
Last edited by Onim; September 5, 2013 at 06:18. Reason: previous message was my misunderstanding. 

September 10, 2013, 06:00 

#8 
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6 
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. http://www.cfdonline.com/Forums/mai...xlimiter.html 

Thread Tools  
Display Modes  


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  SeHee  CFX  2  June 10, 2007 06:29 
Adiabatic and Rotating wall (Convection problem)  ParodDav  CFX  5  April 29, 2007 19:13 