Modeling and Simulation of Axisymmetric Swirl Flow

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 ρv2, 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).

FEATool Taylor-Couette Swirl Flow Model

Originally published at www.featool.com on September 20, 2017.