For a university project, I am trying to see if my system will have choked flow and also plot the resulting pressure spike. I set up the system below to try to model the transient response.

I am able to solve the equations TEq and mEq simultaneously if I use the mDotOutChoked equation in place of the mDotOut equation. But I want to know if it ever even reaches the choked condition, and what the steady-state pressure and temperature will be so I want to start with just using the mDotOut equation.

When I run this as it is below, I get "This computation has exceeded the time limit for your plan" followed by "{sol1} is neither a list or replacement rules nor a valid dispatch table, and so..." for all the graphs.

If I replace the pressure term with just a constant, I can get some limited success depending on the value I use. This makes me think there might be an issue with having the T[t] and P[t] terms in the square root in the mDotOut equation.

Is there some issue in {sol1} or is that issue just a result of me needing to upgrade for more compute power?


Q =  100 ;(*Heat into the shroud in Watts. Based on roughly 1350 W/m^2 from the solar simulator on one face of the shroud*)
QAmb = 0  ;(*Heat loss to ambinet. Zero for now*)
A =  2*1.935*10^-5 ;(*Area of orifice m^2  Based on 1/4 inch pipe with 0.028 inch wall thickness (two outlets)*)
h =  199000 ;(*Heat of vaporization of LN2 J/kg*)
Cp =  1039 ;(*Specific heat of nitrogen J/(kg K) *)
R =  296.8 ;(*Gas constant for nitrogen J/(kg K)*)
γ = 1.40  ;(*Specific heat ratio*)
V = 0.001 ; (*Enclosed volume m^3*)
Pe = 101000 ; (*External pressure in pa*)
ρo = 4 ;(*Approx density of nitrogen at 80K in kg/m^3. This was the lowest temp data I could find*)
tf = 300 ;(*Final time in seconds*)

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

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 of mass*)

icT = T[0] == 77 ;(*initial temp in the cavity in K*)
icm = m[0] == ρo*V ;(*initial mass of the vaporized gas. Assuming it just starts at 77k at 1atm and then adding heat*)

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

Plot[{T[t]/.sol1},{t,0,tf},PlotRange -> Automatic, ImageSize->"Large",PlotLabels->Automatic, AxesLabel -> {Time (s),Temperature (K) }]
Plot[{m[t]/.sol1},{t,0,tf},PlotRange -> Automatic, ImageSize->"Large",PlotLabels->Automatic,  AxesLabel -> {Time (s),Mass (Kg) }]
Plot[P2[t], {t, 0, tf}, PlotRange -> Automatic, ImageSize -> "Large", PlotLabels -> Automatic,  AxesLabel -> {Time (s),Pressure (Pa) }]
  • 3
    $\begingroup$ If you set your tf= 0.1 , you will find that you reach the maximum number of steps very quickly NDSolve::mxst: Maximum number of 99319 steps reached at the point t == 1.0280886376503462*^-7. `. I suspect this could be due to you not handling the complex nonlinearities of density and heat capacity for a multiphase vapor-liquid problem. The answers given here Solving Stefan's solidification problem - for the case of 3 regions, may give you some hints on how to deal with the nonlinear latent heat effects. $\endgroup$ – Tim Laska Dec 16 '20 at 5:00
  • 2
    $\begingroup$ There is a complex solution from the first step with the initial data given. If we put m[0]==.014 then we get solution in 0.014s (on my ASUS). $\endgroup$ – Alex Trounev Dec 16 '20 at 14:08
  • $\begingroup$ Ok, thank you both. Yes, I am definitely stretching some of my assumptions. I ended up playing with the initial condition some more and if I set the initial mass to be just above equilibrium condition for 180C and 1atm, I get a solution. It starts out by purging that extra gas in a fraction of a second and then seems to behave as expected... not sure. $\endgroup$ – OK1991 Dec 21 '20 at 4:38

Your Answer

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

Browse other questions tagged or ask your own question.