The following error is found: "Warning: Failure at t=7.068767e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t. > In ode45 (line 308)" Any suggestions as to a fix?

The two files are attached. If you would like a different method then I will certainly provide it. An increase in the F_FEED or decrease in the Dr in the main file both lead to complications with the function file and the error: "Warning: Failure at t=7.068767e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t. > In ode45 (line 308)"

I am not sure how to fix this and it means I am unable to run the simulation as it is required. Thank you!

if true
  % %CE354ASSIGNMENT19792700
%NONISOTHERMAL NONISOBARIC DESIGN OF A PFR

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global F_FEED T0 Ta P0 WEIGHT_cat R U Ar L HeatingArea VOIDAGE ROcat ROgas0 g Dcat VISC VOIDVOL TORTUISITY CONSTRICTION Sa VOIDFRAC; global F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0 C_T0;

%INPUTS: F_FEED = 20; %FLOWRATE OF THE FEED per tube(mol/s) N = 6000; %NUMBER OF TUBES FEEDIN = F_FEED*N; %TOTAL FLOWRATE OF THE FEED (mol/s) Dr = 0.09; %REACTOR DIAMETER (m) T0 = 575; %INLET TEMPERATURE OF THE REACTOR (K) P0 = 101325*2.85; %INLET PRESSURE OF THE REACTOR (Pa) WEIGHT_cat = 5; %CATALYST WEIGHT (kg) R = 8.314; %Ideal gas constant (J/molK) Ta = 575; %Heat Transfer outer temp K 300 U = 00; %Heat Transfer co-eff W/m^2K 300

%REACTOR PROPERTIES: VOIDAGE = 0.4; %BED VOIDAGE ROcat = 825*1/(1-VOIDAGE); %CATALYST BULK DENSITY (kg/m^3) ROgas0 = 0.5883; %DENSITY OF THE GAS PHASE (kg/m^3) ********************variable? http://www.peacesoftware.de/einigewerte/calc_dampf.php7 Ar = pi*(Dr^2)/4; %CROSS-SECTIONAL AREA OF REACTOR L = (WEIGHT_cat/(Ar*ROcat*(1-VOIDAGE))) %REACTOR LENGTH (m) V = Ar*L*N; %VOLUME OF REACTOR (m^3) HeatingArea = pi*L*Dr; %Heat Transfer area m^2

                                                                            % V = WEIGHT_cat*(1/825);    
                                                                            % Weight = (1-VOIDAGE)*V*ROcat;  %TOTAL WEIGHT OF CATALYST IN BED
                                                                            % u = VOLUMETRICFLOWRATEin*(1/Ar);           %GAS VELOCITY (m/s)

%PARAMETERS R = 8.314; %GAS CONSTANT (J/mol*K) g = 9.81; %GRAVITAQTIONAL ACCELERATION (m^2/s) Dcat = 0.005; %CATALYST DIAMETER (m) (SPHERICAL) ********************variable? DAVE = 0.0000000015; %CATALYST PORE DIAMETER (m) SSA = 377000; %SPECIFIC SURFACE AREA (m^2/kg) VOIDVOL = 0.29*(1/100^3)*1000; %VOID VOLUME (m^3/kg) LIFESPAN = 47300000; %CATALYST LIFETIME (s- equivalent to 1.5 years) VISC = 0.0000213891; %VISCOSITY OF THE GAS PHASE (Pa*s) ********************variable? http://www.peacesoftware.de/einigewerte/calc_dampf.php7 TORTUISITY = 3; %Taken from Fogler pg816 CONSTRICTION = 0.8; %Taken from Fogler pg816 Sa = 337000; %SPECIFIED SURFACE AREA PER KG CAT (m^2/kgcat) VOIDFRAC = 0.29/0.825; %VOID FRACTION OF CATALYST

%CALCULATIONS Feed = zeros(11,1);

Feed(1,1) = 0.014*F_FEED; %INLET FLOWRATE OF G (mol/s) Feed(2,1) = 0.634*F_FEED; %INLET FLOWRATE OF H2O (mol/s) Feed(3,1) = 0.352*F_FEED; %INLET FLOWRATE OF N2 (mol/s) Feed(8,1) = T0; %INITIAL TEMPERATURE OF FEED (K) Feed(9,1) = P0; %INITIAL PRESSURE OF FEED (Pa)

    %INITIALLY

