Two fluid reservoirs, say of oil and water (with different specific densities)
A valve selecting which of the reservoirs is connected to the “in” port of the reactor
A pump that ensures the fluids flow from the reservoir to the reactor at a constant rate (m^3/s)
A control mechanism that for our purpose at each “tick” (say every 10 seconds) selects the valve’s state (i.e. whether oil or water will flow into the reactor)
Now, the “plug flow” assumption means that if the valve switches at time 10 from oil to water, the “in” port will “see” water flowing in immediately, but until a volume of V was pumped into the reactor, the “out” port will still see oil flowing out of the reactor.
I’d like to model the specific density (i.e. water vs. oil) of the fluid exiting the reactor at the port “out”.
Is something like this possible with ModelingToolkit?
I can see (more or less) how to do the control part using DiffEq callbacks (Callback Library · DifferentialEquations.jl) which I think are available at the ModelingToolkit level.
I’m less sure about the “plug flow”: there is no simple equation that describes the state of the reactor.
in the general case, there are several “plugs” traveling down the reactor which need to be kept track off. This can be easily done via Julia code, of course…
In its simplest form this is a 1D partial differential equation (advection equation); so you would ultimately have to decide how to discretize your space domain (length along reactor) and a numerical scheme (finite diff, finite element, finite volume, etc.)
Oh, I think I misspoke: I’d like (at least for now), to model the reactor as a 0D (“lump”) component.
All it does is delay the flow: the output at time t is the same as the input at time t-v/V (where V is the volume of the reactor and v is the flow rate; obviously, if v is not constant, an integral would be needed).
Many applications involve the modelling of variable-speed transport of properties. One option to model this infinite-dimensional system is to approximate it by an ODE, but this requires a large number of state variables and might introduce either numerical diffusion or numerical oscillations. Another option is to use a built-in operator that keeps track of the spatial distribution of z(x, t), by suitable sampling, interpolation, and shifting of the stored distribution. In this case, the internal state of the operator is hidden from the ODE solver.
Seems to result in a delay differential equation (DDE), for which there is support in DifferentialEquations.jl but I think at the moment not in MTK. However you can discretize the advection equation via a method of lines approach and solve the resulting ODE. This can be done in MTK.
using ModelingToolkit, DifferentialEquations
function advection(;name, c=1, L=10, N=10, x0=0.0, xL=1)
@parameters t
sts = @variables x[1:N](t) = fill(x0, N)
Δz = L / N
Dt = Differential(t)
eqs = [
Dt(x[1]) ~ -c * (-xL + x[1]) / Δz
[Dt(x[i]) ~ -c * (-x[i-1] + x[i]) / Δz for i in 2:N]...
]
return ODESystem(eqs, t, vcat(sts...), []; name)
end
@named model = advection(N=10)
sys = structural_simplify(model)
prob = ODAEProblem(sys, Pair[], (0.0, 50.0))
sol = solve(prob, Rodas5())
However, you will get numerical diffusion and for changing velocities you have to use an upwind scheme. There are more sophisticated discretization approaches as well.
Yes, I understand that I can split the reactor into N parts. I was trying to avoid that.
At any rate, I think that a simple enough model for my problem would be to assume each “part” is perfectly mixed. I think this can be done with instream but I can’t find documentation as to how it works.
The stream connectors work in such a way that the stream variable is the value of the outflowing (from the view of the component, i.e. if the flow variable, in this case v is negative). In Modelica these are often called ..._outflow, to make that clear. So the temperature that is flowing out of the ports of the volume can only be the fluid temperature inside the volume itself.
The energy balance seems incorrect as well. The energy content from port b is missing.
eqs = [
port_b.p - port_a.p ~ Δp
port_a.T_outflow ~ T
port_b.T_outflow ~ T
D(T) ~ 1 / V * (port_a.v * actualstream(port_a.T_outflow) + port_b.v * actualstream (port_b.T_outflow) # energy balance
port_b.v + port_a.v ~ 0 # mass balance
]
Depending on what v is, a heat capacity is missing from the energy balance equation as well.
I was assuming port_a is always the “in” port.
In that case, the energy flowing out via port_b is: port_b.v/V*T == -port_a.v/V*T. (v obviously is the fluid’s velocity)
Yeah though MethodOfLines is like 3 months old so it will need some time to be complete enough for all of the cases. For this kind of case it may need WENO discretizations