# Plug flow with ModelingToolkit?

**URL:** <https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [June 10, 2022, 7:20am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547 "2022-06-10T07:20:29Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 10, 2022, 7:20am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/1 "2022-06-10T07:20:29Z")

</div>

Hi,

Let’s say I have the following system:

- A reactor with volume V and two ports (in & out)
- 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?

Thanks,  
DD

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [June 10, 2022, 7:14pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/2 "2022-06-10T19:14:43Z")

</div>

> [@dodoplus](#):
>
> :
> 
> - A reactor with volume V and two ports (in & out)
> - 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?
> 
> Thanks,  
> DD

as far as i know, seems doable. let me check how can this be done in MTK. i know that this seems dumb, but are you asumming 1-phase flow, right?

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 11, 2022, 8:05am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/3 "2022-06-11T08:05:03Z")

</div>

yes: 1D flow.

I can see (more or less) how to do the control part using DiffEq callbacks ([Callback Library · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/features/callback_library/)) 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…

Thanks,  
DD

---

<div class="post-metadata">

**Author:** ![vettert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vettert/32/30599_2.png) [@vettert](https://discourse.julialang.org/u/vettert)\
**Post date:** [June 11, 2022, 9:39am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/4 "2022-06-11T09:39:17Z")

</div>

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.)

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 11, 2022, 11:30am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/5 "2022-06-11T11:30:37Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 11, 2022, 3:58pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/6 "2022-06-11T15:58:54Z")

</div>

If I understand, Modelica has a built-in `spatialDistribution` operator (see [3 Operators and Expressions‣ Modelica® - A Unified Object-Oriented Language for Systems Modeling Language Specification Version 3.4](https://specification.modelica.org/v3.4/Ch3.html#spatialdistribution)).

They explain:

> 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.

This is used to implement Plug Flow ([Buildings.Fluid.FixedResistances.BaseClasses](https://simulationresearch.lbl.gov/modelica/releases/v5.1.0/help/Buildings_Fluid_FixedResistances_BaseClasses.html#Buildings.Fluid.FixedResistances.BaseClasses.PlugFlow))

How would you implement this with ModelingToolkit?

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 26, 2022, 9:20am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/7 "2022-06-26T09:20:40Z")

</div>

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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 26, 2022, 12:11pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/8 "2022-06-26T12:11:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 27, 2022, 3:49pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/9 "2022-06-27T15:49:48Z")

</div>

Anyway,  
I’m going with the following (here with temperature, which is simpler):

```julia
@connector function FlowPort(;name, p=101325.0, v=0.01, T=293.15)
    sts = @variables p(t)=p v(t)=v [connect=Flow] T(t)=T [connect=Stream]
    ODESystem(Equation[], t, sts, []; name=name)
end

function MixingVolume(; name, V=1, delta_p=100, T=293.15)
    @named port_a = FlowPort()
    @named port_b = FlowPort()
    
    ps = @parameters V=V T₀=T Δp = delta_p
    sts = @variables T(t) = T₀
    eqs = [
        port_b.p - port_a.p ~ Δp
        port_a.T ~ instream(port_a.T)
        port_b.T ~ T
        D(T) ~ port_a.v/V*(port_a.T-T)
        port_b.v ~ -port_a.v
        ]
    compose(ODESystem(eqs, t, sts, ps; name=name), [port_a, port_b])
end

function Reactor(; name, V=1, delta_p=100, n = 10, T=293.15)
    @named port_a = FlowPort()
    @named port_b = FlowPort()
    V_seg = V/n
    Δp_seg = delta_p/n
    
    segs = [MixingVolume(name=Symbol(:seg_, i), V=V_seg, delta_p=Δp_seg, T=T) for i in 1:n]
    
    csegs = [connect(segs[i-1].port_b, segs[i].port_a) for i in 2:n]
    
    eqs = vcat(crolls, [
            connect(port_a, segs[1].port_a),
            connect(port_b, segs[n].port_b)])
        
    compose(ODESystem(eqs, t, [], []; name=name), vcat([port_a, port_b], segs))
end

```

DD

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 27, 2022, 6:34pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/10 "2022-06-27T18:34:26Z")

</div>

I think there a multiple issues.

- 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.

```julia
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.

See:  
[https://www.google.com/url?sa=t&source=web&rct=j&url=https://build.openmodelica.org/Documentation/Modelica%25204.0.0/Resources/Documentation/Fluid/Stream-Connectors-Overview-Rationale.pdf&ved=2ahUKEwiD1efooc74AhWmVfEDHdW5CrwQFnoECBIQAQ&usg=AOvVaw2raTHlKBwh9fEBQLMJSmm4](https://www.google.com/url?sa=t&source=web&rct=j&url=https://build.openmodelica.org/Documentation/Modelica%25204.0.0/Resources/Documentation/Fluid/Stream-Connectors-Overview-Rationale.pdf&ved=2ahUKEwiD1efooc74AhWmVfEDHdW5CrwQFnoECBIQAQ&usg=AOvVaw2raTHlKBwh9fEBQLMJSmm4)

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [June 27, 2022, 6:44pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/11 "2022-06-27T18:44:02Z")

</div>

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)

BTW, there’s no `actualstream` in MTK

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 27, 2022, 6:48pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/12 "2022-06-27T18:48:04Z")

</div>

```julia
actualstream(a) = IfElse.ifelse(a.v > 0, instream(a.T_outflow), a.T_outflow)

```

where `a` is a connector.

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [July 3, 2022, 6:48am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/13 "2022-07-03T06:48:43Z")

</div>

Modelica has a `spatialDistribution` [operator](https://specification.modelica.org/v3.4/Ch3.html) which I think deals with this kind of discretization.  
Maybe add something similar to MTK?

Also: if `x` is temperature, how would you incorporate (conductive) heat loss to the environment assuming a pipe with heat resistance `R`?

Thanks,  
DD

---

<div class="post-metadata">

**Author:** ![dodoplus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dodoplus/32/30148_2.png) [@dodoplus](https://discourse.julialang.org/u/dodoplus)\
**Post date:** [July 3, 2022, 6:58am UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/14 "2022-07-03T06:58:55Z")

</div>

Actually, `spacialDistribution` may not be needed:  
[https://mtk.sciml.ai/stable/systems/PDESystem/](https://mtk.sciml.ai/stable/systems/PDESystem/)  
[http://methodoflines.sciml.ai/dev/](http://methodoflines.sciml.ai/dev/)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 3, 2022, 1:18pm UTC](https://discourse.julialang.org/t/plug-flow-with-modelingtoolkit/82547/15 "2022-07-03T13:18:32Z")

</div>

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