Message Boards Message Boards

0
|
3407 Views
|
5 Replies
|
3 Total Likes
View groups...
Share
Share this post:

Problem with NDSolve

Posted 10 years ago

Hi. I'm solving a PDE that is partially elliptic. So I need to solve the equation iteratively. My problem is at the end of the code. I'm doing interactions manually in order to evaluate the code. At the second interaction, I don't know why NDSolve doesn't accept the new boundary condition. Please, find comments at the body of the code. Could some one help me?

ClearAll["Global`*"]
(*All here are input like properties*)
vm = 5;
r0 = 127.25;
di = r0*2;
de = 273;
v = 0.000000291;
[Rho]f = 957.8;
[Rho]t = 7854;
pr = 1.76;
kf = 0.68;
kt = 57.7;
cf = 4217;
ct = 472.7;
re = (vm*(2*r0/1000)/v)
e = 0.00015;
f = 1.325 (Log[0.27*(e/(2*r0/1000)) + 5.74*(1/re)^0.9])^(-2)
tinf = 298;
tme = 373;
[Alpha]f = kf/([Rho]f*cf)
c1 = (-f*[Rho]f*vm^2/(4*r0/1000))
g = 9.81;
[Beta] = 1/((tme + tinf)/2)
var = 0.00001921;
prar = 0.7023;
[Alpha] = 0.0000337;
kar = 0.0289;

(*algebric equation for the turbulent model*) 
ym[r_] := (((r0/1000) - r)*vm*(f/8)^(0.5))/v;
[CurlyEpsilon]m[r_] := v*(0.4*ym[r]/6)*(1 + (r/(r0/1000)))*(1 + 2*(r/(r0/1000))^2);

(*this solve the momentum equation giving the velocity profile. there is a discontinuity at r=0,that's why the equation is  integrated from r=0.0000001*)
momento = NDSolve[{(1/r)*D[(v + [CurlyEpsilon]m[r])*r*D[u[r], r], r] == (1/[Rho]f)*c1, u'[0.0000001] == 0, u[(r0/1000)] == 0}, u[r], {r, 0.0000001, (r0/1000)}, MaxSteps -> 10000]; 

mom[r_] = u[r] /. First[momento];

(*this calculates the average velocity*)
umed = (2/((r0/1000)^2))*NIntegrate[(mom[r]*r), {r, 0.0000001, r0/1000}]

(*this calculates the internal convection*)
hi = (kf/(2*r0/1000))*((re - 1000)*pr*(f/8))/(1 + 12.7*(f/8)^0.5*(pr^(2/3) - 1))

(*this calculates the external convection*)

ra = (g*[Beta]*(tme - tinf)*(de/1000)^3)/(var*[Alpha]);
he = (kar/(de/1000))*(0.6 + (0.387*ra^(1/6))/(1 + (0.559/prar)^(9/16))^(8/27))^2

(*here are the integration limits and initial values*)
l = 160;
tempo = 80;
Clear[tmed];
Clear[tmedpred];
Clear[ttf];
tmed[x_, t_] = 310;
tmedpred[x_, t_] = tinf;
ttf[x_, r_, t_] = tinf;
subrel = 0.2;
erro = 10;
eps = 10^-4;

(*While[erro[GreaterEqual]eps,*)

enf1 = NDSolve[{([Rho]f*cf*(D[tf[x, r, t], t] + mom[r]*D[tf[x, r, t], x])) == ((1/r)*D[(r*([Alpha]f + [CurlyEpsilon]m[r])*D[tf[x, r, t], r]), r]), tf[x, r, 0] == (tme - tinf)*Exp[-1000*x] + tinf, tf[0, r, t] == tme, (D[tf[x, r, t], r] /. r -> 0.0000001) == 0, -kf*(D[tf[x, r, t], r] /. r -> (r0/1000)) == (1 - Exp[-1000*x])*(1 - Exp[-1000*t])*hi*(tmed[x, t] - ttf[x, r0/1000, t])}, tf[x, r, t], {x, 0, l}, {r, 0.0000001, (r0/1000)}, {t, 0, tempo}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 13}}];
tff[x_, r_, t_] = tf[x, r, t] /. First[enf1];

Clear[med]

(*this calculates the average temperature*)

med[x_?NumericQ, t_?NumericQ] := (NIntegrate[(mom[r]*tff[x, r, t]*r), {r, 0.0000001, r0/1000}])*(2/((r0/1000)^2*umed));
(*A plot of med[x,t] show me that the function is ok. med[x,t] came from a numeric integration*)

(*an evaluation of med[x,t] at x=0 shows that med is numeric at this position. this is valid for any time*)
med[0, tempo/2]
373.

(*so, it is created a new function that shall be numerically equal to med[x,t]*)
tmed1[x_?NumericQ, t_?NumericQ] := Evaluate[med[x, t]]

(*so, the code returns to NDSolve. The only thing that has changed here is tmed to tmed1. Why the mistake??*)

enf1 = NDSolve[{([Rho]f*cf*(D[tf[x, r, t], t] + mom[r]*D[tf[x, r, t], x])) == ((1/r)*D[(r*([Alpha]f + [CurlyEpsilon]m[r])*D[tf[x, r, t], r]), r]), tf[x, r, 0] == (tme - tinf)*Exp[-1000*x] + tinf, ttf[0, r, t] == tme, (D[tf[x, r, t], r] /. r -> 0.0000001) == 0, -kf*(D[tf[x, r, t], r] /. r -> (r0/1000)) == (1 - Exp[-1000*x])*hi*(tmed1[x, t] - ttf[x, r0/1000, t])}, tf[x, r, t], {x, 0, l}, {r, 0.0000001, (r0/1000)}, {t, 0, tempo}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 13}}];
tff[x_, r_, t_] = tf[x, r, t] /. First[enf1];
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>
5 Replies
Posted 10 years ago

In your code

In[52]:= ?ttf

Global`ttf

ttf[x_,r_,t_]=298

you have defined your ttf function to have the value 298, no matter what the values of x, r or t have.

Then in your NDSolve you use your ttf function in your system of equations

enf1 = NDSolve[{... ttf[0, r, t] == tme, ...}, ...]

and 298==373 will always be false.

I do not think this is a question of using or not using NumericQ. I think this is a misunderstanding of what you want your ttf function to be. I cannot guess from what you have shown what it is that you are trying to do or what the solution should be.

POSTED BY: Bill Simpson

You are right. It should be tf, not ttf. I'll correct this. Thanks!

Doing

?NumericQ[tmed1]

, the output is

False

What is the smart way to solve this? I don't know.

Posted 10 years ago

If I scrape his code, paste it into a notebook, restore all the missing \ in front of his [Rho], [Alpha], etc, (which I've noticed get eaten by the posting software if I attempt to do any editing of a previous post, it would be very very nice if the editing process would not do that), check the results of each step for any obvious warning signs and peek at the system that he will give to NDSolve to find enf1 the second time, just before I actually try that NDSolve, then I notice one of his items in that system is

In[50]:= ttf[0, r, t] == tme

Out[50]= False

Giving NDSolve a system that includes False as one of the items would be the first thing I would look at if I were trying to figure out why it wasn't working.

POSTED BY: Bill Simpson

You should try isolate your problem to a minimum amount of code. People usually do not have time to go through that much of code with so many definitions.

POSTED BY: Sam Carrettie
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract