Welcome to our community

Be a part of something great, join today!

[SOLVED] setting up IC and BC

dwsmith

Well-known member
Feb 1, 2012
1,673
One can show that mass diffusion without chemical reactions obeys the same basic equation as heat conduction.
$$
D_{AB}\frac{\partial^2 C_A}{\partial x^2} = \frac{\partial C_A}{\partial t}
$$
where $C_A$ is the concentration of the species $A$ diffusing into a medium $B$ and $D_{AB}$ is the mass diffusivity.
Consider steel carburization.
In this high-temperature process carbon is diffused into the steel to achieve certain desirable material properties (e.g., increase tensile strength).
When the carburization is performed at a temperature $1200^{\circ}$C the value of the diffusion coefficient is approximately $D_{AB}\approx 5.6\times 10^{-10}$ $(\text{m}^2/\text{sec})$.
Suppose that a 1cm thick slab initially has a uniform carbon concentration of $C_A = 0.2\%$.
How much time is required to raise the carbon concentration at a depth of 1-mm below the slab surface to a value of 1%?

Is the IC $C_A(x,0) = .2$? What are the boundary conditions?
If I can identify those, I can solve it.
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
I don't think we can say what the BC's are at this point. Something is missing: the concentration of the source of carbon, as well as whether or not the carbon is going to be applied on both sides of the plate, or only one side. Does the problem not say this?
 

dwsmith

Well-known member
Feb 1, 2012
1,673
I don't think we can say what the BC's are at this point. Something is missing: the concentration of the source of carbon, as well as whether or not the carbon is going to be applied on both sides of the plate, or only one side. Does the problem not say this?
[Reminder: The transient solution requires homogeneous boundary conditions; the easiest way to achieve this is by rescaling/redefining the concentration variable.]

This is all that is left.
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
[Reminder: The transient solution requires homogeneous boundary conditions; the easiest way to achieve this is by rescaling/redefining the concentration variable.]

This is all that is left.
Then the concentration at the boundaries is zero for all time: you have perfect sink conditions.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
How do I write that?

$C_A(0,t) = C_A(1,t) = 0$.

So my IC right in the first post?
 
Last edited:

dwsmith

Well-known member
Feb 1, 2012
1,673
So if my IC and BC are correct, I have the solution as
$$
C_A(x,t) = \frac{.008}{\pi}\sum_{n = 1}^{\infty}\frac{\sin(2n-1)\pi x}{2n-1}\exp\left[-D_{AB}\pi^2 t (2n-1)^2\right]
$$
Correct?

How would I solve for the time it rakes to raise the carbon content to 1%?
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
How do I write that?

$C_A(0,t) = C_A(1,t) = 0$.

So my IC right in the first post?
The IC's in the OP are correct, and these BC's are correct.

So if my IC and BC are correct, I have the solution as
$$
C_A(x,t) = \frac{.008}{\pi}\sum_{n = 1}^{\infty}\frac{\sin(2n-1)\pi x}{2n-1}\exp\left[-D_{AB}\pi^2 t (2n-1)^2\right]
$$
Correct?

How would I solve for the time it rakes to raise the carbon content to 1%?
I have not checked this solution, but the form of it is plausible. I would first plug in your $x$ value. What's left might be able to be simplified by a CAS (maybe it's a geometric series, e.g.). If not, you could try taking the first 20 terms or so and solving numerically.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
The IC's in the OP are correct, and these BC's are correct.



I have not checked this solution, but the form of it is plausible. I would first plug in your $x$ value. What's left might be able to be simplified by a CAS (maybe it's a geometric series, e.g.). If not, you could try taking the first 20 terms or so and solving numerically.
What is my x value? If it is one, then the series is 0 for all t.
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
What is my x value? If it is one, then the series is 0 for all t.
Your slab is 1 cm thick, and you're interested in the concentration at 1 mm getting to 1%. So your $x=0.1$.

However, there is something puzzling me. Did the original question read "raise the carbon concentration to 1%"? Because, if it is correct that the boundary conditions are homogeneous perfect sinks, with a zero concentration, then I doubt the concentration will rise to 1%. The only way it could do that would be if there is enough residual concentration on the corresponding half of the slab so that the concentration near the middle pushes out to the depth $x=0.1$ much faster than the concentration near $x=0.1$ pushes out to the sink. That would be highly unlikely in an isotropic medium. Are you sure there isn't a source of carbon at the boundaries? Something with a concentration higher than 1%?
 

dwsmith

