# Question:Solving a combined system of differential and partially non-linear equations with dsolve/solve

## Question:Solving a combined system of differential and partially non-linear equations with dsolve/solve

Maple 18

Dear Maple enthusiasts,

I am unable to find a working method to solve a system of 8 equations, of which 4 are differential equations. The system contains 8 unknown variables and the goal is to find an expression for each of these variables as a function of the time t. I have attached the code of my project at the bottom of this message.

I have tried the following:

1. Using solve/dsolve to solve all 8 equations at once. This results in Maple eating up all of my memory and never finishing its calculations.
2. First using solve to solve the 4 non-differential equations so that I get 4 out of 8 variables as a function of the 4 remaining variables. This results in an expression containing RootOf() for each of the 4 veriables I'm solving for, which prevents me from using these expressions in the 4 remaining differential equations.
3. First using dsolve to solve the differential equations, which gives once again an expression for 4 variables as a function of the 4 remaining variables. I then use solve to solve the 4 remaining equations with the new found expressions. This results in an extremely long solution for each of the variables.

The code below contains the 3rd option I tried.

Any help or suggestions would be greatly appreciated. I have been scratching my head so much that I'm getting bald and whatever I search for on google or in the Maple help, I can't find a good reference to a system of differential equations together with other equations.

 > restart:

PARK - Mixed control

Input parameters

Projected interface area (m²)

 > A_int:=0.025^2*Pi:

Temperature of the process (K)

 > T_proc:=1873:

Densities (kg/m³)

 > Rho_m:=7000: metal
 > Rho_s:=2850: slag

Masses (kg)

 > W_m:=0.5: metal
 > W_s:=0.075: slag

Mass transfer coefficients (m/s)

 > m_Al:=3*10^(-4):
 > m_Si:=3*10^(-4):
 > m_SiO2:=3*10^(-5):
 > m_Al2O3:=3*10^(-5):

Weight percentages in bulk at t=0 (%)

 > Pct_Al_b0:=0.3:
 > Pct_Si_b0:=0:
 > Pct_SiO2_b0:=5:
 > Pct_Al2O3_b0:=50:

Weight percentages in bulk at equilibrium (%)

 > Pct_Al_beq:=0.132:
 > Pct_Si_beq:=0.131:
 > Pct_SiO2_beq:=3.13:
 > Pct_Al2O3_beq:=52.12:

Weight percentages at the interface (%)

Constants

Atomic weights (g/mol)

 > AW_Al:=26.9815385:
 > AW_Si:=28.085:
 > AW_O:=15.999:
 > AW_Mg:=24.305:
 > AW_Ca:=40.078:

Molecular weights (g/mol)

 > MW_SiO2:=AW_Si+2*AW_O:
 > MW_Al2O3:=2*AW_Al+3*AW_O:
 > MW_MgO:=AW_Mg+AW_O:
 > MW_CaO:=AW_Ca+AW_O:

Gas constant (m³*Pa/[K*mol])

 > R_cst:=8.3144621:

Variables

 > with(PDEtools): declare((Pct_Al_b(t),Pct_Al_i(t),Pct_Si_b(t),Pct_Si_i(t),Pct_SiO2_b(t),Pct_SiO2_i(t),Pct_Al2O3_b(t),Pct_Al2O3_i(t))(t),prime=t):

Equations

4 rate equations

 > Rate_eq1:=diff(Pct_Al_b(t),t)=-A_int*Rho_m*m_Al/W_m*(Pct_Al_b(t)-Pct_Al_i(t));

 > Rate_eq2:=diff(Pct_Si_b(t),t)=-A_int*Rho_m*m_Si/W_m*(Pct_Si_b(t)-Pct_Si_i(t));

 > Rate_eq3:=diff(Pct_SiO2_b(t),t)=-A_int*Rho_s*m_SiO2/W_s*(Pct_SiO2_b(t)-Pct_SiO2_i(t));

 > Rate_eq4:=diff(Pct_Al2O3_b(t),t)=-A_int*Rho_s*m_Al2O3/W_s*(Pct_Al2O3_b(t)-Pct_Al2O3_i(t));

