CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > System Analysis

Sods problem: Oscillations with WENO schemes

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 31, 2015, 17:57
Default Sods problem: Oscillations with WENO schemes
  #1
New Member
 
Join Date: May 2015
Posts: 1
Rep Power: 0
Ironiker is on a distinguished road
Hey everyone!

When applying my WENO code for the Sods Shock Tube problem, I get some minor oscillations in the solutions:


I am using a central compact WENO scheme with SSP time stepping, which works well in the scalar case (Transport, Burgers, etc.). Initial data and numerical constants are taken from the paper
http://epubs.siam.org/doi/abs/10.1137/S1064827599359461

Thus the mistake must be in the numerical process of cell updating/ calculation the right hand side. I am just copying the necessary code fragments to maintain clarity.

- U is a matrix evaluated rho, m and E at every cell midpoint.
- The function WENO gives back the reconstructed value at the left side and at the right side for each cell for u1, u2, and u3.
- To calculate solutions at cell boundaries I am shifting the solution to the right, calculate the numerical flux and then shift that to the left due to conservation.
- F and the Jacobian matrix are written down explicitly and checked twice.

Code:
U = [rho; m; E]; 

%%%%% Calculating righthandside for time stepping %%%%%

function du = rhs(u)
vL2,vR2 = WENO(U)
z = f_num(shmatrixperiodl(vR2),vL2);
y = shmatrixperiodr(z);
du = -(y-z);
end

%%%%% Shift function routines %%%%%

function v = shmatrixperiodl(U)
v = [U(:,1) U(:,1:(end-1))];
end


function v = shmatrixperiodr(U)
v = [U(:,2:end) U(:,end)];
end

%%%%% Calculating numerical flux %%%%%

function Z = f_num(A,B)
% Numerical flux function
% Here: local lambdax-Friedrichs flux
[m,n]=size(A);
C=sparse(m,n);
for i=1:n
C(:,i) = max(abs(eig(Jacobmatrix(A(:,i)))), abs(eig(Jacobmatrix(B(:,i)))) );
end
Z = (F(A)+F(B)-C.*(B-A))./2;

%%%%% Flux Function %%%%%

function Y = F(U)
global gamma
Y = [U(2,:); 
     0.5*(3-gamma)*U(2,:).^2./U(1,:)+(gamma-1)*U(3,:); 
     gamma*U(3,:).*U(2,:)./U(1,:)+0.5*(1-gamma)*U(2,:).^3./(U(1,:).^2)]; 
end


%%%%% Jacobian %%%%%%

function Y = Jacobmatrix(U)
global gamma
Y = [ 0                                      1                                    0 
    -0.5*(3-gamma)*U(2).^2./(U(1).^2)      (3-gamma)*U(2)./U(1)                  (gamma-1);
    -gamma*U(3).*U(2)./(U(1).^2)+(gamma-1)*U(2).^3./(U(1).^3)  gamma*U(3)./U(1)+1.5*(1-gamma)*U(2).^2/(U(1).^2)     gamma*U(2)./U(1)];
end
I would be glad for some helpful advices!
Ironiker is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Adiabatic and Rotating wall (Convection problem) ParodDav CFX 5 April 29, 2007 20:13
Colocated schemes : decoupling problem Sylvain Main CFD Forum 0 October 18, 2003 13:10
extremely simple problem... can you solve it properly? Mikhail Main CFD Forum 40 September 9, 1999 10:11
Is this problem well posed? Thomas P. Abraham Main CFD Forum 5 September 8, 1999 15:52
Higher-order bounded convection schemes for pure advection with discontinuity Anthony Main CFD Forum 3 June 13, 1999 03:36


All times are GMT -4. The time now is 08:07.