INTRODUCTION


The purpose of this work is to study the Burgers' equation:

wpe2.jpg (3053 octets)

This nonlinear equation, very similar to the Navier-Stokes equation, is a useful model for numerical experiments. An interesting test case with shock formation is provided by the time evolution of a sinusoidal wave profile.To solve the Burgers' equation, the finite-difference method is used and is programmed in FORTRAN.

We can also explore an unstable behaviour by adding terms to this equation.

 

PHYSICAL BEHAVIOUR OF BURGERS' EQUATION

The Burgers' equation could serve as a nonlinear analog of the Navier-Stokes equations. This single equation have a convective term, a diffusive term and a time-dependent term.

 

eq.gif (3027 octets)

Burgers' equation is parabolic when the viscous term is included. If the viscous term is neglected, the remaining equation is hyperbolic.

If the viscous term is dropped from the Burgers' equation the nonlinearity allows discontinuous solutions to develop. The way that this can occur is illustrated schematically in the following figure. A wave is convecting from left to right and solutions for succesive times are indicated. Points on the wave with larger values of u convect faster and consequently overtake parts of the wave convecting with smaller values of u. It is necessary to postulate a shock accross which u change discontinuously to have a unique solution and so a physically result.

ibe.gif (4963 octets)

The comparable wave development for the 'viscous' Burgers' equation is shown in the following figure. The effect of the viscous term is twofold. First , it reduces the amplitude of the wave for increasing time. Second, it prevents multivalued solutions from developping.

be.gif (4144 octets)

The above features make Burgers' equation a very suitable model for testing computational algorithms for flows where severe gradients or shocks are anticipated.

 

DISCRETIZATION METHODS

I choose to solve the Burgers' equation using a finite-difference technique.  I used first a FTCS scheme obtained by applying forward-time and centered-space differences. This explicit scheme is very easy to program but fails to give a correct solution when the viscosity is too low. Indeed stability conditions need to be respected. To avoid unstable solutions, I programed an implicit Cranck-Nicholson scheme. Of course this scheme can't be used for very low viscosties.

Explicit scheme

A FTCS finite difference representation of Burgers' equation is:

wpe4.jpg (6328 octets)

There is an alternative way of handling the nonlinear convective term. This requires rewriting the equation in conservation form,

wpe7.jpg (2900 octets)   where   wpe6.jpg (1291 octets)

The equivalent FTCS representation is:

wpe5.jpg (5721 octets)

 

Implicit scheme

Applying implicit schemes to nonlinear equations is not quite as straightforward as for linear equations. A Crank-Nicolson implicit formulation is:

wpe8.jpg (5059 octets)

where  wpe17.jpg (1859 octets) and  wpe9.jpg (3893 octets)    

To make use of the very efficient Thomas algorithm, it is necessary to have a system of linear tridiagonal equations for the solutions . The appearance of the nonlinear implicit term  poses a problem. However, this problem can be overcome by using a Taylor series expansion.

wpeB.jpg (5044 octets)

      wpeC.jpg (2756 octets)                          wpeD.jpg (2405 octets)

The following tridiagonal algorithm can be constructed:

wpeF.jpg (5277 octets)

The tridiagonal form can be written as

wpe10.jpg (3040 octets)

where

wpe11.jpg (2326 octets) 

wpe12.jpg (1314 octets) 

wpe13.jpg (2416 octets)

 wpe14.jpg (3358 octets)

wpe15.jpg (1571 octets) 

This implicit scheme is unconditionally stable in the Von Neumann sense and has a truncation error of

wpe16.jpg (1727 octets).

 

PROGRAM

Here is the listing of the program in FORTRAN 90. The user can use both methods described above.

ccccccccccccccccccccccccccccccccccccccccccccccccccc

                         Burgers' equation routine

ccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
REAL U1,U2,r,dt,dx,nu,tf,t,c,X,a,b,d,c2,d2,k,cl,lambda,alpha
integer n,i,j,scheme
ALLOCATABLE U1(:),U2(:),X(:),a(:),b(:),c(:),d(:),c2(:),d2(:)
CHARACTER*80 DON
CHARACTER*80 param
WRITE(*,*) 'Fichier de donnees ? '
READ(*,'(A)') param
OPEN(4,FILE=param)
READ(4,*) n
READ(4,*) dt
READ(4,*) tf
READ(4,*) nu
READ(4,*) scheme
READ(4,*) cl
READ(4,*) lambda
READ(4,*) alpha
CLOSE(4)
IF(scheme.EQ.1) THEN
WRITE(*,*) 'Explicit scheme'
ELSE
WRITE(*,*) 'Implicit scheme'
ENDIF
IF(cl.EQ.0) THEN
WRITE(*,*) 'cl rigide '
ELSE
WRITE(*,*) 'cl libre'
ENDIF
WRITE(*,*) 'n =',n
WRITE(*,*) 'dt =',dt
WRITE(*,*) 'tf =',tf
WRITE(*,*) 'nu =',nu
WRITE(*,*) 'lambda =',lambda
WRITE(*,*) 'alpha =',alpha