3 mass balance equations

 > Mass_eq1:=0=(Pct_Al_b(t)-Pct_Al_i(t))+4*AW_Al/(3*AW_Si)*(Pct_Si_b(t)-Pct_Si_i(t));

 > Mass_eq2:=0=(Pct_Al_b(t)-Pct_Al_i(t))+4*Rho_s*m_SiO2*W_m*AW_Al/(3*Rho_m*m_Al*W_s*MW_SiO2)*(Pct_SiO2_b(t)-Pct_SiO2_i(t));

 > Mass_eq3:=0=(Pct_Al_b(t)-Pct_Al_i(t))+2*Rho_s*m_Al2O3*W_m*AW_Al/(Rho_m*m_Al*W_s*MW_Al2O3)*(Pct_Al2O3_b(t)-Pct_Al2O3_i(t));

1 local equilibrium equation

Gibbs free energy of the reaction when all of the reactants and products are in their standard states (J/mol). Al and Si activities are in 1 wt pct standard state in liquid Fe. SiO2 and Al2O3 activities are in respect to pure solid state.

 > delta_G0:=-720680+133*T_proc:

Expression of mole fractions as a function of weight percentages (whereby MgO is not taken into account, but instead replaced by CaO ?)

 > x_Al2O3_i(t):=(Pct_Al2O3_i(t)/MW_Al2O3)/(Pct_Al2O3_i(t)/MW_Al2O3 + Pct_SiO2_i(t)/MW_SiO2 + (100-Pct_SiO2_i(t)-Pct_Al2O3_i(t))/MW_CaO); x_SiO2_i(t):=(Pct_SiO2_i(t)/MW_SiO2)/(Pct_Al2O3_i(t)/MW_Al2O3 + Pct_SiO2_i(t)/MW_SiO2 + (100-Pct_SiO2_i(t)-Pct_Al2O3_i(t))/MW_CaO);

Activity coefficients

 > Gamma_Al_Hry:=1: because very low percentage present  during the process (~Henry's law)
 > Gamma_Si_Hry:=1: because very low percentage present  during the process (~Henry's law)
 > Gamma_Al2O3_Ra:=1: temporary value!
 > Gamma_SiO2_Ra:=10^(-4.85279678314968+0.457486603678622*Pct_SiO2_b(t)); very small activity coefficient? plot(10^(-4.85279678314968+0.457486603678622*Pct_SiO2_b),Pct_SiO2_b=3..7);

Activities of components

 > a_Al_Hry:=Gamma_Al_Hry*Pct_Al_i(t); a_Si_Hry:=Gamma_Si_Hry*Pct_Si_i(t); a_Al2O3_Ra:=Gamma_Al2O3_Ra*x_Al2O3_i(t); a_SiO2_Ra:=Gamma_SiO2_Ra*x_SiO2_i(t);

Expressions for the equilibrium constant K

 > K_cst:=exp(-delta_G0/(R_cst*T_proc));
 > Equil_eq:=0=K_cst*a_Al_Hry^4*a_SiO2_Ra^3-a_Si_Hry^3*a_Al2O3_Ra^2;

Output

 > with(ListTools): dsys:=Rate_eq1,Rate_eq2,Rate_eq3,Rate_eq4: dvars:={Pct_Al2O3_b(t),Pct_SiO2_b(t),Pct_Al_b(t),Pct_Si_b(t)}: dconds:=Pct_Al2O3_b(0)=Pct_Al2O3_b0,Pct_SiO2_b(0)=Pct_SiO2_b0,Pct_Si_b(0)=Pct_Si_b0,Pct_Al_b(0)=Pct_Al_b0: dsol:=dsolve({dsys,dconds},dvars):
 > Pct_Al2O3_b(t):=rhs(select(has,dsol,Pct_Al2O3_b)[1]); Pct_Al_b(t):=rhs(select(has,dsol,Pct_Al_b)[1]); Pct_SiO2_b(t):=rhs(select(has,dsol,Pct_SiO2_b)[1]); Pct_Si_b(t):=rhs(select(has,dsol,Pct_Si_b)[1]);
 > sys:={Equil_eq,Mass_eq1,Mass_eq2,Mass_eq3}: vars:={Pct_Al2O3_i(t),Pct_SiO2_i(t),Pct_Al_i(t),Pct_Si_i(t)}: sol:=solve(sys,vars);

,