1D Diffusion Equations (Numerical Computing) [kind of long]

[Jelfish]

I am doing research with a professor and he asked me to graph some solutions to a 1D diffusion equations of the following form:

$$C_e \frac{\partial T_e}{\partial t}=\frac{\partial}{\partial z}\left( \kappa \frac{\partial T_e}{\partial z} \right) - g(T_e-T_i)+\alpha (1-R)If(t)\exp (-\alpha z)$$
$$C_i \frac{\partial T_i}{\partial t}=-g(T_e-T_i)$$

Let me describe to you what this equation is used for. My professor wishes to find the heat distribution of a piece of metal after a laser pulse hits it.

This equation deals with one dimension, meaning it only takes into account the depth of the metal slab, not the height or width.

Bad Visual Ahead:

|==    ^          o-o-o-o-o-o-o-o-o-o-o-o

The first object (on the left) is the laser. The "^" is a laser pulse. The line of o-o-o's is a latice. Note the 1 dimensioned nature.

The first differential equation deals with the temperature of the electrons and the second is the temperature of the lattice (the equations are coupled). This is called the "two-temperature model" because it accounts for the time delay between the transfer of energy between the electrons, which initially absorb the energy of the laser, and the lattice, which sends a vibration throughout the crystal.

Here's an explanation of the constants (i.e. mathematically, you can imagine them all to be, say, 0.5):

The electron heat diffusion coefficient is given by $$\kappa$$. The heat capacities of the electron and lattice are given by $$C_e$$ and $$C_i$$, respectively. The coupling strength is given by $$g$$. The damping coefficient and reflectivity of the metal in question is given by $$\alpha$$ and $$R$$, respectively. The maximum intensity of the pulse is given by $$I$$

Now, I'll explain the last term in the first equation, $$\alpha (1-R)If(t)\exp (-\alpha z)$$.

The function f(t) is the intensity of the pulse at a given instant. Here's a visual: Imagine a line of cars at a traffic light. Instead of the light turning green in an instant, the green light's intensity slowly increases, reaches a peak, and then dims again (like a traveling Gaussian distro. wave). This is how the laser pulse is modeled - therefore, f(t) is given by the laser used, but it will be some sort of Gaussian with respect to time.

The exponential function with respect to depth describes the intensity decay due to atoms being "blocked" by the atoms in front of them. Again, using the (poor) line of cars analogy, imagine that your visibility of the green light is partially blocked because of the cars in front of you and said blocking is dependant on how deep you are in the line of cars.



Now for the numerical computing part . I'm no expert at this at all so I've been trying to research this a bit. I've solved the diffusion equation through Fourier series, but this is different in that the boundary conditions are not what I'm use to.

Assumptions that could make no sense:
T(0) and T(infinity) = room temperature (assuming one pulse)
T'(0) and T'(infinity) = 0 (assuming one pulse)
These assumption are for both electron and lattice temperatures.


