|
FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
EX_SWIRL_FLOW3 2D Axisymmetric Taylor-Couette (swirl) flow.
[ FEA, OUT ] = EX_SWIRL_FLOW3( VARARGIN ) Axisymmetric Taylor-Couette swirl flow in a tubular region where the inner cylindrical wall is rotating. Time dependent solution with periodic top and bottom boundaries.
Accepts the following property/value pairs.
Input Value/{Default} Description
-----------------------------------------------------------------------------------
rho scalar {1} Density
miu scalar {1} Molecular/dynamic viscosity
omega scalar {300} Maximum angular velocity (of inner wall)
tmax scalar {3} Maximum time
nstep scalar {300} Number of time steps
ri scalar {1.0} Inner radius
ro scalar {1.5} Outer radius
h scalar {3} Height of cylinder
sf_u string {sflag1} Shape fcn for velocity
sf_p string {sflag1} Shape fcn for pressure
iplot scalar 0/{1} Plot solution (=1)
.
Output Value/(Size) Description
-----------------------------------------------------------------------------------
fea struct Problem definition struct
out struct Output struct
cOptDef = { 'rho', 1;
'miu', 1;
'omega', 300;
'tmax', 3;
'nstep', 300;
'ri', 1.0
'ro', 1.5;
'h', 3;
'sf_u', 'sflag1';
'sf_p', 'sflag1';
'iphys', 1;
'iplot', 1;
'fid', 1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid = opt.fid;
% Geometry and grid generation.
fea.sdim = {'r' 'z'};
ri = opt.ri; % Inner radius.
ro = opt.ro; % Outer radius.
h = opt.h; % Height of cylinder.
fea.geom.objects = { gobj_rectangle(ri,ro,-h/2,h/2) };
n = 16;
fea.grid = rectgrid( n, 6*n, [ri ro; -h/2 h/2] );
% Equation definition.
if ( opt.iphys==1 )
fea = addphys(fea,@swirlflow);
fea.phys.sw.eqn.coef{1,end} = { opt.rho };
fea.phys.sw.eqn.coef{2,end} = { opt.miu };
fea.phys.sw.sfun = [ repmat( {opt.sf_u}, 1, 3 ) {opt.sf_p} ];
fea.phys.sw.bdr.sel = [3 1 3 2];
fea.phys.sw.bdr.coef{2,end}{2,4} = sprintf( '%g*t/%g', opt.omega*ri, opt.tmax );
fea = parsephys(fea);
else
opt.sf_u = 'sflag2';
opt.sf_p = 'sflag1';
fea.dvar = { 'u', 'v', 'w', 'p' };
fea.sfun = [ repmat( {opt.sf_u}, 1, 3 ) {opt.sf_p} ];
c_eqn = { 'r*rho*u'' - r*miu*(2*ur_r + uz_z + wr_z) + r*rho*(u*ur_t + w*uz_t) + r*p_r = r*Fr - 2*miu/r*u_t + p_t + rho*v*v_t';
'r*rho*v'' - r*miu*( vr_r + vz_z) + miu*v_r + r*rho*(u*vr_t + w*vz_t) + rho*u*v_t = r*Fth + miu*(v_r - 1/r*v_t)';
'r*rho*w'' - r*miu*( wr_r + uz_r + 2*wz_z) + r*rho*(u*wr_t + w*wz_t) + r*p_z = r*Fz';
'r*ur_t + r*wz_t + u_t = 0' };
fea.eqn = parseeqn( c_eqn, fea.dvar, fea.sdim );
fea.coef = { 'rho', opt.rho ;
'miu', opt.miu ;
'Fr', 0 ;
'Fth', 0 ;
'Fz', 0 };
% Boundary conditions.
fea.bdr.d = { [] 0 [] 0 ;
[] 0 [] sprintf( '%g*t/%g', opt.omega*ri, opt.tmax ) ;
[] 0 [] 0 ;
[] [] [] [] };
fea.bdr.n = cell(size(fea.bdr.d));
end
% Fix pressure at p([r,z]=[ro,0]) = 0.
[~,ix_p] = min( sqrt( (fea.grid.p(1,:)-ro).^2 + (fea.grid.p(2,:)-0).^2) );
fea.pnt = struct( 'type', 'constr', ...
'index', ix_p, ...
'dvar', 'p', ...
'expr', '0' );
% Parse and solve problem.
fea = parseprob( fea );
if ( opt.iphys==1 )
fea.bdr.n{1}{1} = @periodic_vel_bc_1_3;
fea.bdr.n{2}{1} = @periodic_vel_bc_1_3;
fea.bdr.n{3}{1} = @periodic_vel_bc_1_3;
else
[fea.bdr.n{1:3,1}] = deal( @periodic_vel_bc_1_3 ); % Set periodic BCs.
end
fea.sol.u = solvetime( fea, 'ischeme', 1, 'tstep', opt.tmax/opt.nstep, 'tmax', opt.tmax, 'tolchg', inf, 'fid', fid );
% Postprocessing.
if( opt.iplot )
subplot(1,2,1)
postplot( fea, 'surfexpr', 'sqrt(u^2+w^2)', 'isoexpr', 'sqrt(u^2+w^2)', ...
'arrowexpr', {'u' 'w'}, 'arrowcolor', 'w', 'arrowspacing', [8 48] )
title('In-plane velocity')
subplot(1,2,2)
postplot( fea, 'surfexpr', 'v', 'isoexpr', 'v' )
title('Azimuthal velocity')
end
out = [];
if( nargout==0 )
clear fea out
end
%