Solving a free-boundary PDE by mapping to a fixed domain

Specific Domains Modelling & Simulations
, ,
1

Hello everyone,

I’m trying to solve a 1-D Free Boundary PDE using Julia and it is my first time manipulating the language for such a complex problem so I apologize if this sounds trivial to most of you.

Basically, the model represents a visco-elastic annulus which is going to shrink due to a contraction potential.

#Variables and unknown functions
@parameters ξ t
@variables U(..),Ρ(..),R(..)#Reference frame
# u, ρ, R in the moving frame

# ---- Approach for the free-boundary: Mapping to a fixed domain ----
# We initally wrote the equation for u(r,t) where r is the radius from origin function and r ∈ (R(t),R_max)
# We map (r,t) → ( r-R(t) / R_max-R(t) , t ) = (ξ,t) using the map function ℳ.
# U(ξ,t) = u(ℳ^-1(r,t)) = u(r,t)
ℳ(ξ,t) = R(t) + ξ*(R_max-R(t)) # To retrieve r

#Consider f a function of (r,t) in the moving frame, we note F a function of (ξ,t), its association in the reference frame
#This gives the following expression of the derivatives in the moving frame in the reference frame
#∂f/∂r = (1/(R_max-R(t)))*(∂F/∂ξ)
#∂f/∂t = ∂F/∂t - (dR/dt)*( (1-ξ)/(R_max - R(t)) )*(∂F/∂ξ)

#Our classic system of equations where we inject this is the following

# 1.0(Conservation of momentum) au - E ∂rr(u) - η ∂rrt(u) - ∂r(ζ(ρ)) = 0 becomes
# aU - (1/(R_max-R(t))^2) ( E ∂ξξ(U) - η ( ∂ξξt(U) - R'(t) (1-ξ)/(R_max - R(t)) ∂ξξξ(U) ) ) - 1/(R_max - R(t)) ∂ξ(Ρ) ζ'(Ρ) = 0

# 1.01(Neumann boundary condition for stress on free boundary) (E∂r(u) + η∂rt(u))(R(t),t) = -ζ(ρ)(R(t),t)
# ( (E ∂ξ(U) )/(R_max - R(t)) + ( η ∂ξt(U) )/(R_max - R(t)) - ( ηR'(t) (1-ξ) ∂ξξ(U) )/( (R_max - R(t))^2 ) )(0,t) = -ζ(Ρ)(0,t)

# 1.02(Dirichlet boundary condition on fixed boundary) u(R_max,t) = 0
# U(1,t) = 0

# 2.(Conservation of mass) ∂t(ρ) = -∂r(ρu)
# ∂t(Ρ) - R'(t) (1-ξ)/(R_max - R(t)) ∂ξ(Ρ) = -1/(R_max - R(t))( ∂ξ(U)Ρ + ∂ξ(Ρ)U )

# 2.1(No flux of mass on boundaries) ∂r(ρ) = 0 on R(t) and R_max
# ∂ξ(Ρ) = 0 in (0,t) and (1,t)

# 3.(Free Boundary evolution) R'(t) = ∂t(u)(R(t),t)
# R'(t) (1 + (1-ξ)/(R_max - R(t)) )∂ξ(U)(0,t) = ∂t(U)(0,t)

# This is further completed by inital conditions on each of these quantities

# ---- Code ----
#Differential operators
Dξξ = Differential(ξ)^2
Dξξξ = Differential(ξ)^3
Dξ = Differential(ξ)
Dξt = Differential(ξ)*Differential(t)
Dξξt = Differential(ξ)^2*Differential(t)
Dt = Differential(t)

#Constants
a = 1.0 #Friction
E = 1.0 #Elasticity
η = 1.0 #Viscosity
T_max = 1.0 #Maximum Time
ρ_init = 1.0 #Initial concentration
R_max = 1.0 #Maximum Radius
R_0 = 0.1 #Initial position of the free boundary

#Function
ζ(ρ) = ρ^2 + 3*ρ + 1 # Non linear Contraction potential
ζ_prim(ρ) = 2*ρ + 3 # Contraction potential derivative

#Conservations Equation
eqs = [
        a*U(ξ,t) - (1.0/(R_max-R(t))^2)*( E*Dξξ(U(ξ,t)) - η*( Dξξt(U(ξ,t)) - Dt(R(t))*(1.0-ξ)/(R_max - R(t))*Dξξξ(U(ξ,t)) ) ) - 1.0/(R_max - R(t))*Dξ(Ρ(ξ,t))*ζ_prim(Ρ(ξ,t)) ~ 0.,
        Dt(Ρ(ξ,t)) - Dt(R(t))*(1.0-ξ)/(R_max - R(t))*Dξ(Ρ(ξ,t)) ~ -1.0/(R_max - R(t))*( Dξ(U(ξ,t))*Ρ(ξ,t) + Dξ(Ρ(ξ,t))*U(ξ,t) )
        ]

#Initial and boundary conditions
bcs = [
        # Initial Conditions
         U(ξ,0) ~ 0.,      # No initial Displacement
         Ρ(ξ,0) ~ ρ_init,  # Initial concentration
         R(0) ~ R_0,       # Initial position of interface

         # Boundary Conditions at ξ=1 (fixed boundary)
        U(1,t) ~ 0.,
        Dξ(Ρ(1,t)) ~ 0.,

        # Boundary Conditions at ξ=0 (free boundary)
        # Neumann momentum condition, evaluated at ξ=0
        (E*Dξ(U(0,t)))/(R_max - R(t)) + (η*Dξt(U(0,t)))/(R_max - R(t)) - (η*Dt(R(t))*Dξξ(U(0,t)))/((R_max - R(t))^2) ~ -ζ(Ρ(0,t)),

        # No mass flux
        Dξ(Ρ(0,t)) ~ 0.,

        # Free Boundary evolution equation
        Dt(R(t)) * (1 + (Dξ(U(0,t)))/(R_max - R(t))) ~ Dt(U(0,t))
       ]

#Domains
domains = [
ξ ∈ (0,1), #Projected reference frame
t ∈ (0.0,T_max)
]

#Declaration of the system
@named pde_system = PDESystem(eqs,bcs,domains,[ξ,t],[U,Ρ,R])


#Discretization
N = 32 #Number of points
order = 4
discretization = MOLFiniteDifference([ξ=>N], t, approx_order=order)

@time prob = discretize(pde_system,discretization)

But I’m actually encountering an error during the discretization process.
ERROR: LoadError: Sym doesn't have an operation or arguments!

I cannot find what I am doing wrong, I’ve try to debug this with different LLMs but cannot find a fix that works.

Does anyone have encountered such problems ? Any help will be greatly appreciated…

Thank you,

Naoufel