dx=1./float(n)
r=nu*dt/(dx*dx)
k=dt/(2*dx)

ALLOCATE (U1(n+1),U2(n+1),X(n+1),a(n+1),b(n+1),c(n+1),
& d(n+1),c2(n+1),d2(n+1))
DO i=1,n+1
X(i)=(i-1)*dx
a(i)=0
b(i)=0
c(i)=0
d(i)=0
END DO
c initialisation
DO i=1,n+1
U1(i)=SIN(2*3.1415927*(i-1)/n)
END DO

t=0
DO WHILE (t.LT.tf)
t=t+dt
c conditions limites
U1(1)=0
U1(n+1)=0
U2(1)=0
U2(n+1)=0


c boucle explicite
IF(scheme.EQ.1) THEN
DO i=2,n
U2(i)=r*U1(i+1)+(1-2*r)*U1(i)+r*U1(i-1)
& -k*U1(i)*(U1(i+1)-U1(i-1))
END DO
DO i=1,n
U1(i)=U2(i)
END DO

c boucle implicite
ELSE

DO j=2,n+1
a(j)=-0.25*dt*U1(j-1)/dx-0.5*r
END DO
DO j=1,n+1
b(j)=1+r
END DO
DO j=1,n
c(j)=0.25*dt*U1(j+1)/dx-0.5*r
END DO
d(1)=(1-r)*U1(1)+0.5*r*U1(2)
d(n+1)=0.5*r*U1(n)+(1-r)*U1(n+1)
DO j=2,n
d(j)=0.5*r*U1(j-1)+(1-r)*U1(j)+0.5*r*U1(j+1)
& +dt*(lambda*U1(j)+alpha*U1(j)*U1(j)*U1(j))
END DO

c algorithme de Thomas
c2(1)=c(1)/b(1)
d2(1)=d(1)/b(1)
DO i=2,n+1
c2(i)=c(i)/(b(i)-a(i)*c2(i-1))
d2(i)=(d(i)-a(i)*d2(i-1))/(b(i)-a(i)*c2(i-1))
END DO
U1(n+1)=d2(n+1)
DO i=1,n
U1(n+1-i)=d2(n+1-i)-U1(n-i+2)*c2(n+1-i)
END DO
ENDIF
END DO


WRITE(*,*) 'Fichier de resultat ? '
READ(*,'(A)') DON
OPEN(5,FILE=DON)
DO i=1,n+1
WRITE(5,*) X(i),U1(i)
END DO
CLOSE (5)

END

NUMERICAL SOLUTION

As said previously,an interesting test case for unsteady flows with shock formation is provided by the time evolution of a sinusoidal wave profile. u is set to the initial condition sin(2px), x included in [0,1]. Rigid boundary conditions, i-e u(0)=u(1)=0, are used. Viscosity is set to 1e-3 m^2/s, the grid to 500 points and the time step to 1e-4 s.  Only the Implicit sheme is used to avoid instability problems.

The graph gives the evolution of the solution. Each point propagates with a different speed and leads to the formation of a shock. This shock is then damped owing to the dissipative term.

evo_t.gif (11377 octets)

If the viscosity is too low ( n<1e-4 ), wiggles appears and eventually causes the solution to blow up. For lower viscosities, others numerical schemes must be used. In the following figure n=1e-4 m^2/s: at the top and the bottom of the shock wiggles appears.

VISCOSITY VARIATION

By increasing the viscosity, the expected shock fades. The following graph shows the  solution obtained for different viscosities at time t=0.3 s.

p2.gif (14428 octets)

 

We can show that the distance e between the two crests follows the law:

epsilon.gif (5691 octets)

In logarithm representation, we have the following linear evolution:

p1.gif (5179 octets)

 

INSTABILITY STUDY

We can carry out a instability study by adding others terms to the Burgers' equation.

Let's take unstable values for l=4 and a=1 (n=500,dt=1e-4,n=1e-3 m^2/s). As expected we observe that the shock is no more damped and the solution blows up very quickly.

stab.gif (9108 octets)

 The following graph (t=0,1s, n=500,dt=1e-4, n=1e-2, a=1) gives the numerical solutions for two different values of l: l>0 (unstable) and l<0 (stable). In black the solution of the Burger' equation is given to highlight this difference.

CONCLUSION

The Implicit Cranck-Nicholson scheme proves to be efficient to solve burgers' equations. However for too low viscosities others schemes must be used. This work shows that instabilities can have numerical and physical reasons. Paying attention to both problems is then required.

The Burgers' equation allows to understand the process of shocks formation. The effect of the viscous term is to reduce the amplitude of the wave for increasing time and then can prevent shocks from forming.