Method of characteristics Implicit method

I am trying to solve a system of four equations simultaneously. The original equations are partial differential equations as shown below:

e*dCa_RK/dt + u*dCa_RK/dz = -4*pi*(Rs^2)*Np*Ca_RK*kdash (1)

dXb/dt = ((3*Mb*b)/(Rs*rho_b*a))*Ca_RK*kdash (2)

e*dTf_RK/dt + u*dTf_RK/dz = (4*pi*(Rs^2)*Np*h*(Ts - Tf_RK))/Cf (3)

Rs*Cs*dTs/dt = -3*h*(Ts - Tf_RK) - 3*deltaH*Ca_RK*kdash (4)

where

kdash = 1./(((1/kto)*(exp(Er./(R*(Tf_RK))))) + ((Rs/Deo)*(exp(Ede./(R*(Tf_RK)))).*(((1 - (Xb)).^(-1/3)) - 1)) + ((1/ko)*(exp(Ea./(R*(Ts)))).*((1 - (Xb)).^(-2/3))))

and the above differentials are partial.

Using the method of characteristics, a new time variable is introduced: theta = t - e*(z/u) (5) therefore dtheta/dz = -(e/u) and dtheta/dt = 1 (6)

using (6) and substituting into (1) to (4) gives the following ODEs:

u*dCa_RK/dz = -4*pi*(Rs^2)*Np*Ca_RK*kdash (7)

dXb/dtheta = ((3*Mb*b)/(Rs*rho_b*a))*Ca_RK*kdash (8)

u*dTf_RK/dz = (4*pi*(Rs^2)*Np*h*(Ts - Tf_RK))/Cf (9)

dTs/dtheta = (-3*h*(Ts - Tf_RK))/(Rs*Cs) - (3*deltaH*Ca_RK*kdash)/(Rs*Cs) (10)

where

kdash = 1./(((1/kto)*(exp(Er./(R*(Tf_RK))))) + ((Rs/Deo)*(exp(Ede./(R*(Tf_RK)))).*(((1 - (Xb)).^(-1/3)) - 1)) + ((1/ko)*(exp(Ea./(R*(Ts)))).*((1 - (Xb)).^(-2/3)))).

Equations 7 - 10 are stiff ODEs to be solved implicitly.

I have tried to solve the equations using the code below but it takes too long (25 days) to run the full time required. There might be some in built function to make it faster or a better written up code and would appreciate any help.

the code is as below: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc clear all

%Length of reactor step dz = 0.001;

%Time step dtheta = 0.1; %Local time step

%Total length of reactor z = 0:dz:10; %Length of reactor in m

%Parameters of fluid A and solid B in the reactor a = 2; %Stoichiometric coefficient of fluid A b = 1; %Stoichiometric coefficient of solid B rho_b = 2960; %Density of solid B in kg/m^3 Rs = 1*10^-4; %Initial radius of solid B in m Cao = 1; %Inlet concentration of fluid A in kmol/m^3 Xbo = 0; %Initial conversion of solid B Mb = 84.3; %Molecular mass of solid B e = 0.4; %Void fraction u = 0.25; %Superficial velocity in m/s kto = 4.142*10^-3; %Fluid to solid mass transfer coefficient in m/s Deo = 4.142*10^-8; %Effective diffusion coefficient of fluid A in m/s ko = 1.19*10^11; %First order surface rate constant in m/s Np = (3*(1 - e))/(4*pi*(Rs^3)); %Number of particles suspended in a unit volume of reactor To = 316.7; %Initial bed temperature in K Er = 1*10^7; %Activation energy in J/kmol Ea = 1*10^8; %Activation energy in J/kmol Ede = 1*10^7; %Activation energy in J/kmol for Deo R = 8314; %Universal rate constant in J/kmolK h = 24000; %Specific heat of enthalpy in W/mK Cf = 4.18*10^6; %Specific heat capacity of fluid in J/Km^3 Cs = 2.977*10^6; %Specific heat capacity of solid in J/Km^3 deltaH = -52.2*10^6; %Enthalpy of reaction in J/kmol Vr = (b*Cao*u*Mb)/(a*(1 - e)*rho_b); %Reaction zone velocity in m/s Vh = u/(1 + ((1 - e))*(Cs/Cf)); %Heat zone velocity in m/s Ti = To + (((rho_b/Mb)*(a/b)*(-deltaH))/((((rho_b/Mb)*(a/b)*(1/Cao)) - Cs/Cf)*Cf)); %Intermediate temperature in K

