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

In summary, the conversation discusses the 1D diffusion equations used to graph the heat distribution of a piece of metal after a laser pulse hits it. The equations account for the temperature of the electrons and the lattice, and the time delay between energy transfer. The constants in the equations represent the electron heat diffusion coefficient, heat capacities of the electron and lattice, coupling strength, damping coefficient, reflectivity of the metal, and maximum intensity of the laser pulse. The final term in the first equation models the intensity of the pulse over time and depth. The conversation also discusses the use of numerical methods, specifically the method of finite difference, to solve the equations. The boundary conditions for the equations are set to room temperature and thermal insulation. The implementation
  • #1
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 [tex]\kappa[/tex]. The heat capacities of the electron and lattice are given by [tex]C_e[/tex] and [tex]C_i[/tex], respectively. The coupling strength is given by [tex]g[/tex]. The damping coefficient and reflectivity of the metal in question is given by [tex]\alpha[/tex] and [tex]R[/tex], respectively. The maximum intensity of the pulse is given by [tex]I[/tex]

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

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 :cry: . 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),

[tex]\frac{\partial T}{\partial t}=\frac{\partial^2 T}{\partial z^2}[/tex],

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

[tex]\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}[/tex].

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

[tex]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}[/tex].

By substituting into the original equations:

[tex]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][/tex].

And for the second one:

[tex]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)[/tex]

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!
Physics news on Phys.org
  • #2
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 :biggrin: .
  • #3
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
  • #4
Jelfish said:
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 :wink: .
  • #5
Nice problem. Sure you will have lot of fun integrating it :biggrin: .

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.
  • #6
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.
  • #7
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.
  • #8
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

[tex]\frac{\partial}{\partial z}\left( \kappa \frac{\partial T_e}{\partial z} \right)[/tex]

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.
Last edited:
  • #9
Okkey. Some more info on that. Assuming that you do not have constant diffusion coefficient, i.e. [tex]\kappa[/tex] 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 [tex]n[/tex] discrete parts and on the halfway to each discrete ends (o's), put (x's) as shown in the configuration below.
Let the distance between consecutive o's is [tex]\Delta x[/tex] 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
[tex]\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}[/tex]
If that would help...

1. What is a 1D diffusion equation?

A 1D diffusion equation is a mathematical model that describes the process of diffusion in one dimension. It is used to understand the spread of particles or substances in a medium, such as the movement of heat or chemicals in a solid or gas.

2. How is a 1D diffusion equation solved numerically?

A 1D diffusion equation can be solved numerically using various techniques, such as finite difference methods or finite element methods. These methods involve discretizing the equation into smaller parts and using iterative calculations to approximate the solution.

3. What are the key assumptions made in a 1D diffusion equation?

The key assumptions in a 1D diffusion equation include a constant diffusion coefficient, no external forces acting on the system, and the absence of any significant chemical reactions or reactions with other substances.

4. What are some applications of 1D diffusion equations?

1D diffusion equations have various real-world applications, such as modeling the spread of pollutants in the environment, predicting the behavior of heat and mass transfer in industrial processes, and understanding the diffusion of drugs in biological systems.

5. How do 1D diffusion equations differ from other types of diffusion equations?

1D diffusion equations are specific to a one-dimensional system, meaning that the diffusion only occurs in one direction. Other types of diffusion equations, such as 2D and 3D, take into account diffusion in multiple directions. Additionally, 1D diffusion equations are simpler and easier to solve compared to higher-dimensional equations.

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