FEATool  v1.7.1
Finite Element Analysis Toolbox
 All Files Functions Pages
ex_swirl_flow3.m File Reference

Description

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.

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 {sflag2}        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

Code listing

 cOptDef = { 'rho',      1;
             'miu',      1;
             'omega',    300;
             'tmax',     3;
             'nstep',    300;
             'ri',       1.0
             'ro',       1.5;
             'h',        3;
             'sf_u',     'sflag2';
             'sf_p',     'sflag1';
             '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 -h/2 ro h/2] );


% Equation definition.
 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*Ft + 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 ;
              'Ft',  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));


% Save pediodic boundary function to file.
 s_file = fileread( [mfilename('fullpath'),'.m'] );
 ix = strfind( s_file, 'function' );
 s_fcn = s_file(ix(end):end);
 ix1 = find( s_fcn == '=', 1 );
 ix2 = find( s_fcn == '(', 1 );
 s_fname = strtrim(s_fcn(ix1+1:ix2-1));
 fid_tmp = fopen( [s_fname,'.m'], 'w' );
 fprintf( fid_tmp, '%s', s_fcn );
 fclose( fid_tmp );

 [fea.bdr.n{1:3,1}] = deal( ['solve_hook_',s_fname] );   % Set periodic BCs.


% 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 );
 fea.sol.u = solvetime( fea, 'ischeme', 1, 'tstep', opt.tmax/opt.nstep, 'tmax', opt.tmax, 'tolchg', inf, 'fid', fid );
 delete( [s_fname,'.m'] )


% 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


%