Skip to main content
Mathematica

Return to Question

Corrected spelling
Link
Michael E2
Michael E2
  • 244.7k
  • 18
  • 351
  • 774

Using Solution to Differential EqutionEquation in New Function

added 7 characters in body
Source Link
OK1991
OK1991
  • 37
  • 5

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

Here is where I am currently at with the code:

I am not getting any errors but I just got get any output on the pressure plot (the last one).

Here where I am currently at with the code:

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

Here is where I am currently at with the code:

added 87 characters in body; edited title
Source Link
Αλέξανδρος Ζεγγ
Αλέξανδρος Ζεγγ
  • 9.9k
  • 3
  • 20
  • 43

Using SoltuionSolution to Differential Eqution in New Function

I am using NDSolveNDSolve 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 NDSolveNDSolve 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.

(* 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)*)
\[Gamma]γ = 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[\[Gamma]*Sqrt[γ/R]*((\[Gamma]+1γ+1)/1)^((\[Gamma]+1γ+1)/(2\[Gamma]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]

Using Soltuion to Differential Eqution in New 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.

(* 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)*)
\[Gamma] = 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[\[Gamma]/R]*((\[Gamma]+1)/1)^((\[Gamma]+1)/(2\[Gamma]-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]

Using Solution to Differential Eqution in New 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.

(* 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]
Source Link
OK1991
OK1991
  • 37
  • 5