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.`. >>
