I am currently working on a mathematical model to simulate the transient behaviour of a steam generator in a nuclear reactor, but I have stumbled across some issues when it comes to solving the governing equations.
My model is based on the so called Moving Boundary Models for Dynamic Simulations of Two-Phase Flows in which conservation of mass and energy equations are formulated for the primary fluid in the three different flow regions; sub-cooled liquid, two-phase flow and super-heated liquid. Added to these equations are also equations describing heat transfer in the wall separating the steam and the secondary fluid, as well as conservation of energy equations for the secondary fluid.
I end up with a set of 12 equations, 9 dependent variables and 1 independent variable (time). I will not write the complete system, but as an example, a small part of my system looks like this (conservation eqs. for super-heated region):
\begin{align} &\alpha(T_3)\frac{dL_1}{dt}+\beta(T_3)\frac{dL_2}{dt}+L_3\gamma(T_3)\frac{dh_{out}}{dt}=0\\ & \delta(T_3)h_{out}\left(\frac{dL_1}{dt}+\frac{dL_2}{dt}\right)+\epsilon(T_3)h_{out}L_3\frac{dh_{out}}{dt}=\\ &\dot{m}(h_g-h_{out})+HTC(T_3)L_3\zeta\cdot(T_{3, wall}-T_3) \end{align}
Where $\alpha,\beta,\gamma,\delta,\epsilon$ are dependent on the average temperature $T_3$ in the region (calculated from the average enthaply $h_{avg} = 0.5(h_g+h_{out})$), $L_3$ is the length of the super-heated region and it is calculated as $L_3 = L_{tot}-L_1-L_2$, where $L_{tot}$ is the total heat exchanging length and $L_1$, $L_2$ are the sub-cooled and two-phase lengths respectively. $h_{out}$ is outlet enthalpy and $HTC(T_3)$ is the convective heat transfer coefficient which in turn is calculated from its own set of equations based on the temperature $T_3$.
To summarise, we have in this small example the dependent variables $L_1,L_2,h_{out}$ and their time derivatives written in a coupled system.
With the twelve equations I can formulate a matrix equation on the following form:
\begin{equation} [A]_{12x9}\cdot[x]_{9x1}=[b]_{12x1} \end{equation} where the vector $[x]_{9x1}$ contains all derivatives as
$$[x]=\begin{bmatrix}\frac{dL_1}{dt}\\ \frac{dL_2}{dt}\\ \frac{dh_{out}}{dt}\\\vdots\\ \end{bmatrix}$$
and the matrix $[A]$ and vector $[b]$ contains constants and dependent variables.
Now that my problem has been introduced let me briefly talk about the approach I've used thus far.
My approach has been to write a RK4 solver, in later iterations it was changed into a RKF45 solver, and to recalculate all of the variables which depends on the system state eg. $\alpha(T_3)$ or $HTC(T_3)$, based on the system state in the previous iteration. However all RK/RKF solvers I've seen thus far requires the system to be on the form
$$\dot{\vec{y}}=A\vec{y}.$$
However, that's not how my system looks. My idea was then to use a least squares algorithm to solve for the derivative vector in my matrix equation and then use that as a basis for the RKF45 solver. I do suspect that this method is incorrect. I've been searching all around the internet for a potential solution to a similar problem, but almost everything assumes that the searched for time derivatives are all placed on one side of the equal sign, and that the dependent variables are multiplied with a matrix A.
My Question: What I would like to know is if anyone has some experience of solving a problem of this sort, and to confirm if my method involving least squares is incorrect. I would be happy to be pointed in the right direction if it is obvious where I'm doing something wrong. Also if anyone is aware of some good literature in this topic I would deeply appreciate it.
I hope that this long post makes sense, and that someone has the time/ will to guide me in the right direction.