0
$\begingroup$

I am trying to feed the solutions from my differential equations into another function.

I am using NDSolve to solve two differential equations simultaneously. It seems to be working ok based on the plots it is outputting (I need to work on the physical interpretation, but at this point, I am just trying to get this to output pressure as a function of time on the plot).

Now I want to take the result of the NDSolve and plug it back into the pressure equation to see how pressure changes over time. I assumed I could just put in the interpolated numerical solution directly into the equation and then plot it. I have tried a lot of different stuff but I am new to Mathematica and am probably doing a bunch of stuff wrong.

I am not getting any errors but I am just not getting any output on the pressure plot (the last plot in the code).

Note that the pressure equation P2 (the one I am trying to evaluate using diff eq sol's) is actually used as P to solve the differential equations in the first place (i.e. P2 = P).

Here is where I am currently at with the code:

(* This calculation is intended to check the pressure build up inside the LN2 cavity*)

ClearAll["Global`*"]

Q =  100 ;(*Heat into the shroud W/m^2*)
QAmb = 0  ;(*Heat loss to ambinet*)
A =  1.935*10^-5 ;(*Area of orifice m^2  Based on 1/4 inch pipe with 0.028 inch wall thickness*)
h =  199 ;(*Heat of vaporization of LN2 kJ/kg*)
Cp =  1.039 ;(*Specific heat of LN2 kJ/(kg K) *)
R =  0.2968 ;(*Gas constant for nitrogen kJ/(kg K)*)
γ = 1.40  ;(*Specific heat ratio*)
V = 0.001 ; (*Enclosed volume m^3*)
icT = T[0] == 300 ;(*initial temp in the cavity*)
icT2 = T'[3600] == 0 ;(*After a long time we will reach steady state and T will stop changing*)
icm = m[0] == 0.001 ;(*initial mass of the vaporized gas*)
tf = 10 ;(*Final time in seconds*)

P = m[t]*R*T[t]/V ; (*Pressure term*)
mDotEvap = Q/h ; (*rate of evap*)
mDotOut = (A*P/Sqrt[T[t]])*Sqrt[γ/R]*((γ+1)/1)^((γ+1)/(2γ-2)) ; (*mass flow out of the orifice*)

TEq = T'[t] == 1/(m[t]*Cp)(mDotEvap*h-mDotOut*Cp*T[t]-QAmb) ; (*Diff Eq for Temperature in the cavity*)
mEq = m'[t] == mDotEvap - mDotOut ; (*Conservation*)

sol1 = NDSolve[{TEq,mEq, icT, icm}, {T[t], m[t]},{t, 0, tf} ] ;
 P2[t] = ( m[t]/.sol1)*R*(T[t]/.sol1)/V  ; (*Plugging back to get shroud pressure as functon of time*)


Plot[{T[t]/.sol1},{t,0,tf}, 
PlotRange -> {{0,tf}, {800, 0}}, 
ImageSize->"Large",
PlotLabels->Automatic]

Plot[{m[t]/.sol1},{t,0,tf}, 
PlotRange -> {{0,tf}, {1, 0}}, 
ImageSize->"Large",
PlotLabels->Automatic]

Plot[P2[t],{t,0,tf}, 
PlotRange -> Automatic, 
ImageSize->"Large",
PlotLabels->Automatic]
$\endgroup$
2
  • $\begingroup$ I think NDSolveValue is almost always easier to use than NDSolve since it takes out the whole messing around with replacements. $\endgroup$ – Sjoerd Smit Nov 18 '20 at 9:52
  • $\begingroup$ @SjoerdSmit, thanks for the suggestion. I did skim the NDSolveValue documentation but I did not immediately understand its benefit over NDSolve. But sounds like I should dig a little deeper to understand it better. Thanks! $\endgroup$ – OK1991 Nov 18 '20 at 18:53
2
$\begingroup$
P2[t_] = m[t]*R*T[t]/V /. sol1;
Plot[P2[t], {t, 0, tf}, PlotRange -> Automatic, ImageSize -> "Large", 
 PlotLabels -> Automatic]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.