Schemes and Shocks
- The Cole-Hopf transformation
As we now have developed an algorithm designed to solve Burgersí equation, it is useful to have some test cases to validate the code's behavior. Yet Burgersí non linearity is an obstacle for analytical resolution and exact solution are often unknown. To overcome this difficulty it is possible to use a little trick: the Hopf-Cole transformation. This transformation gives exact solution for Burgersí equation. Therefore it is a useful way to validate the routine and to check the correct evolution of the numerical results.
The Hopf cole transformation introduces a function F that is related to u by the following relation:
As a result Burgers' equation can be rewritten using F :
We can recognize the diffusion equation. For appropriate boundary conditions an exact solution of Burgers is therefore known.
It can be shown that the Cole Hopf transformation can also be applied to the 2D equation. In that case the number and the diversity of solution makes it particularly attractive.
In the 1D case if we consider , we have:
The diffusion equation can be written as:
with k=n and j=-1/(2n)
Using these limitations it is obvious that the solution of the diffusion equation is a solution for Burgers. It is also clear that a flux (Neumann) boundary condition for F is transformed in a Dirichlet condition for u.
In 1D steady solution of the diffusion equation are:
For this study we used two settings:
We obtained the following evolution using these boundary conditions and using as initial condition u(x)=-2n. The result are quite satisfactory as the solution converges towards the expected value.
In that case too we have convergence towards the exact solution.
Remark: these computation were done using the generalized Crank Nicholson
implicit scheme with d=1/6 and q=0,5. Further in this report it is shown
that it is a very good choice, but for these computation all implicit settings
are giving good results.
The propagating shock is a good start. The initial shock moves as a result of the advection term. The diffusion term tends to smoothen the solution finally leading to to decaying of the shock.
The discontinuity makes this initialization particularly suitable for the comparing of the schemes setting. The numerical dispersion of each scheme can be pointed out as it leads to the apparition of wiggles. It seems quite clear that with an increasing Reynolds the shock remains sharper because of the small diffusion and is therefore more difficult to have a correct evolution of the solution because of the small number of point available to describe the strong gradients. It must be observed the the Re is the critical number as it compares the advection term to the diffusion term.
Four main settings are used:
It appears that two scheme are giving outstanding results the 4 pts upwind and the generalized one with a small advantage to the later one.
For this setting of dt (0.001s) and dx (0.0033m) we can roughly consider as minimum acceptable values of viscosity:
Using the sinusoidal initialization we now expect the formation of a
shock at x=0,5.
We do observe a shock building up for all except the highest values of the viscosity. In this case too it is obvious that when the Re gets critically high the right tuning of the scheme is necessary to get an accurate solution.
After the shock has formed u decays towards 0 and the maximum are moving
away from 0.5 (the curve for t=4s).
Evolution with the viscosity
As we can see the viscosity is the key parameter for the solution behavior in particular considering the shocks formation and shape.
The shock width in respect with time is expected to be:
The shock width being the distance between maximums at a given time.
And a trendline:
Remark: these trend lines are only verified on a small portion
of the viscosity domain (viscosity between 0.01 and 0.001).