F_G0 = Feed(1); F_H2O0 = Feed(2); F_N20 = Feed(3); F_A0 = Feed(4); F_AC0 = Feed(5); F_AA0 = Feed(6); F_FA0= Feed(7); F_T0 = F_FEED;

[W,F] = ode45('ReactorsafterDeletingR4P2HOME' , [0 WEIGHT_cat], Feed);

    %FINAL CONCENTRATIONS
F_T = F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7);       
                                                                            %VOLUMETRICFLOWRATEout = F_T*R*F(8)/P;                  
                                                                            %C_TOTAL = F_T/VOLUMETRICFLOWRATEout;

C_G = C_T0.*(F(:,1)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of G C_H2O = C_T0.*(F(:,2)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of H2O C_N2 = C_T0.*(F(:,3)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of N2 C_A = C_T0.*(F(:,4)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of A C_AC = C_T0.*(F(:,5)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of AC C_AA = C_T0.*(F(:,6)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of AA C_FA = C_T0.*(F(:,7)./F_T).*(T0./F(:,8)).*(F(:,9)./P0); %Outlet concentration of FA

                                                                            %selectivity = F(:,4) - (F(:,5)+F(:,6)+ F(:,7));
selectivity = F(:,10);

secondsinayear = 300*24*60*60; %assuming 300 days per year ThousandTonsA_peryear_nocat = (((F(end,4).*0.05606).*secondsinayear)./(1000*1000))*N %a = @(t) exp(-0.000000041*t); with_Cat = @(t) (((((F(end,4).*exp(-0.000000041*t)).*0.05606))./(1000*1000))); ThousandTonsA_peryear_cat = integral(with_Cat, 0 ,secondsinayear)*N

figure subplot(4,1,1) plot(W,F(:,9)) xlabel('W') ylabel('Pressure')

% plot(W,F(:,1),W,F(:,2),W,F(:,3)) % xlabel('W') % ylabel('F_G,F_H_2_O,F_N_2') % legend('F_G','F_H_2_O','F_N_2')

subplot(4,1,2) plot(W,F(:,4),W,F(:,5),W,F(:,6),W,F(:,7)) xlabel('W') ylabel('F_A,F_A_C,F_A_A,F_F_A') legend('F_A','F_A_C','F_A_A','F_F_A')

subplot(4,1,3) plot(W,F(:,8)) xlabel('W') ylabel('Temperature')

subplot(4,1,4) plot(W,selectivity) xlabel('W') ylabel('S') legend('S')

% t = linspace(0,25920000); % figure % plot(t,with_Cat(t)) % xlabel('time (s)') % ylabel('Flowrate out of single tube (mol/s)') % legend('Acrolein')

% t = linspace(0,38880000); % figure % plot(t,a(t)) % xlabel('time (s)') % ylabel('Catalyst Activity') % legend('1st order activity')

figure plot(W,F(:,11)) xlabel('W') ylabel('C_w_p')

end

Function:

if true
  % function [ dFdW ] = ReactorsafterDeletingR4P2HOME( W,F)
%FUNCTION TO CALCULATE RHS OF FLOWRATE DIFFERENTIAL EQUATIONS

dFdW = zeros(11,1);

global F_FEED T0 Ta P0 WEIGHT_cat R U Ar L HeatingArea VOIDAGE ROcat ROgas0 g Dcat VISC VOIDVOL TORTUISITY CONSTRICTION Sa VOIDFRAC; global F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0 C_T0;

%dFdW = [r_G] % [r_H2O] % [r_N2] % [r_A] % [r_AC] % [r_AA] % [r_FA] % [T] Temperature % [P] Pressure % [SELECTIVITY]

%INITIALLY % ONWARDS IN FUNCTION %F_G0 = F(1); % F(1)= F_G ; %F_H2O0 = F(2); % F(2)= F_H2O ; %F_N20 = F(3); % F(3)= F_N2 ; %F_A0 = F(4); % F(4)= F_A ; %F_AC0 = F(5); % F(5)= F_AC ; %F_AA0 = F(6); % F(6)= F_AA ; %F_FA0= F(7); % F(7)= F_FA ; %F_T0 = sum(F_i); %T0 = F(8); % F(8)= T ; % F(9) = P; % F(10) = selectivity;

%REACTION Hf_G = -577900; %Heat of formation: G (J/mol) Hf_H2O = -241818; %Heat of formation: H2O (J/mol) Hf_N2 = 0; %Heat of formation: N2 (J/mol) Hf_A = -65000; %Heat of formation: A (J/mol) Hf_AC = -370060; %Heat of formation: AC (J/mol) Hf_AA = -166100; %Heat of formation: AA (J/mol) Hf_FA = -108600; %Heat of formation: FA (J/mol)

%Cp FUNCTIONS Cp_G = -35.3267443 + 0.645106658*F(8) - 0.00061741*F(8)^2 + 0.000000239229*F(8)^3; Cp_H2O = 32.95266198 - 0.00391977*F(8) + 0.0000222314*F(8)^2 - 0.0000000104953*F(8)^3; Cp_N2 = 28.883 -0.00157*F(8)+ 0.00000808*F(8)^2 -0.000000002871*F(8)^3; Cp_A = -14.9617445 + 0.373240705*F(8) - 0.00034559*F(8)^2 + 0.000000129748*F(8)^3; Cp_AC = -11.7479865 + 0.402868987*F(8) - 0.00031782*F(8)^2 + 0.000000109142*F(8)^3; Cp_AA = -5.82005966 + 0.23787576*F(8) - 0.00017362*F(8)^2 + 0.0000000543415*F(8)^3; Cp_FA = 19.10519375 + 0.048680586*F(8) + 0.00000878405*F(8)^2 - 0.0000000150776*F(8)^3;

%DELTACp FUNCTIONS integraldeltaCp_rxn_1 = 2*(32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)) + (-14.9617445*(F(8)-298) + 0.5*0.373240705*(F(8)^2-298^2) - (1/3)*0.00034559*(F(8)^3-298^3) + (1/4)*0.000000129748*(F(8)^4-298^4)) - ((-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4))); integraldeltaCp_rxn_2 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (-11.7479865*(F(8)-298) + 0.5*0.402868987*(F(8)^2-298^2) - (1/3)*0.00031782*(F(8)^3-298^3) + (1/4)*0.000000109142*(F(8)^4-298^4)) - (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4))); integraldeltaCp_rxn_3 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4))); % integraldeltaCp_rxn_4 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));

