FEATool Multiphysics
v1.11 Finite Element Analysis Toolbox |

Taylor-Couette Swirling Flow

Fluid flows with swirl effects can occur in rotationally symmetric geometries with non-zero azimuthal or angular velocity. In this case one must generally solve the full 3D Navier-Stokes 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 Taylor-Couette flow with periodic and in-plane vortices.

This model is available as an automated tutorial by selecting **Model Examples and Tutorials...** > **Fluid Dynamics** > **Taylor-Couette Swirling Flow** from the **File** menu. Or alternatively, follow the step-by-step instructions below.

- To start a new model click the
**New Model**toolbar button, or select*New Model...*from the*File*menu. Click on the

**Axi 2D***Space Dimension*selection button in the*New Model*dialog box, and select**Swirl Flow**from the*Select Physics*drop-down 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.

The geometry of consists of a rectangular cross section of the cylinder (axisymmetric geometries must be defined in the *r>0* positive half plane).

- Select
**Rectangle**from the*Geometry*menu. - Enter
`1`

into the*x*edit field._{min} - Enter
`1.5`

into the*x*edit field._{max} - Enter
`-1.5`

into the*y*edit field._{min} - Enter
`1.5`

into the*y*edit field._{max} Press

**OK**to finish and close the dialog box.- Switch to
**Grid**mode by clicking on the corresponding*Mode Toolbar*button.

The default grid may be too coarse ensure an accurate solution. Decreasing the grid size and generating a finer grid can resolve the expected vortices better.

- Enter
`0.05`

into the*Grid Size*edit field. Press the

**Generate**button to call the automatic grid generation algorithm.- Switch to
**Equation**mode by clicking on the corresponding*Mode Toolbar*button. 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.

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.

- Press the
**Constants***Toolbar*button, or select the corresponding entry from the*Equation*menu, and enter the following variables in the*Model Constants and Expressions*dialog box. Press*Enter*after the last expression or use the**Add Row**button to expand the expression list.

Name | Expression |
---|---|

rho | 1 |

miu | 1 |

omega | 50 |

ri | 1 |

ro | 1.5 |

v_ref | -omega*ri^2/(ro^2-ri^2)*r+omega*ro^2/(ro^2-ri^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.

- Switch to boundary condition specification mode by clicking on
**Boundary**the*Mode Toolbar*button.

Boundary conditions are defined in *Boundary Mode* and describes how the model interacts with the external environment.

- In the
*Boundary Settings*dialog box, first choose boundaries**1-3**in the left hand side*Boundaries*list box and select the**Wall/no-slip**boundary condition from the drop-down menu. 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*in the tangential-direction._{0}- Press
**OK**to finish the boundary condition specification. - Now that the problem is fully specified, press the
**Solve***Mode Toolbar*button to switch to solve mode. Then press the**=***Tool*button to call the solver with the default solver settings.

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**th-velocity**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`

.

- Select
**Point/Line Evaluation...**from the*Post*menu. - Enter
`1:0.01:1.5`

into the*Evaluation coordinates in r-direction*edit field. - Enter
`0`

into the*Evaluation coordinates in z-direction*edit field. - Select
**th-velocity**from the*Evaluation Expression*drop-down menu. Press the

**Apply**button.- Enter
`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 in-plane vortices will appear.

Open the the **Model Constants and Expressions...** dialog box again and increase **omega** to `175`

.

- Return to
**Solve**mode again, open the**Solver Settings**dialog box. - Enter
`1e-4`

into the*Non-linear stopping criteria for solution defects*edit field. - Enter
`275`

into the*Maximum number of non-linear iterations*edit field. - Press
**OK**to apply the selected solver settings. - This time use the
**Restart**button to start the solver with the previous solution as initial guess. This can help with convergence for highly non-linear problems. Alternatively, using a time-dependent solver is another option. Once the solver has finished and converged, plot and visualize the in-plane velocity

`sqrt(u^2+w^2)`

to clearly see the Taylor vortices.

The *Taylor-Couette 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 m-script 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', (ro-ri)/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 no-slip 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 Navier-Stokes 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 out-flows, 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 non-linear Newton solver is activated when the *jac* solver argument is provided, otherwise standard fix-point 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^2-ri^2); b = omega*ri^2*ro^2 / (ro^2-ri^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 time-dependent Taylor-Couette flow script model is available from the FEATool models and examples directory as the ex_swirl_flow3 model m-script file (note that periodic boundary conditions are implemented with an external *solve_hook* function.