The simplest way to simulate the behavior of a nuclear reactor is through a linearized lumped model, here are shown the basic equations.
In a BWR most of the volume is a mixture of water and steam, for this reason i assumed a constant coolant temperature equal to saturation temperature, moreover the inlet water is in saturation condition.
For neutronic equations i used the simple 1 group model (same equations used for the PWR):
\frac{dP}{dt} = \frac{\rho - \beta}{\Lambda}P +\lambda C
\frac{dC}{dt} = \frac{\beta}{\Lambda}P - \lambda C
For the fuel we have:
c_fM_f\frac{dT_f}{dt} = P - K(T_f - T_{sat})Moreover we can use this relation between fuel-coolant heat exchange and the fraction of steam x exiting the core: K(T_f - T_{sat}) = \Gamma_{total} x (h_g-h_f) where hg and hf are entalpies of steam and liquid in saturation condition. \Gamma_{total} is the total flow rate in the core.
We need x because \rho depends on the void fraction $\alpha$, that is:
\alpha =\frac{1}{x_{out}} \int_0^{x_{out}}\frac{1}{1+\eta\frac{1-x}{x}}dxWith \eta = \frac{\rho_g}{\rho_f}
Now we have to linearize \alpha and we obtain this set of equations:
\frac{d\delta P}{dt} = \frac{\delta \rho}{\Lambda} P^0 -\frac{\beta}{\Lambda}\delta P +\lambda \delta C
\frac{d\delta C}{dt} = \frac{\beta}{\Lambda}\delta P - \lambda \delta C
\delta \rho = (\alpha_f \frac{\Gamma(h_g-h_f)}{K} + \gamma s) \delta x + \alpha_h \delta h s = -\frac{\eta}{1-\eta}\frac{\frac{x_0}{x_0 + \frac{\eta}{1-\eta}} - ln(\frac{x_0 + \frac{\eta}{1-\eta}}{\frac{\eta}{1-\eta}})}{x_0^2}
x0 is the steam flow rate percentage with P0 power.
In the graph you can see the response to a step change in reactivity due to control rods.
download code (matlab) (you need XSteam)