%HEATS OF REACTION Hrxn_1 = 2*Hf_H2O + Hf_A - Hf_G + (integraldeltaCp_rxn_1); %Heat of reaction: reaction 1 Hrxn_2 = 2*Hf_H2O + Hf_AC - Hf_G + (integraldeltaCp_rxn_2); %Heat of reaction: reaction 2 Hrxn_3 = Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_3); %Heat of reaction: reaction 3 % Hrxn_4 = Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_4);%Heat of reaction: reaction 4

%INITIAL CONCENTRATIONS VOLUMETRICFLOWRATEin = (F_FEED*R*T0)/P0; %Initial volumetric flowrate into the reactor (m^3/s) (assuming ideal gas in) NEGLECTING P/P0 to start...

C_G0 = F_G0/VOLUMETRICFLOWRATEin; %Inital concentration of G C_H2O0 = F_H2O0/VOLUMETRICFLOWRATEin; %Inital concentration of H2O C_N20 = F_N20/VOLUMETRICFLOWRATEin; %Inital concentration of N2 C_A0 = 0; %Outlet concentration of A C_AC0 = 0; %Outlet concentration of AC C_AA0 = 0; %Outlet concentration of AA C_FA0 = 0; %Outlet concentration of FA C_T0 = F_FEED/VOLUMETRICFLOWRATEin; %Initial total concentration

%FINAL Concentrations F_T = F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7);

