Unsteady State Diffusion in Slab

...okay. So this is mainly because I forgot EVERYTHING as far as differential equations go. Writing helps me sort things out. So while I go over my notes, I am going to write the steps down.

The first problem is this:

Consider unsteady state diffusion of a species into a slab. The thickness of the slab is b. At time before zero, the slab is at a uniform concentration Co. At time t greater than zero the surfaces of the slab are maintained at concentration C1. The PDE of concern is;

(∂C/∂t)=D(∂2C/∂x2)

where the LHS term accounts for the time dependency and the RHS term for diffusion in the slab. Since the problem is symmetric we can solve for half the slab, using the following BCs, after placing the coordinate system at the middle of the slab:

#1) t=<0, C=Co, all x
#2) t>0, C=C1, x=b/2
#3) any t, (∂C/∂x)=0, x=0

Using separation of variables obtain the solution for C.

The solution?


#1) MAKE DIMENSIONLESS

The first thing you're going to do is to make everything dimensionless. This will make it less of a hassel to solve and clear out some of the terms.

Let u = (C-C1)/(Co-C1)

This is handy in several different ways. When you have a concentration of Co, your value is 1 and when you have a concentration of C1, your concentration is 0. At first glance, this might not be easy to see, but this helps your equation become homogenous, and we love solving homogenous equations. It's a lot easier.

Let X = x/(b/2)

Since our coordinates are in the middle, this makes the middle position equal to 1. Again, it's easier.

Let T = 4D/b^2*t

This can be a little harder to see. Consider that the equation, before the T substitution, looks like this:

(∂u/∂t)=D*(b/2)^2*(∂2u/∂X2)

We just made u and X dimensionless. Now, manipulate t so that it equals the other side. You should get 4D/b^2 on the opposite. Because everything else is dimensionless besides t and these constants, then we can assume that these constants together must be dimensionless.

After all of this, our equation looks like this:

(∂u/∂T)=(∂2u/∂X2)

It's getting simpler! :D


#2) WRITE NEW BOUNDARY CONDITIONS

Now we need new boundary conditions. We write:

#1) T=<0, u=1, all X
#2) T>0, u=0, X=1
#3) any T, (∂u/∂X)=0, X=0

Easy! Just plug and play for this part. And yes, it's really simple, but don't forget this! It's essential for completely the problem satisfactorarily.


#3) SPLIT "U" UP INTO TWO PARTS

This is one place where I seem to get stuck on or forget, though really it's not that difficult.

Let u = F(X)G(T)

This means that

(∂u/∂T)=(∂2u/∂X2)

is now

(∂(FG)/∂T)=(∂2(FG)/∂X2)

where FG = F(X)G(T)

As you probably already guessed by now, on the LHS, F(X) is a constant and G(Y) on the RHS is a constant. So pull it out of the differential.

We now have:

F(∂G/∂T)=G(∂2F/∂X2)

Divide both sides by FG.

We now have:

(1/G)(dG/dT)=(1/F)(d2F/dX2)

Also! Don't forget our favorite little eigenvalue! :D

(1/G)(dG/dT)=(1/F)(d2F/dX2) = -λ^2

Lambda is minus because otherwise this problem would go to infinity and that's really bad. We want it to go down, after all. :P

If you have no idea what I'm talking about, you'll see next section...


#4) SOLVE YOUR ODES!

If you missed the last step or completely blind or whatever, look up! Now, instead of having a weird partial differential equations, we have awesome ordinary differential equations. Awesome? I thought so.

Anyway, then it just becomes a simple process to solve. Let's look at G(T) first! :D

(1/G)(dG/dT) = -λ^2

You've probably memorized the solution, at this point. If you haven't, it's this:

G(T) = k*e^(-λ^2*T)

where k is an arbitrary constant for now (we solve for it later). If you don't know how to get there, it's a relatively simple first-order ODE you can solve for fun, if you want. :P

For F(X), it looks like this:

(1/F)(d2F/dX2) = -λ^2

The solution of this second-order ODE is this:

F(X) = Acos(λX) + Bsin(λX)

Again, this should be fairly simple to derive, if you care. It's not hard in the least.

Remember that we let u = F(X)G(T). Therefore:

u = k*e^(-λ^2*T) * (Acos(λX) + Bsin(λX))

Or:

u = e^(-λ^2*T) * (Pcos(λX) + Qsin(λX))

where P and Q are just arbitrary constants.


#5 APPLY BOUNDARY CONDITIONS

Now we get to make this problem, which has so far been extremely mathematical and not much else, actually worthwhile to the real world! Huzzah! Remember that we are looking at our dimensionless boundary conditions. This will make our lives easier. With that, let's take a look at BC#3!

BC#3

any T, (∂u/∂X)=0, X=0

