Hi,

I have developed a series of two dimension PDE equations (x,t) and I am going to solve them numerically in Maple.

However, after putting all equations, initial conditions and boundary conditions, I faced couple of errors. I would be grateful if you kindly look at the code and let me know how I can solve it.

Thank you

Here is the code:

restart: with(plots): with(LinearAlgebra):

L:=0.003: # thickness

rho_w:=997.77:

rho_s:=1419:

lambda_s:=0.46:

lambda_w:=0.644:

lambda_g:=0.026:

cp_s:=3734:

cp_w:=4183:

cp_v:=1900:

cp_a:=1005.68:

M_v:=18.016:

M_a:=28.966:

R:=8.314:

epsilon:= t -> 0.9*(1-t/10):

Cw0:=6:

Sw0:= Cw0/rho_w/epsilon(0):

pi:=3.1415:

p0:=patm:

patm:=10^5:

T0:=256:

p_air:=0.2*19724:

h_m:=0.017904:

h_T:=16.746:

T_air:=380:

Xdb:= (x,t) -> S(x,t):

Cw:= (x,t) -> rho_w*epsilon(t)*S(x,t):

Cg:= (x,t) -> rho_g(x,t)*epsilon(t)*(1-S(x,t)):

Cv:= (x,t) -> rho_v(x,t)*epsilon(t)*(1-S(x,t)):

Ca:= (x,t) -> rho_a(x,t)*epsilon(t)*(1-S(x,t)):

nw:= (x,t) -> -rho_w* k_rw(x,t)*K(t)/(mu_w(x,t))*(diff(p(x,t),x)-D_c(x,t)*diff(Cw(x,t),x) ):

ng:= (x,t) -> -rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x)):

nv:= (x,t) -> -w_v(x,t)*rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x))-binary(x,t):

na:= (x,t) -> -(1-w_v(x,t))*rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x))+binary(x,t):

M_g:= (x,t) -> M_v*w_v(x,t)+M_a*(1-w_v(x,t)):

rho_g:= (x,t) -> p(x,t)*M_g(x,t)/R/T(x,t):

rho_v:= (x,t) -> rho_g(x,t)*w_v(x,t):

rho_a:= (x,t) -> rho_g(x,t)*(1-w_v(x,t)):

binary:= (x,t) -> rho_g(x,t)*epsilon(t)*(1-S(x,t))*Deff(x,t)*diff(w_v(x,t),x):

Deff:= (x,t) -> 2.3*10^(-5)*p0/p(x,t)*(T(x,t)/T0)^1.81*(epsilon(t)*(1-S(x,t)))^(4/3):

mu_w:= (x,t) -> rho_w*exp(-19.143+1540/T(x,t)):

mu_g:= (x,t) -> 0.017/1000*(T(x,t)/273)^0.65:

p_veq:= (x,t) -> p_vsat(x,t)*a_w(x,t):

p_vsat:= (x,t) -> exp(-5800.2206/T(x,t)+1.3915-0.0486*T(x,t)+0.4176*10^(-4)*T^2-0.01445*10^(-7)*T^3+6.656*ln(T(x,t))):

a_w:= (x,t) -> exp(-0.182*Xdb(x,t)^(-0.696)+0.232*exp(-43.949*Xdb(x,t))*Xdb(x,t)^0.0411*ln(p_vsat(x,t))):

h_fg:= (x,t) -> 3167.2-2.432*T(x,t):

I_vap:= (x,t) -> M_v*K_eff*(p_veq(x,t)-p(x,t))/R/T(x,t):

K_eff:=1000:

rhocp:= (x,t) -> w_v(x,t)*epsilon(t)*(1-S(x,t))*rho_g(x,t)*cp_v+(1-w_v(x,t))*epsilon(t)*(1-S(x,t))*rho_g(x,t)*cp_a+epsilon(t)*S(x,t)*rho_w*cp_w+(1-epsilon(t))*rho_s*cp_s:

lambda:= (x,t) -> epsilon(t)*(1-S(x,t))*lambda_g+epsilon(t)*S(x,t)*lambda_w+(1-epsilon(t))*lambda_s:

ncp:= (x,t) -> nv(x,t)*cp_v+na(x,t)*cp_a+nw(x,t)*cp_w:

k_rw:= (x,t) -> S(x,t)^3:

k_rg:= (x,t) -> 1.01*exp(-10.86*S(x,t) ):

K:= t -> 10^(-10)*(1-t/10):

D_c:= (x,t) -> 10^(-8)*exp(-2.8+2*Xdb(x,t)):

PDE_m1:=diff(Cw(x,t),t)+ diff(nw(x,t),x)=-I_vap(x,t) :

PDE_m2:=diff(Ca(x,t),t)+ diff(na(x,t),x)=0 :

PDE_m3:=diff(Cg(x,t),t)+ diff(ng(x,t),x)=I_vap(x,t) :

PDE_T:=diff(rhocp(x,t),t)+ diff(ncp(x,t)*T(x,t),x)=diff( lambda(x,t)*diff(T(x,t),x),x) -h_fg(x,t)*I_vap(x,t):

IBC_S:={S(x,0)=Sw0, D[1](Cw)(0,t)=0, subs(x=L,nw(x,t))=epsilon(t)*S(L,t)*h_m/R/T(L,t)*(p_veq(L,t)-p_air) }:

IBC_p:={D[1](p)(0,t)=0,p(L,t)=patm,p(x,0)=patm}:

IBC_T:={h_T*(T(L,t)-T_air)+epsilon(t)*S(L,t)*h_m/R/T(L,t)*(p_veq(L,t)-p_air)*h_fg(L,t)=0, T(x,0)=T_air, D[1](T)(0,t)=0}:

IBC_w:={w_v(x,0)=0.0262, subs(x=L,nv(x,t))=epsilon(t)*(1-S(L,t))*h_m/R/T(L,t)*(p_veq(L,t)-p_air), D[1](w_v)(0,t)=0 }:

pds:=pdsolve(PDE_m1,PDE_m2,PDE_m3,PDE_T,{IBC_S,IBC_p,IBC_T,IBC_w},numeric);