Simulating membrane Julia

Specific Domains Modelling & Simulations
1 
using Pkg
Pkg.add("ModelingToolkit")
Pkg.add("MethodOfLines")
Pkg.add("DomainSets")
Pkg.add("OrdinaryDiffEq")

using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq
include("Properties2.jl")

import ModelingToolkit: scalarize, @register_symbolic
import Symbolics: @syms

# Constants

Nc = 1

 # Define Pf as a numeric value
 # Define Pp as a numeric value


@register_symbolic μ_mix(y, Mw, μᵢ)
# Define parameters
@parameters z k=0.1 Tref=300.0 Pref=101325.0 T=298.0 Di=0.15 N=100 D=0.2 Pf_0=1 Pp_L=1 Nf=100 L = 0.6 Pp_value = 101325.0 Pf_value = 101325.0 Pf_value = 101325.0 D0 = 3.0 R = 8.314 
@parameters a_i[1:Nc] Mw_i[1:Nc]

# Define variables with updated array-valued syntax
vars1 = @variables Fₚ(..) Ff(..) Pf(..) Pp(..) 
vars2 = @variables (Fp_i(..))[1:Nc] (Ff_i(..))[1:Nc] (J_i(..))[1:Nc] (x_i(..))[1:Nc] (y_i(..))[1:Nc]
vars3 = @variables (xm_i(z))[1:Nc] (A_i(z))[1:Nc] (N_i(z))[1:Nc] (ρᵢ(z))[1:Nc] Nt(z) X(z)

callable_scalar = [Fₚ, Ff, Pf, Pp]

var_scalars = [var(z) for var in callable_scalar]

var_vectors = [collect(var(z)) for var in vars2]

variables = [var_scalars...; var_vectors...; vars3...]

# Define derivatives with respect to z
Dz = Differential(z)

# Example values for `a_i` and `Mw_i`
a_i_vals = [1.0]  # Example values for `a_i`
Mw_i_vals = [28.0]  # Example molecular weights

# Substitute values for `a_i` and `Mw_i`
param_map = Dict([(a_i[i] => a_i_vals[i]) for i in 1:Nc]...)
param_map = merge(param_map, Dict([(Mw_i[i] => Mw_i_vals[i]) for i in 1:Nc]...))
# Define the equations
eqn = [
    scalarize(J_i(z) .~ a_i .* (Pf(z) .* xm_i .- Pp(z) .* y_i(z)))...
    scalarize(N_i .~ A_i .* (Pf(z) .* xm_i - Pp(z) .* y_i(z)))...
    scalarize(N_i .~ k .* (x_i(z) - xm_i) .* X + x_i(z) .* Nt)...
    scalarize(A_i .~ a_i .* Mw_i ./ ρᵢ)...
    scalarize(ρᵢ .~ (y_i(z) .* Pp(z)) ./ (R * T))...
    Nt ~ scalarize(sum(N_i))
    X ~ ((Pf(z)) / Pref) * (Tref / T)
    #scalarize(sum(x_i(z))) ~ 1  # Constraint on sum of x_i
    #scalarize(sum(y_i(z))) ~ 1  # Constraint on sum of y_i
    #scalarize(sum(xm_i)) ~ 1  # Constraint on sum of xm_i
    scalarize(x_i(z).~ Ff_i(z)/Ff(z))...
    scalarize(y_i(z) .~ Fp_i(z)/Fₚ(z))...
    scalarize(Dz(Ff_i(z)) .~ -J_i(z) .* π .* D0 .* Nf)...  # Differential equation for Ff_i
    scalarize(Dz(Fp_i(z)) .~ J_i(z) .* π .* D0 .* Nf)...    # Differential equation for Fp_i  Counter current flow
    Dz(Pf(z)) ~ (192 * N * D0 * (D + N * D0) * R * T * μ_mix(y_i(z), Mw_i, ρᵢ) * Ff(z)) / (π * (D^2 - N * D0^2)^3 * Pf(z))  # Call μ_mix with arguments
    Dz(Pp(z)) ~ (-128 * R * T * μ_mix(y_i(z), Mw_i, ρᵢ) * Fₚ(z)) / (π * Di^4 * Nf)  # Use Fₚ instead of Fp
    Fₚ(z) ~ scalarize(sum(Fp_i(z)))
    Ff(z) ~ scalarize(sum(Ff_i(z)))
]

println("Number of equations: ", length(eqn))
println("Number of variables: ", length(variables))

# Define boundary conditions
Ix = Integral(z in DomainSets.ClosedInterval(0, L))  # Cumulative sum from 0 to L
init = 0

bcs = [
    Pf(0) ~ Pf_0,
    Pp(L) ~ Pp_L,
    scalarize([Ff_i(0)[i] ~ x_i(0)[i] * Ff(0) for i in 1:Nc])...,
    scalarize([Fp_i(0)[i] ~ Ix(J_i(z)[i]) * π * D0 * Nf for i in 1:Nc])...,
    #scalarize([ρᵢ[i] ~ x_i(0)[i] * 1.0 for i in 1:Nc])...
]

domains = [z ∈ Interval(0.0, L)]

@named pdesys = PDESystem(eqn, bcs, domains, z, variables, param_map)
# Simplify the system
simplified_sys = structural_simplify(pdesys)

# Print the simplified system equations and variables
println("Simplified system equations: ", equations(simplified_sys))
println("Simplified system variables: ", states(simplified_sys))
# Discretizing the domain using method of lines
Nz = 15
dz = L / Nz

discretization = MOLFiniteDifference([z => dz], nothing, approx_order = 4)

prob = discretize(pdesys, discretization, should_transform = false)

Properties2:

import ModelingToolkit: scalarize
#Defining functions for viscosities

#function μ(pars, T)
#    A, B, C, D = pars
#    μ= A*T^B/(1+C/T +D/T^2)
 #   return μ
#end

function μ(pars, T)
    A, B, C, D = pars
    return A * T^B / (1 + C/T + D/T^2)
end

function μ_mix(y, Mw, μᵢ)
    _1_Mw = scalarize(1.0./Mw) 
    num = scalarize(y.*μᵢ)
    sqrt_Mw_mat = scalarize(.√(collect(Mw)*collect(_1_Mw')))
    _den = scalarize(y.*sqrt_Mw_mat)
    den = scalarize(sum([_den[i, :] for i in 1:size(_den, 1)]))
    return scalarize(sum(num./den))
end

I am recieving the following error:
ERROR: MethodError: no method matching PDESystem(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Num, ::Vector{…}, ::Dict{…}, ::Dict{…}, ::Nothing, ::Vector{…}, ::Nothing, ::Nothing, ::Symbol, ::String, ::Nothing, ::Nothing; checks::Bool)

1 Like