Nonlinear heat equation -- Handling the conductivity

[maka89]

Hey! I'm currently solving the heat equation using finite differences. I have a conductivity $k(u)$ that varies greatly with temperature. It even drops to zero at $u=0$.
I have discretized the equations the following way:

$\frac{\partial}{\partial x}\left( k(u) \frac{\partial u}{\partial x}\right) = \alpha \frac{\partial u}{\partial t}$

$\frac{1}{\Delta x^2}\left( k_{i+\frac{1}{2}}(u_{i+1}-u_i) - k_{i-\frac{1}{2}}(u_{i}-u_{i-1}) \right) = \alpha \frac{\partial u}{\partial t}$

My question is basically how to choose $k_{i+1/2}$ in a clever way. One way is to use $k_{i+1/2} = (k_{i}+k_{i+1})/2$(If k is discretized). Another to use $k_{i+1/2} = k(x_i + \Delta_x/2)$.

But the method I am considering is to choose $k_{i+1/2}$ in such a way that the total heat moved between the gridpoints is equal to $k_{i+1/2}(u_{i+1}-u_{i})$, assuming that u(x) is linear between the gridpoints.

Thus we get: $k_{i+1/2}\frac{(u_{i+1}-u_{i})}{\Delta x}\Delta x = \int_{u_i}^{u_{i+1}} k(u) \frac{\partial u}{\partial x} dx = \int_{u_i}^{u_{i+1}} k(u) du$ or $k_{i+1/2} = \frac{1}{(u_{i+1}-u_{i})} \int_{u_i}^{u_{i+1}} k(u) du$.

What do you guys think of the approach? Is it reasonable? This essentially becomes a finite volume method does it not? Will it work even in two dimensions, without further complications?

 

[Chestermiller]

Your method looks kind of interesting, although I'm not sure it will be much more accurate than just using k at the adjacent average temperatures. Don't forget that this would give the exact same result as the value at the average temperature if k varies linearly with u. So, you've got to wonder if it's really worth the effort. Also, it isn't immediately obvious to me how this would be extended to 2D.

Chet

 

[maka89]

Thanks for the answer and the bump.

In reservoir simulation they use a not so different approach to the one i gave above to conclude that the harmonic average should be used instead of the arithmethic for the permeabilities in Darcy's law( I.e. equivalent to the conductivity in the heat equation.) if it varies over space as discrete steps, as these macro-parameters often do in practical use. I agree that it is not very useful in my example.

I haven't pursued this any more, since it was an attempt to to improve on a numerical solution that wasn't working properly, becoming unstable and not conserving heat. I now know that this was due to an advection term not included in the equation above(I wanted to state the problem as generic as possible, little did i know that I was leaving out the problem itself xD). The equation should be:$\frac{\partial}{\partial x}\left( k(u)\frac{\partial(u-c(x,t))}{\partial x}\right) = \alpha\frac{\partial u}{\partial t}$. Choosing $k_{i+1/2} = 0.5(k_i+k_{i+1})$ should give a conservative scheme, so this was not my problem. Rather it was that central differences(which i need to have a conservative scheme) doesn't work well for advection-diffusion equations if the advection flux becomes of the same order of magnitude or larger than the diffusion flux.

I think the right thing for me to do is to include artificial diffusion.

Does anyone have any tips on how to do this should be done, given that this is a 2d problem? Can I just augment all the $k_{i+1/2}$ for the diffusion term, so that I get an acceptable Peclet number (advective flux divided by diffusive flux)? Is it OK to not just add a constant artificial diffusion, but rather let it vary over space?

 

[Chestermiller]

I can give you lots of advice on this, because I have worked on similar problems, both in heat transfer and in reservoir modeling (i.e., groundwater).

Let's start with that transport equation. It doesn't look correct to me. The advection term doesn't look familiar to me. What are you doing there?

Regarding use of artificial diffusion, I would definitely not use that.

Are you solving the time integration fully implicit?

Chet

 

[maka89]

OK. So the equation describes flow of oil, not heat(just wanted to make the question generic, since the math is basically the same). It is a 3D problem that has been reduced to 2D. $k(u)$ is a relative permeability. I.e. When there is no oil, it is zero, when there is only oil and no water, it is equal to the permeability of the rock. The advection term arises from buoyancy of the oil against the top barrier of the layer where the oil flows, $z(x,y,t)$ which is not entirely flat, and changes with time.

I use an ADI scheme. But as of now I keep $k(u)$ explicit. It works well in most cases. However, when one of my parameters becomes small, the advection term becomes stronger compared to the other term, and I get trouble.

Do you think making $k(u)$ implicit as well would make the problem more stable?
I am looking for stability and mass conservation. Accuracy is not that crucial as long as the general features of the solution are intact.

 

[Chestermiller]

OK. So the equation describes flow of oil, not heat(just wanted to make the question generic, since the math is basically the same). It is a 3D problem that has been reduced to 2D. $k(u)$ is a relative permeability. I.e. When there is no oil, it is zero, when there is only oil and no water, it is equal to the permeability of the rock. The advection term arises from buoyancy of the oil against the top barrier of the layer where the oil flows, $z(x,y,t)$ which is not entirely flat, and changes with time.

I use an ADI scheme. But as of now I keep $k(u)$ explicit. It works well in most cases. However, when one of my parameters becomes small, the advection term becomes stronger compared to the other term, and I get trouble.

