Problem-Solving Example #4
Liquid Propellant Rocket Problem
____________
The next Fortran Calculus example illustrates solving implicit differential equations. These equations in Figure 3 are characteristic of a liquid propellant rocket feed system in the presence of a longitudinal vibration that is damping with time.
The solution of these equations is achieved by nesting the implicit equation solver AJAX inside the ODE solver JANUS. Nesting solvers is possible to many levels as long as there is room in ones computer for temporary storage of each 2xNxN matrix where N is the number of variables in the given FIND statement; N=2 in this 'ideq' model.
Figure 3. Calculus Code ... Liquid Propellant Rocket
problem impdes
common xdot, x, ydot, y, t
x=14000 : y=7000 ! initial conditions
xdot=-50 : ydot=-25 ! initial rate guesses
t=0 : dt=.25 : tp=.5
Initiate janus; For ideq; Equations
& xdot / x, ydot / y; Of t; Step dt; To tf;
print *, ' time xdot x'
& , ' ydot y' : tf=tp
do while (tf.le.50)
Integrate ideq; By janus
print '(5e15.6)', t, xdot, x, ydot, y : tf=tf + tp
end do
end
model ideq ! implicit differential equations
common xdot, x, ydot, y, t
Find xdot, ydot; In irate(gx, gy);
& by ajax( acon); to match gx, gy
end
model irate(gx, gy) ! implicit rate equations
common xdot, x, ydot, y, t
gx=xdot + 3.2 * sqrt(1-(xdot + ydot) * exp(-t / 50)
& * (1.15 + 57.5 / (20000 + x + y)))
& * (1 + .1 * exp(-t / 10) * sin(1.5708 * t))
gy=ydot + 1.59 * sqrt(1-(xdot + ydot) * exp(-t / 50)
& * (1.15 + 36.2 / (20000 + x + y)))
& * (1 + .1 * exp(-t / 10) * sin(1.5708 * t))
end
controller acon( ajax)
summary=0
end
Initiate, Integrate, and Find statements are additions to algebra level languages
and allow this implicit system of ODEs to be solved. Initiate states solver
& model names, variable names and their interconnections. In this problem
the solver is Janus; model is IDEQ; and, variables are Xdot / X and Ydot /
Y. This tells the solver that Xdot is the derivative of X with respect to
T, the independent variable appearing in the OF phase. Same relationship exists
for Ydot and Y. Integrate statement integrates equations in model IDEQ from
initial value of T to final value, TF. The Find statement is necessary to
find Xdot & Ydot values such that Gx & Gy equations in IRATE are zero.
Higher order ODEs are solved with the same layout.
Figure 4. (Partial) Output Report ... Liquid Propellant Rocket
TIME XDOT X YDOT Y
0.500000E+00 -0.204710E+02 0.139903E+05 -0.101693E+02 0.699518E+04
0.100000E+01 -0.211458E+02 0.139798E+05 -0.105045E+02 0.698997E+04
0.150000E+01 -0.198543E+02 0.139695E+05 -0.986297E+01 0.698484E+04
0.200000E+01 -0.175393E+02 0.139601E+05 -0.871293E+01 0.698019E+04
o
o
o
0.495000E+02 -0.786710E+01 0.133984E+05 -0.390821E+01 0.670115E+04
0.500000E+02 -0.780443E+01 0.133945E+05 -0.387708E+01 0.669920E+04
ELAPSED TIME = 6.48 SECONDS