%Initial guess for Tf, Ts, Xb and Ca for concentration of fluid A implicit solver DummyTsxvec = To*ones(1,length(z)); DummyXbxvec = Xbo*ones(1,length(z)); DummyTf_RKxvec = To*ones(1,length(z)); DummyCa_RKvec = zeros(1,length(z));

%Implicit Euler for Ca 1st iteration for b1 = 1:length(z) - 1 if b1 == 1 GuessCa_RK = Cao; DummyCa_RKvec(b1) = GuessCa_RK; else GuessCa_RK = (DummyCa_RKvec(b1)); end DummyTsx = (DummyTsxvec(b1)); DummyXb = (DummyXbxvec(b1)); DummyTff_RK = (DummyTf_RKxvec(b1)); DummyCa_RK = fsolve(@(DummyCa_RK) FUNNI4(DummyCa_RK,DummyTsx,DummyXb,DummyTff_RK,GuessCa_RK),GuessCa_RK); DummyCa_RKvec(b1 + 1) = (DummyCa_RK); end

%Store Tf values for Tf matrix Ca_RK = (DummyCa_RKvec);

%Matrix to store Ca 1st iteration S = zeros(1,length(Ca_RK)); S(1,:) = (Ca_RK);

%Initial guess for Ts and Tf for Tf 1st iteration DummyTsxvec1 = To*ones(1,length(z)); DummyTf_RKxvec1 = zeros(1,length(z)); %DummyTf_RKvec = To*ones(1,length(z));

%Implicit solver for Tf 1st iteration for k = 1:length(z) - 1 if k == 1 GuessTf_RK = To; DummyTf_RKxvec1(k) = GuessTf_RK; else GuessTf_RK = (DummyTf_RKxvec1(k)); end DummyTsxx = (DummyTsxvec1(k)); DummyTf_RK = fsolve(@(DummyTf_RK) FUNNI3(DummyTsxx,DummyTf_RK,GuessTf_RK),GuessTf_RK); DummyTf_RKxvec1(k + 1) = (DummyTf_RK); end

%Store Tf values for Tf matrix Tf_RK = (DummyTf_RKxvec1);

%Matrix to store Tf 1st iteration L = zeros(1,length(Tf_RK)); L(1,:) = Tf_RK; Tf_RKK = L(1,:);

%Implicit solver for Xb and Ts 1st iteration %Initial guess for Xb and Ts GuessXbTs = zeros(length(z),2); GuessXbTs(:,1) = Xbo*ones(1,length(z)); GuessXbTs(:,2) = To*ones(1,length(z)); GuessXbTs = GuessXbTs';

%Implicit solver for Xb and Ts 1st iteration DummyC = (Ca_RK); DummyTf_RKK = (Tf_RKK); DummyXbTs = fsolve(@(DummyXbTs) FUNNI2(DummyXbTs(2,:),DummyTf_RKK,DummyXbTs(1,:),DummyC,GuessXbTs(1,:),GuessXbTs(2,:)),GuessXbTs); Xb = (DummyXbTs(1,:)); Ts = (DummyXbTs(2,:));

%Matrix to store Xb 1st iteration M = zeros(1,length(z)); M(2,:) = (Xb);

%Matrix to store Ts 1st iteration N = zeros(1,length(z)); N(1,:) = To*ones(1,length(z)); N(2,:) = (Ts);

%Setting up Ts, Xb and Tf values for implicit solver Tf_old = L(1,:); Tsx_old = N(2,:); Xb_old = M(2,:);

%Store Ts, Tf, Ca and Xb values for Ca 2nd iteration DummyTsxvec1 = (Tsx_old); DummyXbxvec1 = (Xb_old); DummyTf_RKxvec1 = (Tf_old); DummyCa_RKvec1 = zeros(1,length(z));

