|
[Sponsors] |
September 5, 2013, 05:27 |
AUSM+, MUSCL coding problem
|
#1 |
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 12 |
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, 06:33 |
|
#2 |
Senior Member
duri
Join Date: May 2010
Posts: 245
Rep Power: 16 |
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, 06:34 |
|
#3 |
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 12 |
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 07:09. Reason: comment was not up to date |
|
September 5, 2013, 06:38 |
|
#4 |
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 12 |
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, 06:54 |
|
#5 |
Senior Member
duri
Join Date: May 2010
Posts: 245
Rep Power: 16 |
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, 07:04 |
|
#6 |
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 12 |
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, 07:11 |
|
#7 |
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 12 |
I added some definitions to the code.
Last edited by Onim; September 5, 2013 at 07:18. Reason: previous message was my misunderstanding. |
|
September 10, 2013, 07:00 |
|
#8 |
New Member
Minho
Join Date: Sep 2013
Posts: 6
Rep Power: 12 |
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 | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
area does not match neighbour by ... % -- possible face ordering problem | St.Pacholak | OpenFOAM | 10 | February 7, 2024 22:50 |
UDF compiling problem | Wouter | Fluent UDF and Scheme Programming | 6 | June 6, 2012 05:43 |
Problem in implementing cht | tilek | CFX | 3 | May 8, 2011 09:39 |
natural convection problem for a CHT problem | Se-Hee | CFX | 2 | June 10, 2007 07:29 |
Adiabatic and Rotating wall (Convection problem) | ParodDav | CFX | 5 | April 29, 2007 20:13 |