FEATool Multiphysics
v1.16.1
Finite Element Analysis Toolbox

Fluid flows with swirl effects can occur in rotationally symmetric geometries with nonzero azimuthal or angular velocity. In this case one must generally solve the full 3D NavierStokes equations. However, if the azimuthal variations can be neglected, simulations can be limited to axisymmetric 2D domains and thus save significant amounts of computational time.
This example models axisymmetric flow between two concentric cylinders where the inner cylinder is rotating. This leads to swirl effects and so called TaylorCouette flow with periodic and inplane vortices.
This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Fluid Dynamics > TaylorCouette Swirling Flow from the File menu. Or alternatively, follow the stepbystep instructions below. Note that the CFDTool interface differ slightly from the FEATool Multiphysics instructions described in the following.
Click on the Axisymmetry Space Dimension selection button in the New Model dialog box, and select Swirl Flow from the Select Physics dropdown list. Leave the space dimension and dependent variable names to their default values. Finish the physics selection and close the dialog box by clicking on the OK button. (Note that for CFDTool the physics selection is done in the Equation settings dialog box.)
The geometry of consists of a rectangular cross section of the cylinder (axisymmetric geometries must be defined in the r>0 positive half plane).
1
into the x_{min} edit field.1.5
into the x_{max} edit field.1.5
into the y_{min} edit field.1.5
into the y_{max} edit field.Press OK to finish and close the dialog box.
The default grid may be too coarse to ensure an accurate solution. Decreasing the grid size and generating a finer grid can resolve the expected vortices better.
0.05
into the Subdomain Grid Size edit field.Press the Generate button to call the automatic grid generation algorithm.
Equation and material coefficients are specified in Equation/Subdomain mode. In the Equation Settings dialog box enter rho
for the density and miu
for the viscosity. Then press OK to finish with specifications the equation coefficients. (Note that the Equation Settings dialog box may look different for CFDTool.)
The Model Constants and Expressions functionality can be used to define and store convenient expressions which then are available in the point, equation, boundary coefficients, and as postprocessing expressions. Here it is used to define the fluid coefficients, angular velocity, inner and outer radius, and reference solution.
Name  Expression 

