Acausal modeling of a chemical mixing process with heating using ModelingToolkit.jl

Hi, I’m starting to use ModelingToolkit.jl for chemical processes using an acausal approach. Here I’m showing an illustrative example of a simple mixing system with heating. Any feedback to make the model more efficient and general would be appreciated.

Two process streams are mixed to produce one of the feeds for our chemical reactor. After mixing, the blended stream is fed to a heating vessel before being sent to the reactor

Mathematical model

Acausal model

using ModelingToolkit, DifferentialEquations
using ModelingToolkit: t_nounits as t, D_nounits as D

# A port can be an inlet or an outlet
@connector Port begin
    c(t)    # concentration of component A
    T(t)    # temperature
    q(t), [connect = Flow] # heat
    v(t), [connect = Flow] # volumetric flow rate
end

@constants Tref = 25 # Reference temperature, °C

# Components 
# Stream n
@mtkmodel Stream_n begin 
    @parameters begin 
        c         # g/L 
        v         # L/min
        T         # C
        Cp = 1    # cal/g C
        ρ  = 1000 # g/L
    end
    @components begin
        s = Port()
    end
    @equations begin
        s.c ~ c
        s.T ~ T
        s.v ~ -v # flow is leaving
        s.q ~ -ρ*v*Cp*(T - Tref)   
    end
end

# Mixing tee (two inlets, one outlet)
@mtkmodel MixT begin    
    @parameters begin 
        Cp = 1    # cal/g C
        ρ  = 1000 # g/L
    end
    @variables begin
        c(t)    # outlet concentration
        v(t)    # outlet volumetric flow rate
        T(t)    # outlet temperature
        q(t)    # outlet energy flow rate
    end
    @components begin
        in₁ = Port()    # inlet 1
        in₂ = Port()    # inlet 2
        out = Port()    # outlet
    end
    @equations begin
        out.c ~ c
        out.T ~ T
        out.v ~ -v           # flow is leaving
        out.q ~ -q     
        q ~ ρ*v*Cp*(T - Tref)      
        v ~ in₁.v + in₂.v   # overall mass balance
        v*c ~ in₁.v*in₁.c + in₂.v*in₂.c # component mass balance
        q ~ in₁.q + in₂.q   # energy balance        
    end
end

# Tank
@mtkmodel Tank begin
    @parameters begin
        V = 150     # volume, L
        Cp = 1      # cal/g C
        ρ  = 1000   # g/L
        Q  = 1.05e6 # cal/min
    end
    @variables begin
        c(t)
        v(t)
        T(t)
        q(t)
        u(t)
    end
    @components begin
        inlet = Port()
    end
    @equations begin
        inlet.v ~ v # outlet flow rate = inlet flow rate
        q ~ ρ*v*Cp*(T - Tref)  
        u ~ ρ*V*Cp*(T-Tref) # internal energy 
        # ODE system
        V*D(c) ~ inlet.v*inlet.c - inlet.v*c 
        D(u) ~ inlet.q - q + Q        
    end
end

# Connecting and simulating the process 
@mtkmodel MixProcess begin
    @components begin        
        stream1 = Stream_n(c=1,v=20,T=25)
        stream2 = Stream_n(c=4,v=10,T=55)
        mixt = MixT()
        tank = Tank(V=150,Q=1.05e6)
    end
    @equations begin 
        connect(stream1.s,mixt.in₁)
        connect(stream2.s,mixt.in₂)
        connect(mixt.out,tank.inlet)
    end
end

@mtkbuild mix_model = MixProcess()

# Define the ODE (DAE) problem
c0 = 3.0    # initial concentration in the tank
T0 = 80.0   # initial temperature of tank
u0 = [mix_model.tank.c => c0, mix_model.tank.T => T0] 
guess = [mix_model.mixt.c => c0,mix_model.mixt.T => T0]
prob  = ODEProblem(mix_model,u0,(0.0,30.0);guesses = guess)
# Solve problem
sol = solve(prob)  
using Plots
plot(sol;idxs=mix_model.tank.T)
plot(sol;idxs=mix_model.tank.c)


I’m not sure what the question is but this looks right?