I was reading this post online where the person mentioned that using "if statements" and "abs()" functions can have negative repercussions in MATLAB's variable-step ODE solvers (like ODE45). According to the OP, it can significantly affect time-step (requiring too low of a time step) and give poor results when the differential equations are finally integrated. I was wondering whether this is true, and if so, why. Also, how can this problem be mitigated without resorting to fix-step solvers. I've given an example code below as to what I mean:
function [Z,Y] = sauters(We,Re,rhos,nu_G,Uinj,Dinj,theta,ts,SMDs0,Uzs0,... Uts0,Vzs0,zspan,K) Y0 = [SMDs0;Uzs0;Uts0;Vzs0]; %Initial Conditions
options = odeset('RelTol',1e-7,'AbsTol',1e-7); %Tolerance Levels
[Z,Y] = ode45(@func,zspan,Y0,options); function DY = func(z,y) DY = zeros(4,1); %Calculate Local Droplet Reynolds Numbers
Rez = y(1)*abs(y(2)-y(4))*Dinj*Uinj/nu_G; Ret = y(1)*abs(y(3))*Dinj*Uinj/nu_G; %Calculate Droplet Drag Coefficient
Cdz = dragcof(Rez); Cdt = dragcof(Ret); %Calculate Total Relative Velocity and Droplet Reynolds Number
Utot = sqrt((y(2)-y(4))^2 + y(3)^2); Red = y(1)*abs(Utot)*Dinj*Uinj/nu_G; %Calculate Derivatives
%SMD
if(Red > 1) DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ... Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ... (We/Re)*K*(Red^0.5)*Utot*Utot/y(2); elseif(Red < 1) DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ... Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ... (We/Re)*K*(Red)*Utot*Utot/y(2); end %Axial Droplet Velocity
DY(2) = -(3/4)*rhos*(Cdz/y(1))*Utot*(1 - y(4)/y(2)); %Tangential Droplet Velocity
DY(3) = -(3/4)*rhos*(Cdt/y(1))*Utot*(y(3)/y(2)); %Axial Gas Velocity
DY(4) = (3/8)*((ts - ts^2)/(z^2))*(cos(theta)/(tan(theta)^2))*... (Cdz/y(1))*(Utot/y(4))*(1 - y(4)/y(2)) - y(4)/z; end end
Where the function "dragcof" is given by the following:
function Cd = dragcof(Re) if(Re <= 0.01) Cd = (0.1875) + (24.0/Re); elseif(Re > 0.01 && Re <= 260.0) Cd = (24.0/Re)*(1.0 + 0.1315*Re^(0.32 - 0.05*log10(Re))); else Cd = (24.0/Re)*(1.0 + 0.1935*Re^0.6305); end end
Best Answer