Well-known member
Feb 1, 2012
1,673
You should compare this to the analogous heat conduction problem. As
an example problem, consider of steel carburization. In this high-temperature process carbon
is diffused into the steel to achieve certain desirable material properties (e.g., increase
tensile strength). When the carburization is performed at a temperature 1200C the value
of the diffusion coefficient is approximately DAB  5.6 x 10^-10 (m^2=sec). Suppose that a
1cm thick slab initially has a uniform carbon concentration of CA = 0:2%. How much time
is required to raise the carbon concentration at a depth of 1-mm below the slab surface to a value
of 1%? [Reminder: The transient solution requires homogeneous boundary conditions; the
easiest way to achieve this is by rescaling/redefining the concentration variable.]

Here is the question.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
Each side of the slab is exposed to a carbon-rich environment with a constant concentration of CA = 1:5%.

I didn't see this sentence.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
The boundary conditions are $C_A(0,t) = C_A(1,t) = 0.015$ and the initial condition is $C_A(x,t) = .002$.
The solution to
$$
C_{A_{ss}} = 0.015.
$$
We are looking for solutions of the form $C_{A_{\text{trans}}}(x,t) = \varphi(x)\psi(t)$.
We have that
\begin{alignat}{3}
\varphi(x) & = & A\cos x\lambda + B\sin x\lambda\quad\quad\quad (1)
\end{alignat}
and
$$
\psi(t) = C\exp\left[-d_{AB}\lambda^2 t\right]
$$
where $d_{AB} = 5.6\times 10^{-6}$ $(\text{cm}^2/\text{sec})$.
The family of solutions for $\varphi(x)$ can be obtained by
\begin{alignat*}{3}
\varphi_1(0) = 1 & \quad & \varphi_2(0) = 0\\
\varphi_1'(0) = 0 & \quad & \varphi_2'(0) = 1
\end{alignat*}
This leads us to
\begin{alignat}{3}
\varphi(x) & = & A\cos x\lambda + B\frac{\sin x\lambda}{\lambda}\quad\quad\quad (2).
\end{alignat}
To see why this is correct, take $\lambda = 0$ for equation (1).
We would have $\varphi = A$ which is not the general solution of $\varphi'' = 0$.
The general solution is
$$
\varphi = A + Bx.
$$
If we take equation (2), we get
$$
\varphi = A + B\lim_{\lambda\to 0}\frac{x\sin x\lambda}{x\lambda} = A + Bx
$$
as needed.
Now, we can use the boundary conditions.
We have that
$$
\varphi_n(x) = \sin x\lambda = 0,\quad \lambda = \pi n, \quad n\in\mathbb{Z}.
$$
Our general solution is of the form
$$
C_{A_{\text{trans}}}(x,t) = \frac{1}{\pi}\sum_{n = 1}^{\infty}A_n\frac{\sin xn\pi}{n}\exp\left[-d_{AB}\lambda^2 t\right].
$$
Using $-0.013$ as our initial condition, we can solve
$$
A_n = \frac{-0.026}{n\pi}\int_0^1\sin xn\pi dx = \begin{cases}
0, & \text{if n is even}\\
0.052, & \text{if n is odd}
\end{cases}
$$
Therefore, the solution is
$$
C_A(x,t) = C_{A_{ss}} + C_{A_{\text{trans}}} = 0.015 + \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin x(2n - 1)\pi}{2n - 1}\exp\left[-d_{AB}(2n - 1)^2\pi^2 t\right].
$$
I have this now. Is it correct?
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
A few comments:

1. I would substitute $u=C_{A}-0.015$ as the new dependent variable. That will render your boundary conditions homogeneous, as the hint says.

2. I'm not sure your separation scheme was done correctly. Let's suppose we have
$$D\,\frac{\partial^{2}u}{\partial x^{2}}=\frac{\partial u}{\partial t},$$
and we let $u=X(x)\,T(t).$ Then the DE becomes
$$D\,X''\,T=X\,\dot{T}.$$
Dividing through by $XT$ yields
$$D\,\frac{X''}{X}=\frac{\dot{T}}{T}=-\lambda^{2}.$$
Both resulting equations do not need to have $D$ in them. You get
$$X''=-\frac{\lambda^{2}}{D}\,X,\quad \text{and}\quad \dot{T}=-\lambda^{2}T.$$

3. You need to be careful how your IC and BC's translate into the separated ODE's. Double-check how you did that in your last post.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
A few comments:

1. I would substitute $u=C_{A}-0.015$ as the new dependent variable. That will render your boundary conditions homogeneous, as the hint says.
I don't understand this.
 

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
I don't understand this.
Let's suppose that $u=C_{A}-0.015$. I chose $0.015$ because that is the boundary condition on both sides of the slab. If $C_{A}(0,t)=C_{A}(1,t)=0.015$, then $u(0,t)=u(1,t)=0$. See how that works?

