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)