%% MCHI, cv5 2015, file 1 (cv5_2015_1) - Implicit function of 1 variable
%
% Author:       Martin Isoz
% Organisation: ICT Prague
% Date:         27. 09. 2014
%
% License: This code is published under MIT License, please do not abuse
% it.
%

%% Background - Rivulet type flow on a vertical wetted plate
%
% Usually, the fluid flow is described by Navier-Stokes equation (NSE),
% which is a vector partial differential equation (PDE). General form of
% NSE is (in the inertial frame of reference):
%
% $$
% D(\rho \mbf{v})
% --------------- = -\nabla p + \nabla \cdot \mbf{T} + \mbf{f}          (1)
%       DT
% $$
%
% Note: The notation for math is adapted LaTeX syantax
%
% Note: In the Equation (1), rho is the fluid density, \mbf{v} is fluid
% velocity, p is pressure, \mbf{T} is second order tensor standing for
% deviatoric stress component of the total stress tensor and \mbf{f} is the
% external forcing applied on the fluid.
%
% The problem of finding a general solution for NSE is still an open one
% and very probably will not be ever solved. Furthermore, even prove of 
% existence and smoothness of such a solution was not yet published
% (and there is a substantial reward promised to the one who founds it)
%
% In engineering practice, the solution of Equation (1) for more
% complicated flows is approximated by various CFD methods. However, for
% the simplest of the cases, an approximate analytical solution may be
% found with the right set of simplifying assumptions.
%
% One type of such a flow is so called rivulet type flow on a smooth plate.
% In the case of neglecting the liquid spreading (assuming a constat width
% of the rivulet) Equation (1) can be substantially simplified and an
% analytical solution can be obtained (for details see ie Duffy and Moffat
% 1995).
%
% On the other hand, when the liquid spreading is not neglected, a simple
% analytical solution can no longer be obtained. However, in a very special
% case of,
% - wetted vertical plate
% - shallow narrow trickle of liquid
% - small but nonzero contact angles
% - ...
% at least a shape of gas-liquid interface of a spreading rivulet can be
% described by the equation
%
% $$
%
%
%          A                B
% x - ------------[\ln (-----------) - 1/5 ] + C = 0                    (2)
%     \beta^{15/4}      \beta^{3/4}
%  
%     4   A_             \eta        4   C_
% A = -- ------, B = ----------, C = -- ------                          (3)
%     15 \varpi      2 e^{2} l       15 \varpi
%
%              2\mu           27 \eta \mu 
% \vapri = ------------, A_ = -- ---------                              (4)
%          \rho g l^{2}       4   \gamma  
%
%
%      4        A_                       B
% C_ = -- ---------------- [\ln (---------------) - 1/5]                (5)
%      15 \beta_{0}^{15/4}       \beta_{0}^{3/4}
%
%
%           4 \mu Q
% \eta = (-----------)^{1/4}                                            (6)
%         105 \rho g
% $$
%
% Note: In equations (2)-(6), \beta stands for dynamic contact angle, x is
% the plate longitudinal coordinate, \rho is the liquid density, l is the
% thickness of a film layer wetting the plate, \mu is the liquid dynamic
% viscosity, \gamma is the liquid surface tension coefficient, g is the
% gravitational acceleration and Q is the liquid volumetric flow rate.
%
% Note: \beta_{0} is the dynamic contact angle at x = 0
%
% From the physical meaning of constants A, B and C, it can be stated that
%
% $$
% A > 0, B > 0, C > 0.
% $$
%
% Also, from the physical interpretation of the model above, we can neglect
% all the negative distances and put x non-negative (x >= 0).
%
% The gas-liquid interface profile shape itself is then for the combination
% ( x,\beta(x) ) defined by,
%
% $$
%            a \tan{\beta(x)}
% h(\zeta) = ---------------- (1 - \zeta^{2}), \zeta = y/a              (7)
%                 2
% $$
%
% Note: a is the rivulet half width and \zeta is the non-dimensional plate
% transversal coordinate.
%
% References: Will be given upon asking.

%% Questions:
% (i) Does the equation (2) implicitely define dependence
%
% $$
% \beta = \beta(x) \forall x \in \langle 0,\delta), \delta > 0
% $$
%
% (ii) Plot gas-liquid interface profiles at x = [0,0.05,0.1,0.2,0.4] for
% Q = 1e-7 m3/s
% \beta_{0} = 0.05
% \mu       = 1e-3 Pas
% \gamma    = 70 mN/m
% \rho      = 980 kg/m3
% l         = 3e-5 m
% g         = 9.81 m/s2
%                          1
% a         = \eta------------------
%                  \tan{\beta}^{3/4}

%% Solution - (ii)

function caller
% this function is more of a script. it gathers the problem parameters and
% calls the caller to solve it.