C_G = C_T0*(F(1)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of G C_H2O = C_T0*(F(2)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of H2O C_N2 = C_T0*(F(3)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of N2 C_A = C_T0*(F(4)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of A C_AC = C_T0*(F(5)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of AC C_AA = C_T0*(F(6)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of AA C_FA = C_T0*(F(7)/F_T)*(T0/F(8))*(F(9)/P0); %Outlet concentration of FA

%Intrinsic k-VALUES k1 = 23.3244*exp(-45007.1/(R*F(8))) %Reaction rate constant k1 (m^3/kgcat*s) k2 = 2528.8*exp(-79244.89/(R*F(8))); %Reaction rate constant k2 (m^3/kgcat*s) k3 = 2150.595*exp(-73702.78/(R*F(8))); %Reaction rate constant k3 (m^3/kgcat*s) % k4 = 0.146768*exp(-24596.97/(R*F(8))); %Reaction rate constant k4 (m^3/kgcat*s)

%REACTION RATES r_1 = -k1*C_G; %Rate of reaction 1 r_2 = -k2*C_G; %Rate of reaction 2 r_3 = -k3*C_G; %Rate of reaction 3 % r_4 = -k4*C_G; %Rate of reaction 4

r_G = -1*(k1+k2+k3)*C_G; %Rate of formation of G r_N2 = 0; %Rate of formation of N2 r_H2O = k1*(C_G)^2+k2*C_G+(k3)*C_G; %Rate of formation of H2O r_A = k1*C_G; %Rate of formation of A r_AC = k2*C_G; %Rate of formation of AC r_AA = (k3)*C_G; %Rate of formation of AA r_FA = (k3)*C_G; %Rate of formation of FA

%Differential Equations to Solve dFdW(1) = (r_G); dFdW(2) = (r_H2O); dFdW(3) = (r_N2); dFdW(4) = (r_A); dFdW(5) = (r_AC); dFdW(6) = (r_AA); dFdW(7) = (r_FA); dFdW(8) = (U*HeatingArea*(Ta-F(8))+(-r_1)*(-Hrxn_1)+(-r_2)*(-Hrxn_2)+(-r_3)*(-Hrxn_3))/(F(1)*Cp_G + F(2)*Cp_H2O + F(3)*Cp_N2 + F(4)*Cp_A + F(5)*Cp_AC + F(6)*Cp_AA + F(7)*Cp_FA); %dTdW from pg 545 in Fogler; % (U*HeatingArea*(Ta-F(8))+) G = ROgas0*((F_T*R*F(8)*(1/F(9)))/Ar); %SUPERFICIAL MASS VELOCITY (kg/(m^2*s)) B = ((G*(1-VOIDAGE))/(ROgas0*g*Dcat*(VOIDAGE^3)))*(((150*(1-VOIDAGE)*VISC)/Dcat)+1.75*G); dFdW(9) = (-B/(Ar*(1-VOIDAGE)*ROcat))*(P0/F(9))*(F(8)/T0)*(F_T/F_T0); %dP/dW

dFdW(10) = (r_1)./(r_2 + r_3); %selectivity

k1m = k1*ROcat; k2m = k2*ROcat; k3m = k3*ROcat;

Dab = 0.0000105*(F(8)/T0); %Diffusivity % https://www.gsi-net.com/en/publications/gsi-chemical-database/single/11-acrolein.html De = (Dab*VOIDFRAC*CONSTRICTION)/TORTUISITY; %Effective diffusivity

Thiele1 = (Dcat/2)*sqrt(k1m/De); %Thiele modulus for first order reaction Thiele2 = (Dcat/2)*sqrt(k2m/De); Thiele3 = (Dcat/2)*sqrt(k3m/De);

n1 = (3/(Thiele1^2))*(Thiele1*coth(Thiele1)-1); %effectiveness factor n2 = (3/(Thiele2^2))*(Thiele2*coth(Thiele2)-1); n3 = (3/(Thiele3^2))*(Thiele3*coth(Thiele3)-1);

dFdW(11) = n1*(Thiele1^2) %Cwp Weisz-Prater Criterion Cwp2 = n2*(Thiele2^2) Cwp3 = n3*(Thiele3^2)

%((r_A).*n.*ROcat.*(Dcat/2)^2)./((De*C_A)) % https://pubs.acs.org/doi/full/10.1021/acs.iecr.5b02172?src=recsys end

end

[Stuart Hobson]

Apologies for text not submitting correctly.

[Stephan]

Did you try to use ode15s or ode23s like described here?

[Stephen23]

@Stuart Hobson: delete the code from your question and attach it as a file by simply clicking the paperclip button. Putting so much code within the question text hinders understanding of the question.

Answer 1

[John D'Errico]

This is the classic signal of a stiff system. If you needed more than that, the numbers with widely varying orders of magnitude should be a signal. Oh, if you needed more of a signal to that effect, the fact that it is a real world problem, with exponentials in it should have been a clue too. Anytime you have a system with various things happening in it that will have widely different time constants, consider the possibility the system is stiff.

Use a stiff solver. ODE45 is rarely a good choice in these things, even though it is normally the tool of first resort. ode15s or ode23s should be your choice here.