# Modeling and Simulation of Axisymmetric Swirling Flows

# Fluid Flow Simulations with Swirl Effects

Fluid flows with swirl effects can occur in rotationally symmetric geometries where the azimuthal or angular velocity component is non-zero. In this case one must solve the fully three-dimensional Navier-Stokes equations. However, due to the assumption that azimuthal variations can be neglected, it is still sufficient to limit the simulations to axisymmetric 2D domains, thus saving significant amounts of computational time and effort. This post will describe how to set up such a custom equation system in FEATool, and illustrate the modeling process with a selection of swirl flow CFD examples.

## Swirl Flow using the FEATool GUI

To model swirling flows with the FEATool GUI one can use the
axisymmetric Navier-Stokes physics mode. The *u* momentum equation in
the radial direction must be extended with the non-linear source term
*ρv*^{2}, which can be done using the *Edit
Equation* functionality in the *Subdomain Settings* dialog box.

Moreover, an additional equation for the radial velocity component must
also be added. This can be achieved by adding a *Custom Equation*
multiphysics mode with dependent variable named *v*, and corresponding
momentum equation expression defined as

```
r*rho_ns*v' - r*miu_ns*(vr_r + vz_z) + ...
miu_ns*v_r + r*rho_ns*(u*vr_t + w*vz_t) + rho_ns*u*v_t = ...
r*Ft + miu_ns*(v_r - 1/r*v_t)
```

Even though the dependent variables here are split between two different physics modes the solver will automatically solve them as a fully coupled system.

In this implementation approach it is necessary to manually set the
boundary condition for *v* separately, since swirl boundary conditions
are not predefined. Moreover, as the predefined axisymmetric
Navier-Stokes postprocessing expressions, such as velocity magnitude,
do not include the *v* component, one has to re-define these manually
as custom equation expressions if desired.

## Swirl Flow using the CLI

With FEATool one can also conveniently program and script the models on the Matlab and Octave 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.

## Taylor-Couette Flow Simulation

A very interesting phenomena happens when the rotational velocity is
increased, and after a certain point so called Taylor vortices
appear. The flow is still laminar but exhibits lengthwise vortices. As
the velocity is increased further they destabilize and start to travel
in the axial direction. The video below show this behavior, in
particular observe around *t = 10 s* when the stable Taylor vortices
appear (at about Taylor number *Ta = 1500 - 1700* which is in good
agreement with theory), and *t = 20 s* when they start moving length
wise.

The time-dependent Taylor-Couette flow script model can be downloaded
from the link below (note that periodic boundary conditions are
implemented with an external *solve_hook* function which will be
discussed in a future post).