Foam drainage problem using pdepe

The foam drainage equation 𝜕𝜙𝑙 𝜕𝑡 = − 𝜕 𝜕𝑧 𝑘1𝑅 2 𝜂 𝜌𝑔𝜙𝑙 2 − 𝑘2 𝛾 2𝑅 𝜙𝑙 0.5 𝜕 𝜕𝑧 𝜙𝑙 describes the spatio-temporal evolution of liquid in a foam due to gravitation and capillary forces. Herein, 𝜙𝑙(𝑥,𝑡) is the local liquid content, 𝑡 is time, 𝑧 is the spatial coordinate in the direction of the gravitational acceleration 𝑔, 𝜂 is the dynamic liquid viscosity, 𝛾 is the interfacial tension, 𝜌 is the liquid density and 𝑘1 = 0.0066, 𝑘2 = 0.1610.5 are constants. In this task, the effect of bubble size 𝑅 on the drainage of liquid out of the foam shall be studied. You can assume the following conditions: • The height 𝐻 of the foam is constant • The bubble size of constant in the whole foam • At the foam’s bottom edge (𝑧 = 𝐻), the liquid fraction is constant 𝜙𝑙 𝑡, 𝐻 = 0.36 • At the foam’s top edge (𝑧 = 0), the liquid flux is zero • At 𝑡 = 0, assume 𝜙𝑙(0, 𝑧) = 0.36 Multiscale Process Modelling and Process Analysis 2 𝑧 = 0 𝑧 = 𝐻 𝑔 Deliverable Task 02—Foam drainage • Implement the foam drainage equation in MATLAB using the pdepe command with the respective initial and boundary conditions. Use the following parameter values – 𝜂 = 0.001 𝑃𝑎 ⋅ 𝑠, 𝛾 = 0.04 𝑁/𝑚, 𝜌 = 1000 𝑘𝑔/𝑚3 , 𝑔 = 9.81 𝑚 𝑠 2 – 𝐻 = 0.1 𝑚, Δℎ = 0.0001 𝑚, Δ𝑡 = 0.01 𝑠, 𝑇 = 100 𝑠 • Solve the model for 𝑅 = 0.0005 𝑚, 𝑅 = 0.001 𝑚, 𝑅 = 0.005 𝑚 and generate a chart for each solution showing the spatial profile of the liquid fraction at times between 𝑡 = 0 and 𝑡 = 𝑇 in steps of Δ𝜏 = 0.05𝑇 • Compute the cumulated liquid flux out of the foam and compare the time-evolution with respect for each bubble size. – Use the fact that the drained liquid 𝑉𝑙,𝐷 can be computed at any time as the difference of the liquid being initially present in the foam and the liquid in the foam at time 𝑡, thus, 𝑉𝑙,𝐷(𝑡) = 𝑉𝑙,𝐹 𝑡0 − 𝑉𝑙,𝐹(𝑡) – Use the relation 𝑉𝑙,𝐹 = 𝐴𝐻 (integral lower border z=0 upper border z=H) 𝜙𝐿 𝑑𝑧, where 𝐴𝐹 is the foam’s cross-section. The trapz command can be used to calculate the integral (compare to add-on of exercise 6). Assume 𝐴𝐹 = 0.002 m2

Answer 1

[SOUMNATH PAUL]

Hi @Elena Amy ,

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:

https://in.mathworks.com/help/matlab/math/partial-differential-equations.html

Hope it helps!

Regards,

Soumnath