%Implicit solver for Ca 2nd iteration for b2 = 1:length(z) - 1 if b2 == 1 GuessCa_RK1 = Cao; DummyCa_RKvec1(b2) = GuessCa_RK1; else GuessCa_RK1 = (DummyCa_RKvec1(b2)); end DummyTsx1 = (DummyTsxvec1(b2)); DummyXb1 = (DummyXbxvec1(b2)); DummyTf_RK1 = (DummyTf_RKxvec1(b2)); DummyCa_RK1 = fsolve(@(DummyCa_RK1) FUNNI4(DummyCa_RK1,DummyTsx1,DummyXb1,DummyTf_RK1,GuessCa_RK1),GuessCa_RK1); DummyCa_RKvec1(b2 + 1) = (DummyCa_RK1); end

%Store Ca value for matrix Ca_RK1 = (DummyCa_RKvec1);

%Matrix to store Ca 2nd iteration S(2,:) = (Ca_RK1);

%Store values for Tf 2nd iteration DummyTsxvec1 = (Tsx_old); DummyTf_RKvec1 = zeros(1,length(z));

%Implicit solver for Tf 2nd iteration for k1 = 1:length(z) - 1 if k1 == 1 GuessTf_RK1 = To; DummyTf_RKvec1(k1) = GuessTf_RK1; else GuessTf_RK1 = (DummyTf_RKvec1(k1)); end DummyTsx1 = (DummyTsxvec1(k1)); DummyTf_RK1 = fsolve(@(DummyTf_RK1) FUNNI3(DummyTsx1,DummyTf_RK1,GuessTf_RK1),GuessTf_RK1); DummyTf_RKvec1(k1 + 1) = (DummyTf_RK1); end

%Store Tf values for Tf matrix Tf_RK1 = (DummyTf_RKvec1);

%Matrix to store Tf values L(2,:) = Tf_RK1;

%Setting up Ca, Ts and Tf values for Xb iteration Ca_RKold = S(2,:); Tss_old = N(2,:); Tf_RKold= L(2,:);

%Implicit solver for Xb and Ts 2nd iteration %Initial guess for Xb and Ts GuessXbTs1 = zeros(length(z),2); GuessXbTs1(:,1) = Xb_old; GuessXbTs1(:,2) = Tss_old; GuessXbTs1 = GuessXbTs1';

%Implicit solver for Xb 1st iteration DummyC1 = (Ca_RKold); DummyTf_RKK = (Tf_RKold); DummyXbTs1 = fsolve(@(DummyXbTs1) FUNNI2(DummyXbTs1(2,:),DummyTf_RKK,DummyXbTs1(1,:),DummyC1,GuessXbTs1(1,:),GuessXbTs1(2,:)),GuessXbTs1); Xb1 = (DummyXbTs1(1,:)); Ts1 = (DummyXbTs1(2,:));

%Matrix to store Xb 2nd iteration M(3,:) = (Xb1); N(3,:) = (Ts1);

%Loop for temperatre, concentration and conversion profiles for jth iteration for j = 3:10000;

%Setting up Ts, Xb and Tf values for Ca jth iteration Tf_old = L(j - 1,:); Tsx_old = N(j,:); Xb_old = M(j,:);

%Store values Ts, Xb, Ca and Tf values for Ca jth iteration DummyTsxvec1 = (Tsx_old); DummyXbxvec1 = (Xb_old); DummyTf_RKxvec1 = (Tf_old); DummyCa_RKvec1 = zeros(1,length(z));

%Implicit solver for Ca jth iteration for b3 = 1:length(z) - 1 if b3 == 1 GuessCa_RK1 = Cao; DummyCa_RKvec1(b3) = GuessCa_RK1; else GuessCa_RK1 = (DummyCa_RKvec1(b3)); end DummyTsx1 = (DummyTsxvec1(b3)); DummyXb1 = (DummyXbxvec1(b3)); DummyTf_RK1 = (DummyTf_RKxvec1(b3)); DummyCa_RK1 = fsolve(@(DummyCa_RK1) FUNNI4(DummyCa_RK1,DummyTsx1,DummyXb1,DummyTf_RK1,GuessCa_RK1),GuessCa_RK1); DummyCa_RKvec1(b3 + 1) = (DummyCa_RK1); end

