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

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.

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.

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.

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.

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

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.

A FTCS finite difference representation of Burgers' equation is:

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

where

The equivalent FTCS representation is:

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

where and

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.

The following tridiagonal algorithm can be constructed:

The tridiagonal form can be written as

where

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

.

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

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.

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.

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

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

In logarithm representation, we have the following linear evolution:

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.

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.

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.