rho  1 
miu  1 
omega  50 
ri  1 
ro  1.5 
v_ref  omega*ri^2/(ro^2ri^2)*r+omega*ro^2/(ro^2ri^2)/r 
Note that FEATool can work with any unit system, and it is up to the user to use consistent units for geometry dimensions, material, equation, and boundary coefficients.
Boundary conditions are defined in Boundary Mode and describes how the model interacts with the external environment.
Then select the inner rotating boundary (number 4) in the left hand side Boundaries list box and select the Inlet/velocity boundary condition. Enter omega*ri
in the edit field to specify the velocity v_{0} in the tangentialdirection.
After the problem has been solved FEATool will automatically switch to postprocessing mode and display the computed velocity field.
To visualize the azimuthal velocity, open the Plot Options and postprocessing settings dialog box and select to plot and visualize the thvelocity as both surface and contour plots, and also activate the arrow plot.
One can study a section of the velocity profile by using the Point/Line Evaluation... feature from the Post menu. By entering a series of evaluation coordinates, both the evaluated expression and a corresponding cross section plot can be generated.
Create line plots of both the computed azimuthal velocity and the analytic expression previously defined as v_ref
.
1:0.01:1.5
into the Evaluation coordinates in rdirection edit field.0
into the Evaluation coordinates in zdirection edit field.Press the Apply button.
v_ref
into the edit field.Press OK to finish and close the dialog box.
From comparing the curves it should be clear that the simulation produces the expected results.
If the rotational velocity is increased beyond the critical Taylor number periodic inplane vortices will appear.
Open the Model Constants and Expressions... dialog box again and increase omega to 175
.
1e4
into the Nonlinear stopping criteria for solution defects edit field.275
into the Maximum number of nonlinear iterations edit field.Once the solver has finished and converged, plot and visualize the inplane velocity sqrt(u^2+w^2)
to clearly see the Taylor vortices.
The TaylorCouette swirling flow fluid dynamics model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB mscript text file, or GUI script (.fes) file.
With FEATool one can also conveniently program and script the models on the MATLAB command line interfaces (CLI). The following example sets up a model for flow between two cylinders where the inner one is rotating and the outer is stationary. The first step is to define the model parameters, geometry, and grid
% Problem definition. rho = 1.0; % Density. miu = 1.0; % Molecular viscosity. omega = 5.0; % Angular velocity. ri = 0.5; % Inner radius. ro = 1.5; % Outer radius. h = 3.0; % Height of cylinder. % Geometry and grid generation. fea.sdim = {'r' 'z'}; fea.geom.objects = { gobj_rectangle(ri,ro,h/2,h/2) }; fea.grid = gridgen( fea, 'hmax', (rori)/6 );
Next the, governing equations are defined, here directly without using any predefined physics modes
% Equation definition. fea.dvar = { 'u', 'v', 'w', 'p' }; fea.sfun = [ repmat( {'sflag2'}, 1, 3 ) {'sflag1'} ]; 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 =  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 = 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 = 0'; 'r*ur_t + r*wz_t + u_t = 0' }; fea.eqn = parseeqn( c_eqn, fea.dvar, fea.sdim ); fea.coef = { 'rho', rho ; 'miu', miu };
For boundary conditions the inner wall is given a rotational velocity, while the outer is prescribed a noslip wall, and the top and bottom zero normal flows
% Boundary conditions. fea.bdr.d = { [] 0 [] 0 ; [] 0 [] omega*ri ; 0 0 0 [] ; [] [] [] [] }; fea.bdr.n = cell(size(fea.bdr.d));
The incompressible NavierStokes equations typically requires the pressure to be defined and fixed in at least one point. This is usually done with an outflow boundary, but since here there are no in or outflows, the pressure is simply fixed by setting the value of the top outer corner to zero
% Fix pressure at p([r,z]=[ro,h/2]) = 0. [~,ix_p] = min( sqrt( (fea.grid.p(1,:)ro).^2 + (fea.grid.p(2,:)h/2).^2) ); fea.pnt = struct( 'type', 'constr', 'index', ix_p, 'dvar', 'p', 'expr', '0' );
To improve convergence, and compute a stationary solution faster, an analytic Newton Jacobian is defined (the nonlinear Newton solver is activated when the jac solver argument is provided, otherwise standard fixpoint iterations are used)
% Define analytical Newton Jacobian bilinear form. jac.coef = { 'r*rho*ur' 'rho*v' 'r*rho*uz' [] ; 'r*rho*vr+rho*v' [] 'r*rho*vz' [] ; 'r*rho*wr' [] 'r*rho*wz' [] ; [] [] [] [] }; jac.form = cell(size(jac.coef)); [jac.form{~cellfun(@isempty,jac.coef)}] = deal([1;1]);
The complete problem can now be parsed and solved
% Parse and solve problem. fea = parseprob( fea ); fea.sol.u = solvestat( fea, 'jac', jac );
To verify the accuracy of the computed solution a comparison with an analytic solution (valid for small rotational velocities) is made
% Exact (analytical) solution. a =  omega*ri^2 / (ro^2ri^2); b = omega*ri^2*ro^2 / (ro^2ri^2); v_th_ex = @(r,a,b) a.*r + b./r; % Postprocessing. subplot(1,2,1) postplot( fea, 'surfexpr', 'sqrt(u^2+v^2+w^2)', 'isoexpr', 'v' ) subplot(1,2,2) hold on grid on r = linspace( ri, ro, 100 ); v_th = evalexpr( 'v', [r;zeros(1,length(r))], fea ); plot( r, v_th, 'b' ) r = linspace( ri, ro, 10 ); plot( r, v_th_ex(r,a,b), 'r.' ) legend( 'Computed solution', 'Exact solution') xlabel( 'Radius, r') ylabel( 'Angular velocity, v')
From the figure below it is clear that very good agreement with the analytical solution has been achieved.
A timedependent TaylorCouette flow script model is available from the FEATool models and examples directory as the ex_swirl_flow3 model mscript file (note that periodic boundary conditions are implemented with an external solve_hook function.