%Store Ca values for Tf matrix Ca_RK1 = (DummyCa_RKvec1);

%Matrix to store concentration values in the loop S(j,:) = (Ca_RK1);

%Setting up Ca values for Xb jth iteration Ca_RKold = S(j,:);

%Store values Ts and Tf values for Tf jth iteration DummyTsxxvec1 = (Tsx_old); DummyTf_RKvec1 = zeros(1,length(z));

%Implicit solver for Tf in the loop for k2 = 1:length(z) - 1 if k2 == 1 GuessTf_RK1 = To; DummyTf_RKvec1(k2) = GuessTf_RK1; else GuessTf_RK1 = (DummyTf_RKvec1(k2)); end DummyTsxx1 = (DummyTsxxvec1(k2)); DummyTf_RK1 = fsolve(@(DummyTf_RK1) FUNNI3(DummyTsxx1,DummyTf_RK1,GuessTf_RK1),GuessTf_RK1); DummyTf_RKvec1(k2 + 1) = (DummyTf_RK1); end

%Store Tf values for Tf matrix Tf_RK1 = (DummyTf_RKvec1);

%Matrix to store Tf values in the loop L(j,:) = Tf_RK1;

%Setting up Ts and Tf values for Xb jth iteration Tf_RKold = L(j,:); Tss_old = N(j,:); Xb_old = M(j,:);

%Implicit solver for Xb 1st iteration %Initial guess for Xb %GuessXb = zeros(1,length(z)); GuessXbTs1 = zeros(length(z),2); GuessXbTs1(:,1) = Xb_old; GuessXbTs1(:,2) = Tss_old; GuessXbTs1 = GuessXbTs1';

%Implicit solver for Xb 1st iteration DummyC1 = (Ca_RKold); DummyTf_RKK = (Tf_RKold); DummyXbTs1 = fsolve(@(DummyXbTs1) FUNNI2(DummyXbTs1(2,:),DummyTf_RKK,DummyXbTs1(1,:),DummyC1,GuessXbTs1(1,:),GuessXbTs1(2,:)),GuessXbTs1); Xb1 = (DummyXbTs1(1,:)); Ts1 = (DummyXbTs1(2,:));

%Matrix to store Xb 2nd iteration M(j + 1,:) = (Xb1); N(j + 1,:) = (Ts1);

end

