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?

`ClearAll["Global`*"] 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) }] `