% clear the environment
close all
clc

% physical properties
l       = 3e-5;                                                             %wetting liquid layer thickness, m
g       = 9.81;                                                             %gravitational acceleration, m/s2
mu      = 1e-3;                                                             %liquid dynamic viscosity, Pas
gamma   = 10e-3;                                                            %liquid surface tension, N/m
rho     = 980;                                                              %liquid density, kg/m3
Q       = 1e-6;                                                             %liquid volumetric flow rate, m3/s

% initial condition
beta0   = 0.1;                                                             %initial dynamic contact angle, piRad

% define parameters used in solution
eta     = (105/4*mu*Q/(rho*g)).^(1/4);                                      %see Eq. (6)
A_      = 27/4.*eta*mu./gamma;                                              %see Eq. (4)
B       = eta./(2*exp(1)^2*l);                                              %see Eq. (3)
C_      = 4/15.*A_./beta0.^(15/4).*(log(B./beta0^(3/4)) - 1/5);             %see Eq. (5)

varpi   = 2*mu/(rho*g*l^2);                                                 %see Eq. (4)

A       = 4/15.*A_./varpi;                                                  %see Eq. (3)
C       = 4/15.*C_./varpi;

% auxiliary parameters and functions
aFun    = @(beta) eta./(tan(beta).^(3/4));                                  %get rivulet width from its conctact angle
hFun    = @(zeta,beta,a) a.*tan(beta)./2.*(ones(size(zeta)) - zeta.^2);     %function defining the rivulet g-l interface shape  

% create vector of parameters for the model
pars    = [A B C beta0];                                                   %beta0 is passed as initial guess to fsolve

% create vector of distances from the plate top
xVec     = [0.05 0.1 0.2 0.4];                                              %m

%%%%%%%%MODEL%%%%%%%%%%%%%%%
betaVec  = model(xVec,pars);                                                %call the model and get liquid layer thicknesses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% append initial condition
xVec    = [0 xVec]; betaVec = [beta0 betaVec];

% get rivulet half widths
aVec = aFun(betaVec);

% get rivulet g-l interface shape
zetaVec = 0:0.1:1;
hMat    = zeros(numel(xVec),numel(zetaVec));
parfor i = 1:numel(xVec)
    hMat(i,:) = hFun(zetaVec,betaVec(i),aVec(i));                           %can be rewritten w/o for loop, but would be more demanding 
end

figure('Units','Pixels','Position',[100 50 1000 700],'Visible','on');       %open a figure window
% graphical output - dependence of beta on x
subplot(221)
plot(xVec,betaVec,'ko','MarkerFaceColor','Black','MarkerSize',11);
title('\beta = \beta(x)','FontSize',18,'FontWeight','Bold')
xlabel('x, [m]','FontSize',18,'FontWeight','Normal');
ylabel('\beta, [\pi Rad]','FontSize',18,'FontWeight','Normal');

% graphical output - dependence of a on x
subplot(222)
plot(xVec,aVec,'k^','MarkerFaceColor','Black','MarkerSize',11);
title('a = a(x)','FontSize',18,'FontWeight','Bold')
xlabel('x, [m]','FontSize',18,'FontWeight','Normal');
ylabel('a, [m]','FontSize',18,'FontWeight','Normal');

% graphical output - g-l interface shape
subplot(2,2,3:4)
for i = 1:numel(xVec)
    plot(aVec(i).*zetaVec,hMat(i,:),'k','LineWidth',2)
    hold on    
end
title('h_{i} = h_{i}(\zeta), i = 1,2,...,numel(x)','FontSize',18,'FontWeight','Bold')
xlabel('x, [m]','FontSize',18,'FontWeight','Normal');
ylabel('h, [m]','FontSize',18,'FontWeight','Normal');


end


%% Model
function betaVec = model(xVec,pars)
% model that for given set of parameters, distance from the plate top and
% the mesh for interface shape function evaluation returns the shape
% dynamic contact angles along the rivulet

% get the parameters
A  = pars(1);B = pars(2);C = pars(3);beta0  = pars(4);                      %beta0 is an auxiliary parameter to get initial guess for fsolve

% allocate space for the solution
betaVec = zeros(1,numel(xVec));

% create lHS function
lhs = @(x,beta) ...                                                         %implicitly defined function for contact angle in dependence of 
    x - A./beta.^(15/4).*(log(B./beta.^(3/4)) - 1/5) + C;
% auxiliary function for the 

% solve the Eq. (2) to obtain \beta at given x
fopts   = optimset('Display','off','TypicalX',beta0);                       %set conditions for the numeric solver

parfor i = 1:numel(xVec)
    betaVec(i) = fsolve(@(beta) lhs(xVec(i),beta),beta0,fopts);             %find the current beta
end

end