Using basic calculus:

u = e^(-λ^2*T) * (Pcos(λX) + Qsin(λX))

(∂u/∂X) = e^(-λ^2*T) * (-Pλsin(λX) + Qλcos(λX))

For X=0:

e^(-λ^2*T) * (-Pλsin(0) + Qλcos(0)) = 0

sin(0) = 0, cos(0) = 1

Simplifying:

e^(-λ^2*T) * (Qλ) = 0

Therefore, Q = 0.

So now our problem is even simpler! :D Our new equation for u is this:

u = e^(-λ^2*T) * (Pcos(λX))

Now let's apply another boundary condition! How about BC#2?

BC#2:

T>0, u=0, X=1

This time, we don't even have to do calculus! :D

u = e^(-λ^2*T) * (Pcos(λX))

u = e^(-λ^2*T) * (Pcos(λ)) = 0

This puts us in a little bit of a pickle. As you probably know, the cosine function is periodic, which means there are an infinite amount of numbers that make it 0. However, the numbers always fall under the pattern:

λ = π/2, 3π/2, 5π/2...

An easier way to express that is the following:

λ = λn = (1/2)(2n+1)π, where n = 0, 1, 2, 3...

Therefore:

u = e^(-λ^2*T) * (Pncos(λnX))

Now, let's do boundary condition #1! :D

BC#1

T=<0, u=1

This is the trickiest one, but don't freak out! Here's where we use the definition of orthogonality. First, assume that T=0. Then:

u = sum(Pncos(λnX),X,0,∞) = 1

The function cos(λnX) is an orthogonal function.

Dun dun dunnnnnn!

So how are we ever going to solve Pn, if we have an infinite amount of values? Simple! Well... not simple, but it's possible to solve for Pn anyway. The way we are going to do that is to use what we know about orthogonality and solve for it that way.


ORTHOGONALITY AND THE KRONECKER DELTA

This really ties in with the boundary conditions, but for some reason, orthogonality is hard for me to understand, so I'm going to write it out long, which means it gets its own section.

Orthogonality basically says this:

∫(fi(x)*fj(x),x,0,1)

= 1, if i=j

= 0, if i=/=j

Remember, our function is this:

u = sum(Pncos(λnX),X,0,∞) = 1

If this is indeed orthogonal, then:

∫((1)cos(λnX)dX,0,1) = Pn∫((1)(cos(λnX))^2dX,0,1)

Where (1) is the weighting function for cartesian coordinates. Solving for Pn:

Pn = ∫((1)cos(λnX)dX,0,1) / ∫((1)(cos(λnX))^2dX,0,1)

This answer is fine for the test, which is what this is essentially for. However, using several trig identities, this can be further simplified to:

Pn = (-1)^n*(2/λn)

Thus, the solution now looks like this:

u = e^(-λ^2*T) * ((-1)^n*(2/λn)cos(λnX))

This is a much better solution! However, we're still not done.


CONVERT FROM "U" TO "C"

Remember that your original problem had to do with concentration! So now you have to substitute your original values in. Recall:

u = (C-C1)/(Co-C1)

X = x/(b/2)

T = 4D/b^2*t

Now, substitute these values in your solution:

u = e^(-λ^2*T) * ((-1)^n*(2/λn)cos(λnX))

(C-C1)/(Co-C1) = e^((-λ^2*4D/b^2)*t) * ((-1)^n*(2/λn)cos(2λnx/b))

Now you're done.


SO WHAT DID WE JUST DO?

Yeah. So some of you might be saying, "Wow, that math was easy, but what the hell did we just do?" Good question!

The problem asked us to find the PDE solution to a diffusion scenario. This meant that, for this problem, there was a solid box that had something diffuse into it. What the box was and what sort of material that diffused into the box really doesn't matter. That's covered in our constant "D" anyway. What matters to us is how fast this material moved through the box and how far the material moved through the box. This is what the PDE is trying to model. To measure diffusion, we looked at concentration. As the material moves through the box, the concentration of the material changes, depending on time and length.

Looking at our final solution:

(C-C1)/(Co-C1) = e^((-λ^2*4D/b^2)*t) * ((-1)^n*(2/λn)cos(2λnx/b))

...we see that both time and length play an important role into figuring out what concentration is. Intuitively, this is correct. As C goes to C1, we see that time runs to infinity and the change of concentration for the slab essentially stops. This is also intuitive. As the concentration goes to C1, we expect that there will be no change in concentration, unless another factor was there. For x, the variables surrounding that are periodic, which means it changes constantly as it goes through the slab. This also makes sense. In the long run, what determines the change of concentration through the slab is time, not the length of the slab.

This, of course, is a very simplistic problem. But these basics form many other worthier problems in the real world.