Then,
$$\frac{\partial u}{\partial x}=\frac{\partial C_{A}}{\partial x},$$
and you can see that all the derivatives are equal because all I've done is shift by a constant.

Finally, you have to shift the IC as well, and say that the initial concentration in terms of $u$ is $u(x,0)=C_{A}(x,0)-0.015=-0.013.$

So this shift ends up not complicating either the DE or the IC's, and it gives you homogeneous BC's. So it's definitely a gain.

Does that make more sense? Apologies if I didn't explain things clearly before.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
Let's suppose that $u=C_{A}-0.015$. I chose $0.015$ because that is the boundary condition on both sides of the slab. If $C_{A}(0,t)=C_{A}(1,t)=0.015$, then $u(0,t)=u(1,t)=0$. See how that works?

Then,
$$\frac{\partial u}{\partial x}=\frac{\partial C_{A}}{\partial x},$$
and you can see that all the derivatives are equal because all I've done is shift by a constant.

Finally, you have to shift the IC as well, and say that the initial concentration in terms of $u$ is $u(x,0)=C_{A}(x,0)-0.015=-0.013.$

So this shift ends up not complicating either the DE or the IC's, and it gives you homogeneous BC's. So it's definitely a gain.

Does that make more sense? Apologies if I didn't explain things clearly before.
That is what I did though in my solution.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
I will elaborated on what I did.

So we have $u(0,t) = u(1,t) = 0.015$ and $u(x,0) = 0.002$.

$d_{AB} = 5.6\times 10^{-6}$ $(\text{cm}^2/\text{sec})$.

We have
$$
\varphi_n = A\cos\lambda x + B\frac{\sin\lambda x}{\lambda}
$$
and
$$
\psi = C\exp\left[-d_{AB}\lambda_n^2t\right]
$$
Solving for the steady state, we have $\varphi'' = 0$ so $\varphi = ax+b$.
The steady state solution is $C_{A_{ss}} = 0.015$.
Using the homogeneous bc now, we have
$$
\varphi_n = \sin\lambda = 0\Rightarrow \lambda=\pi n
$$
Our general solution is
$$
C_A(x,t) = 0.015 + \frac{1}{\pi}\sum_{n = 1}^{\infty}A_n\frac{\sin\pi nx}{n}\exp\left[-d_{AB}\pi^2n^2t\right]
$$
Using the IC,
$$
0.015 + \frac{1}{\pi}\sum_{n = 1}^{\infty}A_n\frac{\sin\pi nx}{n} = 0.002\iff -0.013 = \frac{1}{\pi}\sum_{n = 1}^{\infty}\frac{\sin\pi nx}{n}
$$
Solving for the Fourier coefficients
$$
-0.013\int_0^1\sin\pi n x dx = A_n\int_0^1\sin^2\pi n x dx\Rightarrow A_n = \begin{cases} 0, & \text{if n is even}\\
-0.052, & \text{if n is odd}\end{cases}
$$
We can disregard the $1/n\pi$ since it was already accounted for.
$$
C_A(x,t) = 0.015 - \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin\pi (2n-1)x}{2n-1}\exp\left[-d_{AB}\pi^2(2n-1)^2t\right]
$$
So that is my solution.

Now (1) is that correct? (2) I am still having trouble solving for what time raise the carbon in the steel to .01 at a depth of 1mm.
I can use software to solve for $t$ but I don't know how.

Suppose it is correct, I would have
$$
C_A(0.1,t) = 0.015 - \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin\frac{\pi}{10} (2n-1)}{2n-1}\exp\left[-d_{AB}\pi^2(2n-1)^2t\right] = 0.01
$$
$$\Rightarrow C_A(0.1,t) = \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin\frac{\pi}{10} (2n-1)}{2n-1}\exp\left[-d_{AB}\pi^2(2n-1)^2t\right] = 0.005
$$
 
