# AUSM+, MUSCL coding problem

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

September 5, 2013, 04:27
AUSM+, MUSCL coding problem
#1
New Member

Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 6

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
 sod.jpg (26.4 KB, 13 views)

 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 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

 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 i-2,i-1 and i-1,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.cfd-online.com/Forums/mai...x-limiter.html

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Wouter Fluent UDF and Scheme Programming 6 June 6, 2012 04:43 St.Pacholak OpenFOAM 9 November 22, 2011 11:02 tilek CFX 3 May 8, 2011 08:39 Se-Hee CFX 2 June 10, 2007 06:29 ParodDav CFX 5 April 29, 2007 19:13

All times are GMT -4. The time now is 01:55.