Code covered by the BSD License

# A cavitation algorithm

25 Apr 2013 (Updated )

Matlab routines solving a linear complementarity problem appearing in lubrication with cavitation

%% A matlab routine solving a linear complementarity problem:
%
% A*u = f + B*eta, u'*eta=0, u>=0, 0<=eta<1
%
% appearing in lubrication with cavitation
%
% % Copyright (c) 2013 Andreas Almqvist, Andrew Spencer and Peter Wall
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1. Redistributions of source code must retain the above copyright notice,
%    this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

clear all;

%% Input parameters
%  Domain and boundary conditions
pc    = 50e3;  % Cavitation pressure [Pa]
pL1   = 100e3;  % Leading edge boundary pressure [Pa]
pL2   = 100e3;  % Trailing edge boundary pressure [Pa]
% Physical quantities
beta = 5e8;     % Bulk modulus [Pa]
V    = 1;       % Speed [m/s]
mu   = 0.01;    % Dynamic viscosity [Pas]
% Number of elements (Nodes = N+1)
N	 = 2^10;
% Defining the computational domain x\in(0,L) for the pocket slider
L  = 20e-3;     % = ap+bp+cp, see below.
dx = L/N;       % Grid size
x  = 0:L/N:L;   % Discrete grid
% Film thickness function representing a pocket slider bearing
ap = 2e-3;      % Inlet zone
bp = 3e-3;      % Pocket width
cp = L-ap-bp;   % Outlet zone
h1 = 10e-6;     % Pocket depth
h0 = 1e-6;      % Bearing clearance
h  = h1.*(ap<x & x<ap+bp) + h0;

%% Preprocessing
% Defining the coefficients and the BC's for the corresponding LCP problem
a = beta*(h.^3).'/12/mu;
b = -h.'*V/2;
F = h.'*V/2;
% Boundary coefficients
uL1	  = exp((pL1-pc)/beta)-1;	uL2 = exp((pL2-pc)/beta)-1;  % BC u(L1) and u(L2)
etaL1 = (uL1==0);	        etaL2	= (uL2==0);     % Complementarity
% Discretising the system
i = 2:N; % Indices for internal nodes
% Discretisation of d/dx(a*d/dx(u)+b*u) :=
% aw_i*u_(i-1)+am_i*u_i+ae_i*u_(i+1) where:
aw =  (a(i-1)+a(i))/(2*dx^2)-(b(i-1))/(2*dx);
am = -(a(i-1)+2*a(i)+a(i+1))/(2*dx^2);
ae =  (a(i)+a(i+1))/(2*dx^2)+(b(i+1))/(2*dx);
% Discretisation of d/dx(eta*F)
% by 1st order backward scheme, i.e. := bm_i*eta_i+bm_(i-1)*eta_(i-1)
bm = -F(i)/(dx);
bw =  F(i-1)/(dx);
be = -F(i+1)/(dx);
% Discretisation of dF/dx
z = (F(i+1)-F(i-1))/(2*dx);
% Create the matrix A and the matrix B
A = diag(aw(2:N-1),-1)+diag(am(1:N-1))+diag(ae(1:N-2),1);
B = diag(bw(2:N-1),-1)+diag(bm(1:N-1),0);
% Create rhs
f = z;
f(1)   = f(1)   - (aw(1)+bw(1))*uL1 + bw(1)*etaL1;
f(N-1) = f(N-1) - (ae(N-1)+be(N-1))*uL2;
% Create q and M in; u = q + M*eta <=> Au = f + B*eta
q = A^-1*f;     M = A^-1*B;

%% Solver
% Call the pivoting algorithm
tic; [u,eta,retcode] = LCPSolve(M,q,1e-15); toc;
u 	= [uL1;u;uL2];  eta = [etaL1;eta;etaL2];
% Compute the pressure and the saturation
p     = pc+beta*log(u+1);   % Converting solution to pressure
theta = 1-eta;

%% Postprocessing and Visualisation
% Derivative of u
ux  = [-3*u(1)+4*u(2)-u(3); u((2:N)+1)-u((2:N)-1);...
u(N-1)-4*u(N)+3*u(N+1)]/2/dx;
% The mass flow
q   = h.'.*u*V/2 + h.'*V/2-beta*(h.^3).'.*ux/(12*mu) - eta.*h.'*V/2;
legend('Bearing geometry (\mum)','Fluid pressure (bar)')