Numerical solution of non-linear diffusion equation

[maka89]

Hello. I have some questions regarding the equation:
$k\frac{\partial}{\partial x}\left( u\frac{\partial (u-r(x,t))}{\partial x} \right) = \kappa \frac{\partial u}{\partial t}$.
u is positive. r(x,t) is given as an input.

I have implemented this non-linear diffusion equation using backward Euler time stepping and a spatial scheme as if it were a standard diffusion equation : $\frac{\partial}{\partial x}\left( u\frac{\partial u}{\partial x} \right) \approx \frac{1}{\Delta x}\left[u_{i+1/2}\frac{u_{i+1}-u_i}{\Delta x}- u_{i-1/2}\frac{u_{i}-u_{i-1}}{\Delta x}\right]$.

Now, my numerical solution works. But not as good as I was hoping. I want it to work over longer time-periods than I am allowed to with my numerical solution within any reasonable time. When trying to use the solution for these kinds of time-periods, the integral over u is not conserved(In my case: The total amount of oil changes during simulation). I hypotize that it is due to the oil front(where u goes from u>0 to u=0). At each timestep I put all the nodes with negative u values to 0 to get a stable and non-negative solution. Some times, in order to make my solution work better, I have to increase the number of spatial nodes, to get a better resolution for the front.

I am thus looking for better schemes to implement my equation. Accuracy is not important, but stability and conservation of u is(At least more or less. The integral of u is allowed to vary with say 10-20% if it is a really long simulation)

Questions:
1. Is there any way to change the problem into one that is more convenient to solve numerically by, say, a change of variables, or applying a transform to the equation? I have been looking at the Kirchhoff transform, to transform the diffusion equation to a linear one. Would that be applicable to this problem?

2. Which numerical scheme would you use to solve the equation, if solving it as is?

Thanks :)

 

[Chestermiller]

Hi. I see you're still working on this.

It might help to know what r looks like. Maybe you can provide a graph. I would like to see how r varies with x and t over the domain of interest. For example, it would be of interest to know if r is always positive, and to know how r compares in magnitude with u. Maybe seeing some results would be helpful.

Do you also treat r implicitly in the difference formulation.

Chet

 

[maka89]

Hey! Yes I am still working on it ^^
r(x,t) is an arbritrary function that the user inputs to my simulation function. It can be anything, given that $\frac{\partial r}{\partial x}$ is not too steep at any point(This will give problems). The user gives two arrays representing $r(x,t=0)$ and $r(x,t=t_{max})$, representing the start and the end of the simulation. In between, r(x,t) is assumed to change linearly with time at each node. I treat r(x,t) implicitly as well. Only the gradient of r is used so hard to say anything about its magnitude. Which is larger of the gradients depends on where you look on the solution, but a lot of my problem is that the gradient of r(x,t) is too large. If I use a flat or nearly flat r(x,t), the solution works for longer time periods

Try e.g. : $r(x,t) = 100sin(pi \frac{x}{x_{max}})-30exp(-(\frac{(x-x_{max})/2}{x_{max}/8})^2) - 50(\frac{x}{x_{max}})(\frac{t}{t_{max}})$. Unit meters. The domain is x from 0 to xmax and t from 0 to tmax. In this case use e.g. xmax somewhere between 10^3 and 10^4 m.

 

[bigfooted]

The backward Euler method is still first order accurate in time. If you combine it with the forward Euler method, you get the Crank-Nicolson method, which is second order accurate in time. It was developed specifically for diffusion equations:
https://en.wikipedia.org/wiki/Crank–Nicolson_method

 

[maka89]

Yes, I will try the CN scheme and see what happens. That's an easy change to implement too! But I was thinking more if there were other spatial discretizations or fancy methods using e.g. flux limiters that are suitable for this kind of problem? Hard to find any good literature that treats similar equations and arent too mathematical for me to bother ^^

 

[Chestermiller]

Hey! Yes I am still working on it ^^
r(x,t) is an arbritrary function that the user inputs to my simulation function. It can be anything, given that $\frac{\partial r}{\partial x}$ is not too steep at any point(This will give problems). The user gives two arrays representing $r(x,t=0)$ and $r(x,t=t_{max})$, representing the start and the end of the simulation. In between, r(x,t) is assumed to change linearly with time at each node. I treat r(x,t) implicitly as well. Only the gradient of r is used so hard to say anything about its magnitude. Which is larger of the gradients depends on where you look on the solution, but a lot of my problem is that the gradient of r(x,t) is too large. If I use a flat or nearly flat r(x,t), the solution works for longer time periods