%Plotting temperature, concentration and conversion profiles figure plotyy(z,(Ca_RK1/Cao),z,(Tf_RK1)); hold on plotyy(z,(Xb1),z,(Ts1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNNI4 %Function to solve for concentration of fluid A implicitly function FNI4 = FUNNI4(Ca_RK,Ts,Xb,Tf_RK,Ca_RKold)

%Parameters of fluid A and solid B in the reactor Rs = 1*10^-4; %Initial radius of solid B in m e = 0.4; %Void fraction dz = 0.001; u = 0.25; %Superficial velocity in m/s kto = 4.142*10^-3; %Fluid to solid mass transfer coefficient in m/s Deo = 4.142*10^-8; %Effective diffusion coefficient of fluid A in m/s ko = 1.19*10^11; %First order surface rate constant in m/s Np = (3*(1 - e))/(4*pi*(Rs^3)); %Number of particles suspended in a unit volume of reactor Er = 1*10^7; %Activation energy in J/kmol Ea = 1*10^8; %Activation energy in J/kmol Ede = 1*10^7; %Activation energy in J/kmol for Deo R = 8314; %Universal rate constant in J/kmolK kdash = 1/(((1/kto)*(exp(Er/(R*(Tf_RK))))) + ((Rs/Deo)*(exp(Ede/(R*(Tf_RK))))*(((1 - (Xb))^(-1/3)) - 1)) + ((1/ko)*(exp(Ea/(R*(Ts))))*((1 - (Xb))^(-2/3)))); %Overall reaction rate

%Function to be minimised to near zero FNI4 = (((-4*pi*(Rs^2)*(Np/(u)))*(kdash)*(Ca_RK)))*dz - (Ca_RK) + (Ca_RKold);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNNI3

%Function to solve for temperature of fluid implicitly function FNI3 = FUNNI3(Ts,Tf_RK,Tf_old)

%Parameters of fluid A and solid B in the reactor Rs = 1*10^-4; %Initial radius of solid B in m dz = 0.001; h = 24000; %Specific heat of enthalpy in W/mK Cf = 4.18*10^6; %Specific heat capacity of fluid in J/Km^3 u = 0.25; %Superficial velocity in m/s e = 0.4; %Void fraction Np = (3*(1 - e))/(4*pi*(Rs^3)); %Number of particles suspended in a unit volume of reactor

%Function to be minimised to near zero FNI3 = ((4*pi*(Rs^2)*(Np/(u)))*(h/(Cf))*((Ts) - (Tf_RK)))*dz - (Tf_RK) + (Tf_old);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNNI2

%Function to solve for conversion and temperature of solid implicitly function FNI2 = FUNNI2(Ts,Tf_RK,Xb,Ca_RK,Xb_old,Ts_old)

%Parameters of fluid A and solid B in the reactor a = 2; %Stoichiometric coefficient of fluid A b = 1; %Stoichiometric coefficient of solid B rho_b = 2960; %Density of solid B in kg/m^3 Mb = 84.3; %Molecular mass of solid B Rs = 1*10^-4; %Initial radius of solid B in m kto = 4.142*10^-3; %Fluid to solid mass transfer coefficient in m/s Deo = 4.142*10^-8; %Effective diffusion coefficient of fluid A in m/s ko = 1.19*10^11; %First order surface rate constant in m/s Er = 1*10^7; %Activation energy in J/kmol Ea = 1*10^8; %Activation energy in J/kmol Ede = 1*10^7; %Activation energy in J/kmol R = 8314; %Universal rate constant in J/kmolK h = 24000; %Specific heat of enthalpy in W/mK Cs = 2.977*10^6; %Specific heat capacity of solid in J/Km^3 deltaH = -52.2*10^6; %Enthalpy of reaction in J/kmol dtheta = 0.1; %Time step

kdash = 1./(((1/kto)*(exp(Er./(R*(Tf_RK))))) + ((Rs/Deo)*(exp(Ede./(R*(Tf_RK)))).*(((1 - (Xb)).^(-1/3)) - 1)) + ((1/ko)*(exp(Ea./(R*(Ts)))).*((1 - (Xb)).^(-2/3)))); %Overall reaction rate

%Function to be minimised to near zero FNI2(:,1) = (((((3*Mb*b)/(Rs*rho_b*a))*(kdash).*(Ca_RK)))*dtheta) - (Xb) + (Xb_old); FNI2(:,2) = ((((((-3*h)/(Rs*Cs))*((Ts) - (Tf_RK)))) - (((3*deltaH)/(Rs*Cs))*(kdash).*(Ca_RK)))*dtheta) - (Ts) + (Ts_old);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNNI4 , FUNNI3 and FUNNI2 have to be saved as separate scripts to the main body of code. dz and dtheta have been reduced to such a small value because it gives more accurate results but takes a much longer period of time. If there is anyway to speed up the code or to rewrite it in a different way to solve the equations (7) to (10) please let me know. [e-mail].

[Torsten]

If you have access to the NAG Toolbox,try

http://www.nag.co.uk/numeric/MB/manual64_24_1/html/D03/d03pff.html

It is fast and especially suited for your type of equations (in the original form (1)-(4)).

Best wishes

Torsten.

[Allamin]

Thanks Torsten, I will have a look and try and solve the equations this way.

Best Wishes, Daggash.

Answer 1

[John D'Errico]

I would bet that it is nearly impossible for anyone to make an attempt to answer you. (You might make it easier if you learn to format your code as such. That requires clicking on all of one button when you pasted it in.)

Instead, you can help yourself. After all, you are the only one who will be willing to invest the serious amount of time needed to figure out what it is you have done here, since you are the one who stands to gain from it. Note that to really dig into what you have written here would require one to really understand your model, what it means, where it comes from.

 - So learn to use the profiler tool. It will help you to find bottlenecks in your code.
 - Spend some time with a colleague (or advisor or boss), and explain your problem clearly. Don't just dump a mess of code on them and expect help. In fact, don't show them any code at all at first, just ask them how they might try to solve the problem. You may get ideas about a completely different approach.
 - Find a faster computer.