Let the symbol “a” represent the wave speed. “a” is defined by the derivative

of the flux “f” with respect to the conserved variable “u”, So, for the Burgers

equation, since f = u x u /2, a = u.

At a shock, the characteristics converge and at an expansion, characteristics

diverge. We can use this fact to separate shock and expansion in 1-D and enforce

the entropy fix only for the expansion, without spoiling the shock. Otherwise,

since entropy fix is nothing but dissipation addition, the exact shock capturing

property of Roe scheme will be spoilt.

Let j and j+1 represent the left and right centroids of the cell interface at j+1/2.

Case 1 : If a(j) > 0 and a(j+1) < 0, then it is a shock.

Case 2 : If a(j) < 0 and a(j+1) > 0, then it is an expansion fan.

Enforce entropy fix only at the expansion.

The entropy fix is nothing but the quadratic variation of dissipation with respect to

the wave speed, instead of a linear variation as in Roe scheme. In Roe scheme, as

the wave speed goes to zero, the dissipation also goes to zero and this leads to

central differencing. That is why it gives a problem at zeor wave speeds (sonic

points). Harten’s correction keeps non-zero dissipation at zero wave speeds.

I will post about the roe’s method in the next few days.

The matlab module is as follows.

lamda is delt/delx

**THIS IS FOR 1D BURGER’S EQUATION**

if u(i-1,j-1)~=u(i-1,j)

a1=0.5*(u(i-1,j)+u(i-1,j-1));

else

a1=u(i-1,j-1);

end

if u(i-1,j)~=u(i-1,j+1) %to define the wave velocity flux

a2=0.5*(u(i-1,j+1)+u(i-1,j));

else

a2=u(i-1,j);

end

a1=0.5*(u(i-1,j)+u(i-1,j-1));

else

a1=u(i-1,j-1);

end

if u(i-1,j)~=u(i-1,j+1) %to define the wave velocity flux

a2=0.5*(u(i-1,j+1)+u(i-1,j));

else

a2=u(i-1,j);

end

s2=0.5*(u(i-1,j+1)-u(i-1,j));%defining the free parameter

s1=0.5*(u(i-1,j)-u(i-1,j+1));

if(a1<s1 && a1>-s1) %defining the local value for u.

e1=((a1*a1+s1*s1)/(2*s1));

else

if (a1>0)

e1=a1;

else

e1=-a1;

end

end

if (a2<s2 && a2>-s2)

e2=((a2*a2+s2*s2)/(2*s2));

else

if (a2>0)

e2=a2;

else

e2=-a2;

end

end

u(i,j)=u(i-1,j)-lamda*0.25*(u(i-1,j+1)*u(i-1,j+1)-u(i-1,j-1)*u(i-1,j-1))+

lamda*0.5*(e2*(u(i-1,j+1)-u(i-1,j))-e1*(u(i-1,j)-u(i-1,j-1)));

Advertisements

Please comment your code adequately