values don't change

valliers.jl

# Success[0.1, 0.1, 0.1, 0.1, 0.1, 0.09999999999999999, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.09999999999999999, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.09999999999999999, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]% 
# The problem is probably the pressure tank composition equation.

# challenge Valliers 2004 sweep gas claim
# flanges are a,b,c,d etc.
# where "a" is upstream and "b" is downstream
#

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

using DifferentialEquations: solve

# Stull 1947 https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Mask=4
pw(T) = 760 * 10^(4.6543 - 1435.264 / (T + 273.15 - 64.848))


RH0 = 0.5

# air, water standard cm3
@connector Flange begin
  F(t), [connect = Flow, guess = 0, bounds = (-200, 200)]
  w(t), [connect = Stream, bounds = (0, 1), guess = RH0 * pw(20) / 760]
  P(t), [bounds = (100, 2000), guess = 760]
end

@mtkmodel Ambient begin
  @components begin
  a = Flange()
end
  @parameters begin
    P = 760.0
    T = 20.0
    RH = 0.5
  end
  @equations begin
    #a.w ~ ifelse(a.F > 0, a.w, RH * pw(T) / P)
    a.w ~ RH * pw(T) / P
    a.P ~ P
  end
end

@mtkmodel PressureTank begin
  @components begin
    a = Flange()
# do I have to split this into two flanges?
# input and output?
# or even more because I have to evaluate
# the F*w term?
# Shouldn the connect/flow property of Flange take care of this?
# but I have the case where:
#    connect(p1.a, c.a, ob.b, ot.b, t2.b)
#    p1.a.w ~ c.a.w
#    c.a.w ~ ob.b.w
#    c.a.w ~ ot.b.w
#    c.a.w ~ t2.b.w
# the w equations are not included in the connect statement,
# possibly because it somehow depends on the flow?
# src/systems/connectors.jl seems to take care of this
  end
  @parameters begin
    V = 100.0  # cm3
  end
  @variables begin
    w(t), [guess  = RH0 * pw(20) / 760]
    P(t), [guess = 760 ]
  end
  @equations begin
    D(P) ~ 760. * a.F / V
    D(w)*P ~ ifelse(a.F <= 0, 0, 760 * (a.F * a.w) / V)
    #D(w) ~ 0
    a.P ~ P
    a.w ~ ifelse(a.F <= 0, w, a.w)
  end
end

@mtkmodel Membrane begin
  @components begin
    a = Flange()
    b = Flange()
  end
  @parameters begin
    pad1 = 2000.0e-11 / 0.2e-4 * 1.4e3
    pad2 = pad1 / 6000
  end
  @variables begin
    j1(t), [guess = 0]
    j2(t), [guess = 0]
  end
  @equations begin
    pad1 * (a.w * a.P - b.P * b.w) ~ j1
    pad2 * ((1.0 - a.w) * a.P - b.P * (1.0 - b.w)) ~ j2

    j2 ~ a.F * (1 - a.w)
    j1 ~ a.F * a.w

    0 ~ a.F * a.w + b.F * a.w
    0 ~ a.F + b.F
  end
end

@mtkmodel MembraneCST begin
  @components begin
    ta = PressureTank()
    tb = PressureTank()

    membrane = Membrane()
  end

  @equations begin
    connect(ta.a, membrane.a)
    connect(ta.a, membrane.b)
  end
end

# pressure regulator / throttling valve set right
# such that the b.P pressure is the set pressure (if it can be
# supplied by the "a" side)
@mtkmodel PressureRegulator begin
  @components begin
    a = Flange()
    b = Flange()
  end
  @parameters begin
    P = 760.0
  end
  @variables begin end
  @equations begin
    0 ~ a.F + b.F
    a.w ~ b.w
    b.P ~ min(P, a.P)
  end
end

@mtkmodel PressureRegA begin
  @components begin
    a = Flange()
    b = Flange()
  end
  @parameters begin
    P = 760.0
  end
  @variables begin end
  @equations begin
    0 ~ a.F + b.F
    a.F ~ ifelse(a.P <= P, 0 , ifelse(a.P <= b.P, 0, 20))
    a.w ~ b.w
  end
end

@mtkmodel Orifice begin
  @components begin
    a = Flange()
    b = Flange()
  end
  @parameters begin
    A = 1e-8 # m3
    Cd = 0.6
  end
  @variables begin
    mdot(t), [guess = 0]
  end
  @equations begin
    mdot ~ Cd * A * copysign(sqrt(2 * 100e3 * abs(a.P - b.P) / 760),
             a.P - b.P)
    a.F ~ mdot / 1.22e-6 # converted from kg to standard cm3

    0 ~ a.F + b.F
    a.w ~ b.w
  end
end

# diaphragm pump max 400 mmHg and 33. cm3/s
# different than valliars 2004 adiabatic compression
@mtkmodel Pump begin
  @components begin
    a = Flange()
    b = Flange()
  end
  @parameters begin
    pmax = 400.0
    Fmax = 33.0
  end
  @equations begin
    0 ~ a.F + b.F
    a.w ~ b.w
    max(0, min(1, 1 - (b.P - a.P) / pmax)) * Fmax ~ a.F
  end
end

@mtkmodel Setup begin
  @components begin
    ob = Orifice() # orifice at the bottom
    ot = Orifice() # orifice at the top
    a1 = Ambient(; P=761) # ambient for the bottom
    a2 = Ambient(; P=760) # ambient for everything else
    c = PressureTank() # chamber
    m1 = Membrane()
    p1 = Pump() # feed pump
    p2 = Pump() # permeate pump
    t1 = PressureRegulator(; P=600) # sweep gas throttling
    # t2 = Orifice(; A=1e-11) # retentate throttling
    t2 = PressureRegA(; P=800) # retentate throttling
  end

  @parameters begin
    Pc = 800.0
  end

  @equations begin
    # connections to ambient
    connect(a1.a, ob.a)
    connect(a2.a, ot.a, t1.a, p2.b)

    # connections to chamber
    connect(ob.b, c.a, p1.a, ot.b, t2.b)

    # connections to membrane feed side
    connect(p1.b, t2.a, m1.a)

    # connections to membrane permeate side
    connect(t1.b, p2.a, m1.b)
  end
end


@mtkbuild m = Setup()
# initial flow is zero
prob = ODEProblem(m, [m.c.w => 0.1, m.c.P => 800],  (0.0, 10), [])
sol = solve(prob, Rodas5P(), saveat=0.1)
print(sol.retcode)
print(sol[m.c.w])