 |
FEATool Multiphysics
v1.16.5
Finite Element Analysis Toolbox
|
|
FEATool Multiphysics
v1.16.5
Finite Element Analysis Toolbox
|
EX_NAVIERSTOKES3 2D Example for incompressible stationary flow around a cylinder.
[1] John V, Matthies G. Higher-order finite element discretizations in a benchmark problem for incompressible flows. International Journal for Numerical Methods in Fluids 2001.
[2] Nabh G. On higher order methods for the stationary incompressible Navier-Stokes equations. PhD Thesis, Universitaet Heidelberg, 1998.
Accepts the following property/value pairs.
Input Value/{Default} Description
-----------------------------------------------------------------------------------
rho scalar {1} Density
miu scalar {0.001} Molecular/dynamic viscosity
umax scalar {0.3} Maximum magnitude of inlet velocity
igrid scalar {2} Grid type: >0 regular (igrid refinements)
<0 unstruc. grid (with hmax=|igrid|)
sf_u string {sflag1} Shape function for velocity
sf_p string {sflag1} Shape function for pressure
iphys scalar 0/{1} Use physics mode to define problem (=1)
solver string {} Solver selection default, openfoam, su2, fenics
iplot scalar 0/{1} Plot solution (=1)
.
Output Value/(Size) Description
-----------------------------------------------------------------------------------
fea struct Problem definition struct
out struct Output struct
cOptDef = { ...
'rho', 1;
'miu', 0.001;
'umax', 0.3;
'igrid', 2;
'sf_u', 'sflag1';
'sf_p', 'sflag1';
'iphys', 1;
'solver', '';
'iplot', 1;
'tol', [0.05 0.35 0.1];
'fid', 1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid = opt.fid;
% Model parameters.
rho = opt.rho; % Density.
miu = opt.miu; % Molecular/dynamic viscosity.
umax = opt.umax; % Maximum magnitude of inlet velocity.
umean = 2/3*umax; % Mean inlet velocity.
% Geometry and grid parameters.
h = 0.41; % Height of rectangular domain.
l = 2.2; % Length of rectangular domain.
xc = 0.2; % x-coordinate of cylinder center.
yc = 0.2; % y-coordinate of cylinder center.
diam = 0.1; % Diameter of cylinder.
% Discretization parameters.
sf_u = opt.sf_u; % FEM shape function type for velocity.
sf_p = opt.sf_p; % FEM shape function type for pressure.
% Grid generation.
fea.sdim = { 'x' 'y' };
if( opt.igrid>=1 )
fea.grid = cylbenchgrid( opt.igrid );
fea.grid.s(:) = 1;
else
gobj1 = gobj_rectangle( 0, 2.2, 0, 0.41, 'R1' );
gobj2 = gobj_circle( [0.2 0.2], 0.05, 'C1' );
fea.geom.objects = { gobj1 gobj2 };
fea = geom_apply_formula( fea, 'R1-C1' );
fea.grid = gridgen( fea, 'hmax', abs(opt.igrid), 'fid', fid );
end
if( opt.iphys==1 && strcmp(opt.solver,'fenics') && size(fea.grid.c,1)==4 )
warning( 'Converting quadrilateral mesh to triangular to support FEniCS.' )
fea.grid = quad2tri( fea.grid );
end
n_bdr = max(fea.grid.b(3,:)); % Number of boundaries.
% Boundary conditions.
dtol = sqrt(eps)*1e3;
i_inflow = findbdr( fea, ['x<=',num2str(dtol)] ); % Inflow boundary number.
i_outflow = findbdr( fea, ['x>=',num2str(l-dtol)] ); % Outflow boundary number.
s_inflow = ['4*',num2str(umax),'*(y*(',num2str(h),'-y))/',num2str(h),'^2']; % Definition of inflow profile.
i_cyl = findbdr( fea, ['sqrt((x-',num2str(xc),').^2+(y-',num2str(yc),').^2)<=(',num2str(diam/2+dtol),')'] ); % Cylinder boundary number.
% Problem definition.
if ( opt.iphys==1 )
fea = addphys(fea,@navierstokes); % Add Navier-Stokes equations physics mode.
fea.phys.ns.eqn.coef{1,end} = { rho };
fea.phys.ns.eqn.coef{2,end} = { miu };
fea.phys.ns.sfun = { sf_u sf_u sf_p }; % Set shape functions.
fea.phys.ns.bdr.sel(i_inflow) = 2;
fea.phys.ns.bdr.sel(i_outflow) = 4;
fea.phys.ns.bdr.coef{2,end}{1,i_inflow} = s_inflow; % Set inflow profile.
fea = parsephys(fea); % Check and parse physics modes.
else
fea.dvar = { 'u' 'v' 'p' }; % Dependent variable name.
fea.sfun = { sf_u sf_u sf_p }; % Shape function.
% Define equation system.
cvelx = [num2str(rho),'*',fea.dvar{1}]; % Convection velocity in x-direction.
cvely = [num2str(rho),'*',fea.dvar{2}]; % Convection velocity in y-direction.
fea.eqn.a.form = { [2 3 2 3;2 3 1 1] [2;3] [1;2];
[3;2] [2 3 2 3;2 3 1 1] [1;3];
[2;1] [3;1] [] };
fea.eqn.a.coef = { {2*miu miu cvelx cvely} miu -1;
miu {miu 2*miu cvelx cvely} -1;
1 1 [] };
fea.eqn.f.form = { 1 1 1 };
fea.eqn.f.coef = { 0 0 0 };
% Define boundary conditions.
fea.bdr.d = cell(3,n_bdr);
[fea.bdr.d{1:2,:}] = deal( 0);
fea.bdr.d{1,i_inflow} = s_inflow;
[fea.bdr.d{:,i_outflow }] = deal([]);
fea.bdr.n = cell(3,n_bdr);
end
% Parse and solve problem.
fea = parseprob(fea); % Check and parse problem struct.
if( opt.iphys==1 && strcmp(opt.solver,'openfoam') )
logfid = fid; if( ~got.fid ), fid = []; end
fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid );
fid = logfid;
elseif( opt.iphys==1 && strcmp(opt.solver,'su2') )
logfid = fid; if( ~got.fid ), fid = []; end
fea.sol.u = su2( fea, 'fid', fid, 'logfid', logfid );
fid = logfid;
elseif( opt.iphys==1 && strcmp(opt.solver,'fenics') )
fea = fenics( fea, 'fid', fid );
else
jac.form = {[1;1] [1;1] [];[1;1] [1;1] []; [] [] []};
jac.coef = {[num2str(rho),'*ux'] [num2str(rho),'*uy'] []; [num2str(rho),'*vx'] [num2str(rho),'*vy'] []; [] [] []};
fea.sol.u = solvestat( fea, 'fid', fid, 'nsolve', 2, 'jac', jac, 'nlrlx', '(1+(it>2))/2' ); % Call to stationary solver.
end
% Postprocessing.
s_velm = 'sqrt(u^2+v^2)';
if ( opt.iplot>0 )
figure
subplot(2,1,1)
postplot( fea, 'surfexpr', s_velm )
title( 'Velocity field' )
subplot(2,1,2)
postplot( fea, 'surfexpr', 'p' )
title( 'Pressure' )
end
% Calculate benchmark quantities (line integration method).
s_tfx = ['nx*p+',num2str(miu),'*(-2*nx*ux-ny*(uy+vx))'];
s_tfy = ['ny*p+',num2str(miu),'*(-nx*(vx+uy)-2*ny*vy)'];
s_cd = ['2*(',s_tfx,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
s_cl = ['2*(',s_tfy,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
i_cub = 10;
c_d1 = intbdr(s_cd,fea,i_cyl,i_cub);
c_l1 = intbdr(s_cl,fea,i_cyl,i_cub);
dp = evalexpr('p',[0.15 0.25;0.2 0.2],fea);
% Calculate benchmark quantities (volume integration method).
bdrm = fea.bdr.bdrm{1};
ind_b = [];
ind_bm = [];
for ii=i_cyl
ind_b = [ind_b find(fea.grid.b(3,:)==ii)];
ind_bm = [ind_bm find(bdrm(3,:)==ii)];
end
ind_c = fea.grid.b(1,ind_b);
ind_gdof = bdrm(4,ind_bm);
% Create field 'a' with values one on the cylinder and zero everywhere else.
fea.dvar = [ fea.dvar, {'a'} ];
fea.sfun = [ fea.sfun, fea.sfun(1) ];
fea = parseprob(fea);
n_dof = max(fea.eqn.dofm{1}(:));
u_a = zeros(n_dof,1);
u_a(ind_gdof) = 1;
fea.sol.u= [fea.sol.u;u_a];
fea.eqn = struct;
fea.bdr = struct;
fea = parseprob(fea);
s_tfx = ['ax*p+',num2str(miu),'*(-2*ax*ux-ay*(uy+vx))-(u*ux+v*uy)*a'];
s_tfy = ['ay*p+',num2str(miu),'*(-ax*(vx+uy)-2*ay*vy)-(u*vx+v*vy)*a'];
s_cd = ['2*(',s_tfx,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
s_cl = ['2*(',s_tfy,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
c_d2 = intsubd(s_cd,fea,[],[],3);
c_l2 = intsubd(s_cl,fea,[],[],3);
if( ~isempty(fid) )
fprintf(fid,'\n\nBenchmark quantities:\n\n')
fprintf(fid,'Drag coefficient, cd = %6f (l), %6f (v) (Ref: 5.579535)\n',c_d1,c_d2)
fprintf(fid,'Lift coefficient, cl = %6f (l), %6f (v) (Ref: 0.010619)\n',c_l1,c_l2)
fprintf(fid,'Pressure, dp = %6f (Ref: 0.117520)\n',dp(1)-dp(2))
end
% Error checking.
out.cd = [c_d1 c_d2];
out.cl = [c_l1 c_l2];
out.dp = dp(1)-dp(2);
out.err = [abs(out.cd-5.579535)/5.579535;
abs(out.cl-0.010619)/0.010619;
abs(dp(1)-dp(2)-0.117520)/0.117520 0];
out.pass = (out.err(1,2)<opt.tol(1))&&(out.err(2,2)<opt.tol(2))&&(out.err(3,1)<opt.tol(3));
if ( nargout==0 )
clear fea out
end
%