Defining time every step in the loop

Answer 1

[Alan Stevens]

Are you meant to use one of Matlab's ODE solvers? If so, the following might do:

%% Times

tInitial = 0; %initial time [s]
tFinal = 1000; %final time [s] 
tSteps = 100;
timeStepArray = tInitial:(tFinal-tInitial)/tSteps:tFinal;

%% Initialising

CAD_init = 65; %contact angle [degrees]
CAR_init = CAD_init*pi/180; %contact angle [radians]
IV = 1.47e-9; %initial droplet volume[m^3]
WR = (((6*IV/pi)/...
        tan(CAR_init/2))/...
        (3 + (tan(CAR_init/2))^2))^(1/3); %initial wetted radius [m]
   
%% Calculation

F0 = [IV, CAD_init];

[t, F] = ode45(@lossrate,timeStepArray,F0);

V = F(:,1);
CAD = F(:,2);
Rho = 1000; %density [kg/m^3]
M = V*Rho;

%% Plot results
subplot(3,1,1)
plot(t,M),grid,legend('M')
subplot(3,1,2)
plot(t,V),grid,legend('V')
subplot(3,1,3)
plot(t,CAD),grid,legend('CAD')

%% Function

function dFdt = lossrate(~,F)

        RH = 0.25; %relative humidity 
        LOBF = 0.093; %average contact angle reduction as a function of time from the best fit line
        Rho = 1000; %density [kg/m^3]
        T = 25; %temperature [celsius] 
        D_T = 2.5e-4*exp(-684.15/(T+273.15)); %diffuson coefficient [m^2/s]
        c_sat = (9.99e-7)*T^3 - (6.94e-5)*T^2 + (3.2e-3)*T - 2.87e-2; %saturation concentration[kg/m^3]
     
        V = F(1);
        CAD = F(2);
        
        CAR = CAD*pi/180;
        WR = (((6*V/pi)/...
                tan(CAR)/2))/...
                (3 + (tan(CAR/2))^2)^(1/3);
        
        Vdot = (-pi*WR*D_T*(1 - RH)*c_sat*(0.27*CAR^2+1.30))/Rho;
        CADdot = -LOBF;
       
        dFdt = [Vdot; CADdot];

end

[Alan Stevens]

Are you sure the values for LOBFB and LOBFC are correct? They give rise to a huge change in value for the new contact angle at each time step!

[Alan Stevens]

Ok. The following works. Can't say if the results are sensible!

%% Input Variables

tInitial = 0; %initial time [s]
tFinal = 10000; %final time [s] 
dt = 1;  % [s]

IV = 1.47e-9; %initial droplet volume[m^3]
Rho = 1000; %density [kg/m^3]
T = 23.5; %temperature [celsius] 

%The line of best fit constants are calculated using the function file in
%the regression analysis

LOBFA = -3.3596e-5; %line of best fit constant A
LOBFB = -965.95554; %line of best fit constant B
LOBFC = 96.7624; %line of best fit constant C

CAD_init = 65; %contact angle [degrees]

RH = 0.25; %relative humidity 

%% Initialising

D_T = 2.5e-4*exp(-684.15/(T+273.15)); %diffuson coefficient [m^2/s]
c_sat = (9.99e-7)*T^3 - (6.94e-5)*T^2 + (3.2e-3)*T - 2.87e-2; %saturation concentration[kg/m^3]
CAR_init = CAD_init*pi/180; %contact angle [radians]
IWR = (((6*IV/pi)/...
        tan(CAR_init/2))/...
        (3 + (tan(CAR_init/2))^2))^(1/3); %initial wetted radius [m]
    
%% Calculation
t(1) = 0;
CAD(1) = CAD_init; %initial contact angle [degrees]
CAR(1) = CAR_init; %initial contact angle [radians]
WR(1) = IWR;
V(1) = IV; 
M(1) = IV*Rho;
i = 1;
while V(i) >= 0
            t(i+1) = i*dt;
            
            M_dot(i) = -pi*WR(i)*D_T*(1 - RH)*c_sat*(0.27*CAR(i)^2+1.30); %mass flow rate [kg/s]
            M(i+1) = M(i) + M_dot(i)*dt; %mass loss at each time step [kg]
            V(i+1) = M(i+1)/Rho;
            
            CAD(i+1) = max((LOBFA*(t(i+1)-LOBFB)^2+LOBFC),0); % new contact angle [degrees]
            CAR(i+1)= CAD(i+1)*pi/180; % new contact angle [radians]
            
            WR(i+1) = (((6*V(i+1)/pi)/...
                tan(CAR(i+1)/2))/...
                (3 + (tan(CAR(i+1)/2))^2))^(1/3); %new wetted radius [m]
            
            i = i + 1; %counter which increases for every loop
        
end
i = 1:i-1;
plot(t(i),CAD(i)),grid, legend('CAD')
figure(2)
plot(t(i),V(i)),grid, legend('V')
figure(3)
plot(t(i),M(i)),grid, legend('M')

[Alan Stevens]

max(A, 0) uses the maximum value of A or 0, So, if A goes negative zero is used.

If you don't do that here, the moment CAD goes negative then V becomes complex. I didn't look in detail to see why.

WR stays virtually constant most of the time until it too goes complex and crazy!

This suggests not all of your other equations are correct, or that the simple Euler approach to integration here is causing an unacceptable build up of errors.

[AB]

Thank you for your assistance