%thinfilm_evolution %This is the function used to solve the thin-film equations (2.25)-(2.29) %in Shemilt et al. (2022). %Input arguments are: A (initial perturbation), B (thin-film capillary %Bingham number), Ymin (regularisation parameter), Tend (end-time of %simulation), N (number of finite difference grid-points), L (domain %length), Tol (error tolerance in ode15s), t (vector of time-points for the %solution to be evaluated at). %Outputs are: solH (solution struct), H, Y, pz (pressure gradient), tauw (wall shear stress). %For all simulations: Ymin = 1e-6, L = sqrt(2)*pi, Tol = 1e-11. To produce %data in figures 2(a) and 2(b): A = 0.2, B = 0.05, Tend = 1e4, N = 200. For %Newtonian simulation in figure 2(b): A = 0.2, B = 0, Tend = 1e4, N = %200. To produce data in figure 2(c): A = 0.2, B = 0.05, Tend = 1, N = 700, %t = 1. %This function is also called in the code used to produce the data in %figures 3 and 5. function [solH, H, Y, pz, tauw] = thinfilm_evolution(A,B,Ymin,Tend,N,L,Tol,t) dz = L/(N-1); %finite difference grid spacing z = (0:dz:L)'; %finite difference grid %Initial Conditions H0 = 1-A.*cos(pi.*z./L); %--% Dz and Dzzz are 1st and 3rd derivative matrices which enforce symmetry % boundary conditions. Dz = zeros(N); Dzzz = zeros(N); for n = 2:N-1 Dz(n,n-1) = -1/(2*dz); Dz(n,n+1) = 1/(2*dz); end for n = 3:N-2 Dzzz(n,n-2:n+2)= [-0.5, 1, 0,-1,0.5]/(dz^3); end Dzzz(2,1:4) = [1, -0.5, -1, 0.5]/(dz^3); Dzzz(N-1,N-3:N) = [-0.5, 1, 0.5,-1]/(dz^3); %Dzq is a 1st derivative matrix which enforces the flux boundary condition. Dzq = zeros(N); Dzq(1,2) = 1/dz; for n = 2:N-1 Dzq(n,n-1) = -1/(2*dz); Dzq(n,n+1) = 1/(2*dz); end Dzq(N,N-1) = -1/dz; Dzzz = sparse(Dzzz); Dz = sparse(Dz); Dzq = sparse(Dzq); M = -Dz-Dzzz; whos %Sparsity matrix provides a Jacobian pattern, with '1' where there might %be a non-zero entry in the Jacobian. Sparsity = zeros(N); Sparsity(1:4,1:6) = 1; Sparsity(N-4:N,N-6:N) = 1; for n = 4:N-4 Sparsity(n,n-3:n+3) = 1; end Sparsity = sparse(Sparsity); %Set options for ode15s options = odeset('RelTol',Tol,'AbsTol',Tol, 'Vectorized', 'on', 'JPattern', Sparsity); %Time span of integration tspan = [0 Tend]; %Solve the system of ODEs derived from finite differencing the evolution %equation using ode15s solH = ode15s(@deriv,tspan,H0,options); %Evaluate the solution at the points in t to find H H = deval(solH,t); %Evaluate solution quantities, derived from H Hz = Dz*H; %1st derivative Hzzz = Dzzz*H; %3rd derivative pz = -Hz-Hzzz; %pressure gradient Y = 0*H; Y(2:N-1,:) = H(2:N-1,:) - B./abs(pz(2:N-1,:)); %Y(z,t) Y = max(Y,Ymin); %include regularisation of Y tauw = pz.*H; %wall shear stress %function defining the ODE system function dH = deriv(~,H) %Define presure gradient pz = M*H; %Define Y(z,t) with Jalaal (2016) regularisation Y = 0*H; Y(2:N-1,:) = H(2:N-1,:) - B./abs(pz(2:N-1,:)); Y = max(Y,Ymin); %Define the thin-film flux q = -(1/6).*pz.*(Y.^2).*((3.*H)-Y); %Derivative dH/dt from thin-film evolution equation dH = - Dzq*q; end end