To my understanding of this problem, we need to define the partial differential equation, set the initial and the boundary conditions. In the second step I have solved the equation for the specified bubble sizes and plotted the spatial profile of the liquid fraction over time. Finally, I have computed the cumulated liquid flux out of the foam.
function foam_drainage
% Given parameters
eta = 0.001; % Pa.s
gamma_tension = 0.04; % N/m
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
H = 0.1; % m
A_f = 0.002; % m^2 (cross-section of the foam)
k1 = 0.0066;
k2 = 0.161^0.5;
% Discretization parameters
dz = 0.0001; % m
dt = 0.01; % s
T = 100; % s
dTau = 0.05 * T; % s
% Bubble sizes to study
R_values = [0.0005, 0.001, 0.005]; % m
% Spatial mesh
zmesh = 0:dz:H;
% Time vector
tspan = 0:dt:T;
% Loop over the bubble sizes
for R = R_values
% Solve PDE
sol = pdepe(0, @(z,t,u,dudz) pdefun(z,t,u,dudz,k1,k2,R,eta,rho,g,gamma_tension), ...
@icfun, ...
@bcfun, ...
zmesh, tspan);
% Extract the solution for phi_l
phi_l = sol(:,:,1);
% Plot the spatial profile of the liquid fraction at specified times
figure;
hold on;
for t = 0:dTau:T
[~, tIdx] = min(abs(tspan - t)); % Find the closest time index
plot(zmesh, phi_l(tIdx,:), 'DisplayName', sprintf('t = %.2f s', t));
end
hold off;
title(sprintf('Spatial profile of liquid fraction over time (R = %.4f m)', R));
xlabel('Height z (m)');
ylabel('Liquid fraction \phi_l');
legend('show');
% Compute the cumulated liquid flux out of the foam
V_l_F_initial = A_f * H * trapz(zmesh, phi_l(1,:)); % Initial liquid volume in the foam
V_l_F = A_f * H * trapz(zmesh, phi_l, 2); % Liquid volume in the foam over time
V_l_D = V_l_F_initial - V_l_F; % Drained liquid volume over time
% Plot the time-evolution of the cumulated liquid flux for each bubble size
figure;
plot(tspan, V_l_D, 'DisplayName', sprintf('R = %.4f m', R));
title('Cumulated liquid flux out of the foam over time');
xlabel('Time t (s)');
ylabel('Drained liquid volume V_l_D (m^3)');
legend('show');
end
end
% Define the PDE function
function [c,f,s] = pdefun(z,t,u,dudz,k1,k2,R,eta,rho,g,gamma_tension)
c = 1;
f = -k1*R^2/(eta*rho*g) * u^2 - k2*gamma_tension^2/(2*R) * sqrt(u) * dudz;
s = 0;
end
% Define the initial condition function
function u0 = icfun(z)
u0 = 0.36; % Initial liquid fraction
end
% Define the boundary condition function
function [pl,ql,pr,qr] = bcfun(zl,ul,zr,ur,t)
pl = ul - 0.36; % At the bottom edge (z = H), the liquid fraction is constant
ql = 0;
pr = 0; % At the top edge (z = 0), the liquid flux is zero
qr = 1;
end
You can also see the examples section in below link for further reference:
Hope it helps!
Regards,
Soumnath
