NDSolve

• To: mathgroup at smc.vnet.net
• Subject: [mg91372] NDSolve
• From: Uli Wuerfel <uli.wuerfel at fmf.uni-freiburg.de>
• Date: Tue, 19 Aug 2008 07:13:32 -0400 (EDT)

```Hallo,
I was trying to solve the heat equation with NDSolve(one spatial dimension, x, and time t). The simulated device (centered around x=0) consists of three regions: Two rather big volumina (with thickness each being dg) to the left and right of a thin layer (with thickness dp). A heat source is only present in the thin layer. The outer left boundary of the whole device is open, i.e. the spatial derivative of the temperature is zero there. At the right boundary, there is perfect cooling, i.e. the temperature stays constant.
As the heat is generated only in the thin layer, the change is quite abrupt. What I would like to know is if there is a way to tell Mathematica to place a very fine grid there because the result I get is not correct.

Here is the code:
dp=1e-5;
dg=1e-3;
eta=500/dp;

sol=u/.First[
NDSolve[D[u[t, x],t]==ldrp[x,0.5*dp,0.1*dp]*D[u[t,x],(x,2)]+eta*heat[x,0.5*dp,0.1*dp],
u[0,x]==0,D[u[t,-dg-dp/2],x]==0,u[t,dg+dp/2]==0},u,{t,0,tmax},{x,-dg-dp/2,dg+dp/2},MaxSteps->25000,MaxStepFraction->0.0002,PrecisionGoal->2]]

with
heat[x,a,b]:=-1+1/(Exp[(x - a)/b]+1)+1/(Exp[-(x+a)/b]
being the (steep) function describing the extension of the heat source and
ldrp[x,a,b]:=kg-(kg-kp)*heat[x,a,b]
describing the spatial characteristics of the heat conductivity divided by the product of mass density and heat capacity.
The result depends unfortunately very much on the thickness of the thin layer, although the TOTAL heat generated there is always the same (as eta=500/dp).
Also, the steady-state temperature distribution in the right volume (with the cooling at the outer right boundary) where there is no heat generation should be
just a function of the total generated heat and the heat conductivity of that volume. So the temperature difference
should (in the steady-state) always be:
u[x=dp/2]-u[x=dg+dp/2]=500/lambda_g,
with lambda_g being the heat conductivity of that volume.
But the results produced by Mathematica are different.
Do you know what I am doing wrong or is it just that it
is not possible to handle such abrupt changes?