Hi everyone,
I am trying to solve 2 coupled non-linear PDEs using the MOL with DifferentialEquations.jl and Finite differences first in 1D and I am stuck somewhere.
My two equations describe the evolution of a porosity wave in a two-phase flow and are, in 1D and in dimensionless form:
with De, n and m beeing constants and Φ, the porosity and Pe, the pressure, beeing my two unknows.
My code so far looks like that:
#=
Equations dimensionless:
∂ϕ/∂t = - ∂/∂z( ϕ^n (∂/∂z(Pe) + 1) )
De ∂Pe/∂t = ∂/∂z( ϕ^n (∂/∂z(Pe) + 1) ) - ϕ^m Pe
discretisation:
∂ϕ/∂t = - 1/Δx * ((ϕ_i+1/2) * (((Pe_i+1 - Pe_i) / Δx) + 1)
- (ϕ_i-1/2) * (((Pe_i - Pe_i-1) / Δx) + 1))
∂Pe/∂t = 1/De * (1/Δx * ((ϕ_i+1/2) * (((Pe_i+1 - Pe_i) / Δx) + 1)
- (ϕ_i-1/2) * (((Pe_i - Pe_i-1) / Δx) + 1))
- ϕ_i * Pe_i)
=#
using DifferentialEquations
using Plots
const nx = 501
const Lx = 150
const Δx = Lx / (nx - 1)
# Constants of the problem
const n = 3
const m = 2
const De = 1
# Initial conditions
const Phi0 = 1
const dPhi1 = 8
const dPhi2 = 1
const z1 = 10
const z2 = 40
const lambda = 1
const dPe1 = dPhi1 / De
const dPe2 = dPhi2 / De
ϕ0 = zeros(nx)
Pe0 = zeros(nx)
for i in 1:nx
ϕ0[i] =
Phi0 +
dPhi1 * exp(-((i - z1)^2) / lambda^2) +
dPhi2 * exp(-((i - z2)^2) / lambda^2)
Pe0[i] =
-dPe1 * exp(-((i - z1)^2) / lambda^2) -
dPe2 * exp(-((i - z2)^2) / lambda^2)
end
# Definition of the problem
u0 = [ϕ0 Pe0]
p = [De, n, m]
t = [0, 20]
function porosity_wave(du,u,p,t)
ϕ = @view u[:, 1]
Pe = @view u[:, 2]
De, n, m = p
@inbounds for i in 1:nx
# boundary condition with dirichlet ϕ = 1, Pe = 0 on both sides
if (i == 1)
# porosity
du[i, 1] =
- (((ϕ[i+1]^n + ϕ[i]^n) / 2 * ((Pe[i+1] - Pe[i]) / Δx + 1)
- (ϕ[i]^n + 1^n) / 2 * ((Pe[i] - 0) / Δx + 1)) / Δx)
# pressure
du[i, 2] = 1 / De * (
(((ϕ[i+1]^n + ϕ[i]^n) / 2 * ((Pe[i+1] - Pe[i]) / Δx + 1)
- (ϕ[i]^n + 1) / 2 * ((Pe[i] - 0) / Δx + 1)) / Δx)
- ϕ[i]^m * Pe[i])
elseif (i == nx)
# porosity
du[i, 1] =
- (((1 + ϕ[i]^n) / 2 * ((0 - Pe[i]) / Δx + 1)
- (ϕ[i]^n + ϕ[i-1]^n) / 2 * ((Pe[i] - Pe[i-1]) / Δx + 1)) / Δx)
# pressure
du[i, 2] = 1 / De * (
(((1 + ϕ[i]^n) / 2 * ((0 - Pe[i]) / Δx + 1)
- (ϕ[i]^n + ϕ[i-1]^n) / 2 * ((Pe[i] - Pe[i-1]) / Δx + 1)) / Δx)
- ϕ[i]^m * Pe[i])
else
# porosity
du[i, 1] =
- (((ϕ[i+1]^n + ϕ[i]^n) / 2 * ((Pe[i+1] - Pe[i]) / Δx + 1)
- (ϕ[i]^n + ϕ[i-1]^n) / 2 * ((Pe[i] - Pe[i-1]) / Δx + 1)) / Δx)
# pressure
du[i, 2] = 1 / De *
(((ϕ[i+1]^n + ϕ[i]^n) / 2 * ((Pe[i+1] - Pe[i]) / Δx + 1)
- (ϕ[i]^n + ϕ[i-1]^n) / 2 * ((Pe[i] - Pe[i-1]) / Δx + 1)) / Δx
- ϕ[i]^m * Pe[i])
end
end
end
prob = ODEProblem(porosity_wave, u0, t, p)
# solver for stiff equations
sol = solve(prob, TRBDF2(linsolve=LinSolveGMRES()))
but I don’t find the solution I expect. Here is what ϕ should look like in time:
and here is 1 timestep for example with my code:
So basically it blows up.
As I am quite new to DifferentialEquations.jl and the MOL, did I do something wrong with the way I am discretizing the equations? Or is it the wrong solver for my problem?
Thx a lot in advance,
