Laplace Transform Method:

Diffusion Equation With Neumann Condition

M.M. Yovanovich

LPTDENEU.MWS

Laplace transform of one-dimensional diffusion equation in half-space with Neumann boundary condition at . The PDE is with initial condition , and boundary conditions: and .

> restart:

> with(inttrans):

> pde:= diff(u(x,t),x,x) = 1/alpha*diff(u(x,t),t);

> alias(U(x,s)=laplace(u(x,t),t,s)):

> ode:= laplace(pde,t,s);

> ic:= u(x,0) = 0;

> bc1:= subs(x=0, diff(u(x,t),x)) = -q[0]/k;

> bc2:= u(infinity,t) = 0;

> ode1:= subs(ic,ode);

> ode2:= subs(U(x,s) = U(x), ode1);

> solode2:= dsolve(ode2, U(x));

To have a bounded solution at , eliminate the positive exponential.

> _C1:= 0;

> solode3:= solode2;

Take the Laplace transform of the boundary condition at .

> bc1:= laplace(diff(u(x,t),x),t,s) = - laplace(q[0]/k, t,s);

Obtain the second constant of integration from the Neumann boundary condition.

> rhs(solode3);
_C2:= solve(simplify(subs(x=0, diff(%, x))) =
- q[0]/k/s, _C2);

> solode4:= solode3;

Error, invalid terms in product

The solution of the problem is the inverse Laplace tranform of the above relation.

The Laplace transform tables cannot provide the inverse directly from the above expression. It will be necessary to let before the tables can be used. We will try to use Maple to get the inverse.

> try1:= (q[0]/k)*sqrt(alpha)*exp(-sqrt(s)*x/sqrt(alpha))/
s^(3/2);

> invlaplace(try1,s,t);

Maple cannot provide the inverse of the expression. Use the Table presented to you. In the Table we find that the inverse Laplace transform of is . Comparing the table expression and the expression of the problem, we find that and that the transform is the product in the problem.

The solution can b e written in the form:

> u(x,t):= q[0]/k*(2*sqrt(alpha*t/Pi)*exp(-x^2/4/alpha/t) - x*erfc(x/2/sqrt(alpha*t)));

The solution can be non-dimensionalized by multiplying by and dividing by the product . Also introduce the similarity parameter: .

> ustar:= k*u(x,t)/(2*q[0]*sqrt(alpha*t));

> ustar:= 1/sqrt(Pi)*exp(-eta^2) + eta*erfc(eta);

The solution is now complete.