# Problem with NDSolve

Posted 9 years ago
3166 Views
|
5 Replies
|
3 Total Likes
|
 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
Sort By:
Posted 9 years ago
 In your code In:= ?ttf Globalttf 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 9 years ago
 You are right. It should be tf, not ttf. I'll correct this. Thanks!
Posted 9 years ago
 Doing ?NumericQ[tmed1] , the output is False What is the smart way to solve this? I don't know.
Posted 9 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:= ttf[0, r, t] == tme Out= 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 9 years ago
 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.