Do you think making $k(u)$ implicit as well would make the problem more stable?

I don't think so.

I am looking for stability and mass conservation. Accuracy is not that crucial as long as the general features of the solution are intact.

Is it possible to go fully implicit, or is the problem too large for that?

My experience has been single phase, associated with deepwell injection of hazardous waste. Some of the modeling we did took into account density driven flow in cases where the injectate density (typically dilute aqueous solution) is of lower or higher density than that of the formation fluid.

I have always found advection diffusion equations very nasty, and it is an area where, in my opinion, no satisfactory solution has ever been found. In the work that we did on underground injection, we struggled to deal with it all the time, mostly by going to smaller time steps.

How important is the diffusion? If it can be approximately neglected, it is possible to consider solving the transport by the method of characteristics.

Chet

 

[Strum]

It sounds to me like it is still single phase. Also if you have a nicely resolved lattice perhaps you could solve the outer derivative using a Fourier transforms and get the conservation equation on u that way.

What is c(x,t) ?

 

[maka89]

I don't think so.

Is it possible to go fully implicit, or is the problem too large for that?

How important is the diffusion? If it can be approximately neglected, it is possible to consider solving the transport by the method of characteristics.

I tested going fully implicit, but it didn't help much. The diffusion may be unimportant in some areas but in others, the diffusion and advection/buoyancy need to cancel each other out. If not, all the oil would accumulate in a delta function of a local maximum of the layer, instead of having a horizontal barrier between oil and water. Tried a bit with artificial diffusion, to smooth out the moving front, but it looks like it won't help much either. Guess I need to crank up the number of timesteps :/

It sounds to me like it is still single phase. Also if you have a nicely resolved lattice perhaps you could solve the outer derivative using a Fourier transforms and get the conservation equation on u that way.

What is c(x,t) ?

The equation only looks at the oil phase, yes. But it is derived from two-phase equations, assuming that the pressure fluctuations in the water phase are small, due to a higher permeability. c(x,t) is the gradient of the layer, which creates flow together with the buoyancy of the oil.
This method sounds interesting. I am using a regular 2d grid. Do you have any links to examples or references?

 

[Chestermiller]

Going fully implicit is usually supposed to help, but not if vΔt ≥ Δx. I forgot what this limitiation is called.

I take it back. Fully implicit is supposed to be unconditionally stable irrespective of the step size if the problem is linear (from all the discussions I've read up on). Even if the equations are not linear, I'm surprised that the fully implicit scheme was not stable. Certainly, it should allow you to go to much larger time steps.

Chet

 

[maka89]

Yes, I was already taking rather large timesteps, though ^^ Solution did not conserve mass well enough for practical use, but it did not crash, so I guess that makes it stable?

 

[Chestermiller]

Yes, I was already taking rather large timesteps, though ^^ Solution did not conserve mass well enough for practical use, but it did not crash, so I guess that makes it stable?

Is it possible to set up the spatial finite difference scheme so that it automatically conserves mass?

 

[maka89]

Not to my knowledge... I guess using larger step sizes/control volumes where there is expected to be more rapid flow would help some?(at the price of accuracy)

 

[Chestermiller]

Not to my knowledge... I guess using larger step sizes/control volumes where there is expected to be more rapid flow would help some?(at the price of accuracy)

What I'm saying is that, when I set up finite difference schemes in problems like this, I choose a scheme that I can show automatically conserves mass at the finite difference level (while still being second order). Would you like to try me out on a couple of terms to see how I would do it?

 

[maka89]

Sure!

 

[Strum]

Don't we still need an equation on c(x,t)? Or perhaps it is just specified somehow?

 

[maka89]

The c(x,t) field is one of the input parameters for the model :)

 

[Chestermiller]

The c(x,t) field is one of the input parameters for the model :)

Is that term in the equation presented in post #4 really what you are calling the advection term? It doesn't look like the kind of advection terms I'm used to seeing. How come it is being multiplied by k if it is advection? Shouldn't it involve u?

Chet

 

[maka89]

You're right, the u dependence comes from k(u)... I don't know if that is enough to call it an advection term technically or not... But definately not your standard advection case.

I have made a bit more progress on the problem, by making some simplifications in my model, making k(u) not change so rapidly. I have a bit more flexibility with what parameters to choose now. However, my numerical solution has the dubious property that it sometime needs a lot of gridpoints(not timesteps, but gridpoints) to converge. This strengthens my suspicion that it is the discontinuous edges where u drops to zero that is the problem...

 

[Chestermiller]

OK. So, for whatever reason, your advection term is $\frac{\partial}{\partial x}\left(k(u)\frac{\partial c}{\partial x}\right)=\frac{\partial (f(x,t)k(u))}{\partial x}$. Here's a finite difference representation that automatically conserves mass at the finite difference level:

$$\frac{f(x+(\Delta x /2)\bar{k}_{x+(\Delta x /2)}-f(x-(\Delta x /2)\bar{k}_{x-(\Delta x /2)}}{\Delta x}$$
where
$$\bar{k}_{x+(\Delta x /2)}=\frac{\int_{u(x)}^{u(x+\Delta x)}{k(u)du}}{u(x+\Delta x)-u(x)}$$

Chet

 

[maka89]

Thanks! Will try it out ^^