FEATool Multiphysics  v1.16.6 Finite Element Analysis Toolbox
ex_multiphase3.m File Reference

## Description

EX_MULTIPHASE3 Breaking dam multiphase flow example.

[ FEA, OUT ] = EX_MULTIPHASE3( VARARGIN ) Breaking dam multiphase flow benchmark problem as defined in Martin J, Moyce W. Part IV. An Experimental study of the collapse of liquid columns on a rigid horizontal plane. Phil. Trans. R. Soc. Lond. A 1952; 244(882):312–324, 1952, DOI: 10.1098/rsta.1952.0006. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
rho1        scalar {998}           Density of the water
rho2        scalar {1.205}         Density of the surrounding air
miu1        scalar {0.001}         Viscosity of the water
miu2        scalar {1.983e-5}      Viscosity of the surrounding air
gy          scalar {9.82}          Gravitational constant
a           scalar {0.05715}       Width and height of the water column
igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
hmax        scalar {a/12}          Max grid cell size
dt          scalar {}              Time step size
sf_u        string {sflag2}        Shape function for velocity
sf_p        string {sflag1}        Shape function for pressure
sf_c        string {sflag1}        Shape function for the level set function
iplot       scalar 0/{1}           Plot solution and error (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { ...
'rho1',     998; ...
'rho2',     1.205; ...
'miu1',     0.001; ...
'miu2',     1.983e-5; ...
'gy',       9.82; ...
'a',        0.05715; ...
'igrid',    0; ...
'hmax',     0.05715/18; ...
'dt',       {}; ...
'sf_u',     'sflag2';  ...
'sf_p',     'sflag1'; ...
'sf_c',     'sflag1'; ...
'iplot',    1; ...
'itest',    0; ...
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;
if( ~got.dt )
TSTEP = 0.05/sqrt(opt.gy/opt.a);
else
TSTEP = opt.dt;
end
TMAX= 3.0/sqrt(opt.gy/opt.a);
MAXNIT = 30;
if( opt.itest )
TMAX = TSTEP;
MAXNIT = 1;
end

% Geometry and grid generation.
l = 6*opt.a;
h = 1.5*opt.a;
fea.grid = rectgrid(round(l/opt.hmax),round(h/opt.hmax),[0 l;0 h]);
if( opt.igrid==1 )
end
fea.sdim = { 'x' 'y' };   % Coordinate names.
n_bdr = max(fea.grid.b(3,:)) + 1;   % Increment number of boundaries.
fea.grid.b(3,1) = n_bdr;

% Problem definition.
dh    = num2str(1.5*opt.hmax);      % Smoothing region width.
dw    = ['(c/',dh,'/(sqrt(cx^2+cy^2+eps)))'];
smhs  = ['((0.5*(1+',dw,'+1/pi*sin(pi*',dw,')))*(',dw,'>-1)*(',dw,'<1)+(',dw,'>=1))'];   % Smooth heaviside function.
rho   = [num2str(opt.rho1),'+',num2str(opt.rho2-opt.rho1),'*',smhs];
miu   = [num2str(opt.miu1),'+',num2str(opt.miu2-opt.miu1),'*',smhs,'+2e2*',num2str(opt.hmax),'*sqrt(u^2+v^2+1e-2)'];

% Add Navier-Stokes equations physics mode.
fea.phys.ns.eqn.coef{1,end} = { rho };
fea.phys.ns.eqn.coef{2,end} = { miu };
fea.phys.ns.eqn.coef{4,end} = { ['-(',rho,')*',num2str(opt.gy)] };
fea.phys.ns.bdr.sel(n_bdr)  = 4;   % Set pressure to zero on last boundary segment.
fea.phys.ns.sfun            = { opt.sf_u opt.sf_u opt.sf_p };

fea.phys.cd.sfun            = { opt.sf_c };
fea.phys.cd.eqn.coef{2,4}   = { 0.001 };   % Set diffusion coefficient.
fea.phys.cd.eqn.coef{3,4}   = { 'u' };    % Convection velocity in x-direction.
fea.phys.cd.eqn.coef{4,4}   = { 'v' };    % Convection velocity in y-direction.

% Parse physics modes.
fea = parsephys(fea);

% Implement slip boundary conditions on left and bottom walls.
fea.bdr.d{1}{1} = [];
fea.bdr.d{1}{5} = [];
fea.bdr.d{2}{4} = [];
fea.bdr.n{1}{1} = 0;
fea.bdr.n{1}{5} = 0;
fea.bdr.n{2}{5} = 0;

% Parse problem.
fea = parseprob(fea);

% Call to time-dependent solver.
init = { '0', '0', '0', ['(x-',num2str(opt.a),')*((x-',num2str(opt.a),')>(y-',num2str(opt.a),'))+(y-',num2str(opt.a),')*((x-',num2str(opt.a),')<=(y-',num2str(opt.a),'))'] };
fea.sol.u = solvetime( fea, ...
'fid',     fid, ...
'tmax',    TMAX, ...
'tstep',   TSTEP, ...
'icub',    3, ...
'maxnit',  MAXNIT, ...
'init',    init, ...
'ischeme', 3 );

% Postprocessing.
if ( opt.iplot>0 )
figure('position',[0 0 1920 1080])
subplot(2,2,1)
postplot(fea,'surfexpr','sqrt(u^2+v^2)')
title('Velocity field')
subplot(2,2,2)
postplot(fea,'surfexpr','p','evalstyle','exaxct')
title('Pressure')
subplot(2,2,3)
postplot(fea,'surfexpr','c')
title('Level set field')
subplot(2,2,4)
postplot(fea,'surfexpr',rho)
% postplot(fea,'isoexpr','c','isolev',[0 0],'arrowexpr',{'u','v'})
plot([0 l l 0 0],[0 0 h h 0],'k')
title('Interface')
end

out.err  = [];
out.pass = [];
if ( nargout==0 )
clear fea out
end