Last edited:

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
The boundary conditions are $C_A(0,t) = C_A(1,t) = 0.015$ and the initial condition is $C_A(x,t) = .002$.
The solution to
$$
C_{A_{ss}} = 0.015.
$$
We are looking for solutions of the form $C_{A_{\text{trans}}}(x,t) = \varphi(x)\psi(t)$.
We have that
\begin{alignat}{3}
\varphi(x) & = & A\cos x\lambda + B\sin x\lambda\quad\quad\quad (1)
\end{alignat}
and
$$
\psi(t) = C\exp\left[-d_{AB}\lambda^2 t\right]
$$
where $d_{AB} = 5.6\times 10^{-6}$ $(\text{cm}^2/\text{sec})$.
The family of solutions for $\varphi(x)$ can be obtained by
\begin{alignat*}{3}
\varphi_1(0) = 1 & \quad & \varphi_2(0) = 0\\
\varphi_1'(0) = 0 & \quad & \varphi_2'(0) = 1
\end{alignat*}
This leads us to
\begin{alignat}{3}
\varphi(x) & = & A\cos x\lambda + B\frac{\sin x\lambda}{\lambda}\quad\quad\quad (2).
\end{alignat}
To see why this is correct, take $\lambda = 0$ for equation (1).
We would have $\varphi = A$ which is not the general solution of $\varphi'' = 0$.
The general solution is
$$
\varphi = A + Bx.
$$
If we take equation (2), we get
$$
\varphi = A + B\lim_{\lambda\to 0}\frac{x\sin x\lambda}{x\lambda} = A + Bx
$$
as needed.
Now, we can use the boundary conditions.
We have that
$$
\varphi_n(x) = \sin x\lambda = 0,\quad \lambda = \pi n, \quad n\in\mathbb{Z}.
$$
Our general solution is of the form
$$
C_{A_{\text{trans}}}(x,t) = \frac{1}{\pi}\sum_{n = 1}^{\infty}A_n\frac{\sin xn\pi}{n}\exp\left[-d_{AB}\lambda^2 t\right].
$$
Using $-0.013$ as our initial condition, we can solve
$$
A_n = \frac{-0.026}{n\pi}\int_0^1\sin xn\pi dx = \begin{cases}
0, & \text{if n is even}\\
0.052, & \text{if n is odd}
\end{cases}
$$
Therefore, the solution is
$$
C_A(x,t) = C_{A_{ss}} + C_{A_{\text{trans}}} = 0.015 + \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin x(2n - 1)\pi}{2n - 1}\exp\left[-d_{AB}(2n - 1)^2\pi^2 t\right].
$$
I have this now. Is it correct?
I see what you're doing now. Your transient solution is my $u$, and your steady-state solution is the shift.

I would agree with your solution completely except that I think this statement:

$$
A_n = \frac{-0.026}{n\pi}\int_0^1\sin (n\pi x) dx = \begin{cases}
0, & \text{if n is even}\\
0.052, & \text{if n is odd}
\end{cases}
$$
should be
$$
A_n = -0.026\int_0^1\sin (n\pi x) dx = \begin{cases}
0, & \text{if n is even}\\
-\frac{0.052}{n\pi}, & \text{if n is odd}
\end{cases},
$$
with the resulting solution
$$
C_A(x,t) = 0.015 - \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin [(2n - 1)\pi x]}{2n - 1}\exp\left[-d_{AB}(2n - 1)^2\pi^2 t\right].
$$

Now then: we are required to find the time necessary to raise the concentration at $x=0.1$ to a level of $0.01$. That is, find the time $t_{s}$ such that $C_{A}(0.1,t_{s})=0.01$. So, we require that
$$
0.01=C_A(0.1,t_{s}) = 0.015 - \frac{0.052}{\pi}\sum_{n = 1}^{\infty}\frac{\sin [(2n - 1)\pi (0.1)]}{2n - 1}\exp\left[-d_{AB}(2n - 1)^2\pi^2 t_{s}\right].
$$
I highly doubt that you can find this symbolically. The problem is that the target variable appears in a changing exponent in an infinite series. It's not a geometric sum. I think your best bet is to truncate the series to, say, 20 terms (although check to see how your solution changes if you add a few more terms), and solve numerically.
 

dwsmith

Well-known member
Feb 1, 2012
1,673
Mathematica tells me it cannot do it.

Solve[0.052/Pi*
Sum[Sin[Pi/10*(2 n - 1)]/(2 n - 1)*
E^{-5.6*10^{-6}*Pi^2*(2 n - 1)^2*t}, {n, 1, 20}] == 0.005, t]

Matlab wont solve it as well

solve(0.052/pi*symsum(sin(pi/10*(2*n-1))*exp(-5.6*10.^{-6}*(2*n-1).^2*pi.^2*t),n,1,20)==0.005,t)
It just says undefined variable n.
 
Last edited:

Ackbach

Indicium Physicus
Staff member
Jan 26, 2012
4,191
The Mathematica command you entered is the symbolic one. Try this:

\begin{align*}
&f[t\_]=0.0150-\frac{0.052}{\pi}\sum_{n=1}^{20}\left(\frac{\text{Sin}[(2n-1)\pi(0.1)]}{2n-1}\times \text{Exp}[-5.6\times 10^{-6}(2n-1)^{2}\pi^{2}t]\right)\\
&\text{Plot}[f[t],\{t,0,10000\}]\\
&\text{FindRoot}[f[t]-0.01,\{t,3000,4000\}].
\end{align*}