Defining time every step in the loop
显示 更早的评论
采纳的回答
Alan Stevens
2020-7-23
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
4 个评论
Alan Stevens
2020-7-23
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
2020-7-23
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
2020-7-23
编辑：Alan Stevens
2020-7-23
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
2020-7-23
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Tuning, Analysis, and Validation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