Try e.g. : $r(x,t) = 100sin(pi \frac{x}{x_{max}})-30exp(-(\frac{(x-x_{max})/2}{x_{max}/8})^2) - 50(\frac{x}{x_{max}})(\frac{t}{t_{max}})$. Unit meters. The domain is x from 0 to xmax and t from 0 to tmax. In this case use e.g. xmax somewhere between 10^3 and 10^4 m.

Please check over the equation for r again to make sure it is correct (particularly that middle term and the time factor). I don't want to start thinking about this some more until you're certain.

What are the units on k, kappa, and u? What is the initial condition on u? What are the boundary conditions? I want to reduce the equations to dimensionless form.

 

[maka89]

I actually managed to fix my code, so that the problem now is gone(I hope, need to test some more next week).
I used the following strategy to ensure that u is always positive(if i dont, it gets unstable):
-At each timestep, after calculating the new solution: for all i: u->max(0,u)
This removes the negative values of u at certain nodes. However, if the scheme was conservative to begin with, this will make sure that the mass increases, which was my problem.

I fixed it with the following solution, which basically normalizes the solution after removing the "negative mass":
1. calculate total amount of "negative" mass removed. I.e. mass added.
2. subtract a constant $\Delta u$from all nodes that have positive u values, so that it all adds up to the total added mass.
3. Check if any nodes has negative u values after subtraction. If yes: Repeat 1-3. If no: Done.

Maybe this is a bit shady, but I am not looking for high accuracy, and this makes the solution stable and gives reasonable results as far as I can see. Plus, if I increase the number of timesteps, the solution will go back to the same one as if I didnt do this normalization procedure.

If you still want to try solving this equation, Chestermiller, here is the info you requested:
k is dimensionless, u is meters, kappa is seconds/meters. I use no-flow BC. Initial condition: Try $u(x,0)=\begin{cases} 0, & \text{if }x< x_{max}/2 \\ max(r(x,0) - 72, 0), & \text{if }x\geq x_{max}/2 \end{cases}$.

There was one error in the middle term of r(x,t): $r(x,t) = 100sin(pi \frac{x}{x_{max}})-30exp(-(\frac{(x-x_{max}/2)}{x_{max}/8})^2) - 50(\frac{x}{x_{max}})(\frac{t}{t_{max}})$

 

[Chestermiller]

What are the values of $\kappa$ and k?

I've worked the differential equation into the following form for your consideration:
$$\kappa\frac{\partial y}{\partial t}=\frac{k}{2}\frac{\partial ^2y^2}{\partial x^2}+k\frac{\partial}{\partial x}\left(r\frac{\partial y}{\partial x}\right)+\kappa f(x)$$where y = (u - r) and $f(x)=-\frac{\partial r}{\partial t}$

I'd like to reduce the equation to dimensionless form, so I need kappa and k.

 

[maka89]

$\kappa = 1.53 10^{-7}$ and $k = 10^{-13}$ is a good place to start

 

[Chestermiller]

When you say that there is a no-flow boundary condition, do you mean that du/dx = 0, even though udr/dx is not? Or do you mean du/dx=dr/dx?

 

[maka89]

The last mentioned, I would think. The way I discretized it, mentioned in the first post, I only have to drop the term corresponding to the flux across the boundary. I think you might get problem with mass conservation if using a different spatial discretization scheme than that, but maybe there are other ways around this.

 

[Chestermiller]

How do you solve the problem using backward Euler in time if the equation is non-linear in u? Do you iterate at each time step until a corrector equation is satisfied to within a certain error tolerance?

Do you have access to an automatic ODE solver with the capability of handling stiff sets of ODEs? There are such integrators available for free over the internet such as DASSL. The IMSL library of subroutines also has a Gear integrator. You should really be solving this problem with a stiff integrator using the Method of Lines. Then you wouldn't have all the numerical problems you are encountering. And, it's really easy to program your problem when employing such software.

Chet