(If you're still reading THANK YOU!)

Ok - here's how I started:

To solve for T_e and T_i numerically, the method of finite difference was used (I know, it's not the best, but I'm learning). For a general diffusion equation, T(z,t),

$$\frac{\partial T}{\partial t}=\frac{\partial^2 T}{\partial z^2}$$,

the numerical solution, given in terms of finite differentials, delta t and delta z, and at a given point (i,j) is:

$$\frac{T(i,j+\delta t) - T(i,j)}{\delta t} = \frac{T(i-\delta z,j)-2T(i,j)+T(i+\delta z,j)}{(\delta z)^2}$$.

Therefore, the function T can be found by interating over i and j given some initial boundary conditions:

$$T(i,j+\delta t) = T(i,j) + \delta t \frac{T(i-\delta z,j)-2T(i,j)+T(i+\delta z,j)}{(\delta z)^2}$$.

By substituting into the original equations:

$$T_e(i,j+\delta t)=T_e(i,j) + \frac{\delta t}{C_e} \biggl[\kappa \frac{T(i-\delta z,j)-2T(i,j)+T(i+\delta z,j)}{(\delta z)^2} - g\bigl(T_e(i,j)-T_i(i,j)\bigl)+\alpha (1-R)If(j)\exp (-\alpha i)\biggl]$$.

And for the second one:

$$T_i(i,j+\delta t)=T_i(i,j)-\frac{\delta t}{C_i}g\bigl(T_e(i,j)-T_i(i,j)\bigl)$$


Now, before people start asking for consultation fees, I wanted to know if this much made sense. If so, goodie. This is sort of my num-comp 'theory.' If this is wrong, there's no sense in displaying my also-would-be-wrong MatLab code.


Thank you VERY MUCH in advance!

 

[PerennialII]

Hi Jelfish,

interesting application, squeezing it all out from the heat equation !

The "2-T" approach looks neat, using the lattice as thermal inertia and doing the coupling via "convection", let alone the nice bcs terms.

I'd say for the initial & boundary conditions setting the temperature to room and assuming thermal insulation, other than the actual bcs driving the solution, will work fine.

Your difference derivation looks fine by me, if you get "bored" with the method could do the same with a 1D linear finite element .

 

[Jelfish]

I'll look into that and other methods as well. I'm becoming more and more certain that my troubles lie in my coding, as my graphs seem to show random spikes into infinity. Thanks for the reply.

/me switches back to Matlab

 

[PerennialII]

I'll look into that and other methods as well. I'm becoming more and more certain that my troubles lie in my coding, as my graphs seem to show random spikes into infinity. Thanks for the reply.

/me switches back to Matlab

... yeah, one of the problems with FDMs IMO is their implementation ... not really a fan .

 

[Clausius2]

Nice problem. Sure you will have lot of fun integrating it .

You should be careful about using an explicit or implicit method. I have noticed your method is explicit, which usually will be numerically instable. Try to solve the problem implicitly in time.

 

[chitambira]

Seems Jellyfish you are working on some similar kind of thing to what i want to do. Please could you help me with coding for the finite differences? If you can send me your code for one of your equations that can help me. look at my posts and you can see where i am i am having problems. I haven't coded any FD before so I don't know how to implemenet the iterations in both time and space.

 

[Integral]

chitambira,
Note that the last response before yours in this thread is over 3yrs old. Do not expect anyone to available to continue the discussion with you.

 

[bilgealp]

Hi chitambira,
This is a very popular problem and a somewhat old paper using FD in 2-dimensions may help.
Noye, B.J., Tan, H.H., "Finite Difference Methods for solving Two-Dimensional Advection-Diffusion Equation", International Journal For Numerical Methods in Fluids, vol.9, 75-98, 1989.

This paper was quite useful during a homework a couple of years ago. I have programmed this diffusion problem several times and it will quite be an interesting study for you.
The main difficulty lies in the discretization of the first term on the right hand side, namely

$$\frac{\partial}{\partial z}\left( \kappa \frac{\partial T_e}{\partial z} \right)$$

but it is quite straightforward. You can find the solution in several texts. I was not a FD pro at the time of developing a program for 2D version of this problem but the following book helped me a lot.
Kowalik, Z. and T. S. Murty, 1993. Numerical modeling of ocean dynamics, World Scientific Publ., 481 pp.

 

[bilgealp]

Okkey. Some more info on that. Assuming that you do not have constant diffusion coefficient, i.e. $$\kappa$$ is inside the derivative with respect to z
$$\frac{\partial}{\partial z}\left( \kappa \frac{\partial U}{\partial z} \right)$$
then the best solution would be a staggered grid system in space. Divide the domain into $$n$$ discrete parts and on the halfway to each discrete ends (o's), put (x's) as shown in the configuration below.
o-x-o-x-o
Let the distance between consecutive o's is $$\Delta x$$ and label o's from one end to the other with the sequence 1,2,3,...,i-1,i,i+1,...,n-1,n. Since x's are at the midpoints, label them as 1+1/2,2+1/2,...,i-1/2,i+1/2,...,n-1/2. Then the discrete form for this term is approximately as
$$\frac{\partial}{\partial z}\left( \kappa \frac{\partial U}{\partial z} \right) \approx \frac{\kappa_{i+1/2}\frac{U_{i+1}-U_i}{\Delta x}-\kappa_{i-1/2}\frac{U_i-U_{i-1}}{\Delta x}}{\Delta x}$$
If that would help...