|
FEATool Multiphysics
v1.17.5
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.5
Finite Element Analysis Toolbox
|
EX_NAVIERSTOKES6 2D incompressible time-dependent flow around a cylinder.
[ FEA, OUT ] = EX_NAVIERSTOKES6( VARARGIN ) Benchmark example for time-dependent flow around a cylinder.
[1] John V. Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder. International Journal for Numerical Methods in Fluids 2004; 44:777-788.
Accepts the following property/value pairs.
Input Value/{Default} Description
-----------------------------------------------------------------------------------
rho scalar {1} Density
miu scalar {0.001} Molecular/dynamic viscosity
igrid scalar {2} Grid type: >0 regular (igrid refinements)
<0 unstruc. grid (with hmax=|igrid|)
dt scalar {0.02} Time step size
instatbc boolean {true} Use instationary boundary conditions
ischeme scalar {2} Time stepping scheme
solver string {} Solver selection default, openfoam, su2
sf_u string {sflag2} Shape function for velocity
sf_p string {sf_disc1} Shape function for pressure
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;
'igrid', 2;
'dt', 0.02;
'instatbc' true;
'ischeme' 2;
'solver', '';
'sf_u', 'sflag2';
'sf_p', 'sf_disc1';
'iplot', 1;
'tol', [0.01 0.01 0.1 0.5 0.2];
'fid', 1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid = opt.fid;
if( any(strcmp(opt.solver,{'su2'})) )
if( ~isempty(opt.fid) )
warning('SU2 does not support intationary BCs.')
end
opt.instatbc = false;
end
if( strcmp(opt.solver,'openfoam') && ~got.igrid )
opt.igrid = 3;
end
% 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
% Boundary conditions.
inflow_boundary = findbdr( fea, ['x<=',num2str(sqrt(eps))] ); % Inflow boundary number.
outflow_boundary = findbdr( fea, ['x>=',num2str(2.1)] ); % Outflow boundary number.
if( opt.instatbc )
inflow_bc = '6*sin(pi*t/8)*(y*(0.41-y))/0.41^2';
else
inflow_bc = '6*(y*(0.41-y))/0.41^2';
end
% Problem definition.
fea = addphys( fea, @navierstokes );
fea.phys.ns.eqn.coef{1,end} = { opt.rho };
fea.phys.ns.eqn.coef{2,end} = { opt.miu };
if( any(strcmp(opt.solver,{'openfoam','su2'})) )
fea.phys.ns.sfun = { 'sflag1', 'sflag1', 'sflag1' };
else
fea.phys.ns.sfun = { opt.sf_u, opt.sf_u, opt.sf_p };
end
fea.phys.ns.bdr.sel( inflow_boundary ) = 2;
fea.phys.ns.bdr.sel( outflow_boundary ) = 4;
fea.phys.ns.bdr.coef{2,end}{1,inflow_boundary} = inflow_bc;
fea = parsephys( fea );
% Call solver
if( strcmp(opt.solver,'openfoam') )
if( opt.ischeme==1 )
ddtSchemes = 'backward';
else
ddtSchemes = 'CrankNicolson 0.9';
end
logfid = fid; if( ~got.fid ), fid = []; end
[fea.sol.u,tlist] = openfoam( fea, 'application', 'pimpleFoam', 'maxCo', 1.0, 'ddtSchemes', ddtSchemes, ...
'deltaT', opt.dt, 'endTime', 8, 'nproc', 1, 'fid', fid, 'logfid', logfid );
fea.sol.u = fea.sol.u(:,tlist > 0);
tlist = tlist(tlist > 0);
fid = logfid;
elseif( strcmp(opt.solver,'su2') )
logfid = fid; if( ~got.fid ), fid = []; end
[fea.sol.u,tlist] = su2( fea, 'fid', fid, 'logfid', logfid, 'tstep', opt.dt, 'tmax', 8, 'ischeme', opt.ischeme, 'wrtfreq', floor(8/opt.dt/400) );
fid = logfid;
else
[fea.sol.u,tlist] = solvetime( fea, 'fid', fid, 'tstep', opt.dt, 'tmax', 8, 'maxnit', 5, 'ischeme', opt.ischeme );
end
% Benchmark quantities.
[ c_d, c_l, dp ] = calc_bench_quant( fea, 1 );
[c_d_max, i] = max( c_d );
t_c_d_max = tlist( i );
[c_l_max, i] = max( c_l );
t_c_l_max = tlist( i );
out.c_d = c_d;
out.c_l = c_l;
out.dp = dp;
out.c_d_max = c_d_max;
out.c_l_max = c_l_max;
out.t_c_d_max = t_c_d_max;
out.t_c_l_max = t_c_l_max;
if( opt.instatbc )
[~, i] = min(abs(tlist-8));
dp_t8 = dp(i);
out.dp_t8 = dp_t8;
comp_data = [ t_c_d_max c_d_max t_c_l_max c_l_max dp_t8 ];
ref_data = [ 3.93625 2.950923849 5.6925 0.47834818 -0.11162153 ];
out.err = abs(ref_data-comp_data)./abs(ref_data);
out.pass = all( out.err < opt.tol );
else
found = false;
for j=i+1:length(c_l)
if( j==length(c_l) )
break
end
if( j>1 && c_l(j) > c_l(j-1) && c_l(j) > c_l(j+1) )
found = true;
c_l_max2 = c_l(j);
t_c_l_max2 = tlist(j);
break
end
end
for j=i-1:-1:1
if( j>1 && c_l(j) > c_l(j-1) && c_l(j) > c_l(j+1) )
found = true;
c_l_max2 = c_l(j);
t_c_l_max2 = tlist(j);
break
end
end
f = 1/(t_c_l_max2-t_c_l_max);
St = 0.1*f/1.5;
out.St = St;
t_dp = t_c_l_max + 1/2/f;
[~,j] = find( tlist > t_dp, 1 );
dp_f2 = dp(j-1) + (t_dp - tlist(j-1))/(tlist(j) - tlist(j-1)) * (dp(j)-dp(j-1));
out.dp_f2 = dp_f2;
c_d(tlist<0.5) = 0;
[c_d_max, i] = max( c_d );
comp_data = [ c_d_max c_l_max St dp_f2 ];
ref_data = [ 3.23 1.00 0.3 2.48 ];
out.err = abs(ref_data-comp_data)./abs(ref_data);
out.pass = all( out.err([1,2,4]) < opt.tol );
end
if( ~isempty(fid) )
fmtf = '%12.8f |';
fmts = '%12s |';
fmt = ['| ',repmat(fmtf,1,4+double(opt.instatbc)),'\n'];
fmts = ['| ',repmat(fmts,1,4+double(opt.instatbc)),'\n'];
fmtl = ['|------',repmat('-------------+',1,4+double(opt.instatbc))];
fmtl = [fmtl(1:end-1),'|\n'];
fprintf( fid, '\n\n' );
fprintf( fid, fmtl );
if( opt.instatbc )
fprintf( fid, fmts, 't(cd_max)', 'cd_max', 't(cl_max)', 'cl_max', 'dp(t=8)' );
fprintf( fid, fmtl );
fprintf( fid, fmt, t_c_d_max, c_d_max, t_c_l_max, c_l_max, dp_t8 );
else
fprintf( fid, fmts, 'cd_max', 'cl_max', 'St', 'dp' );
fprintf( fid, fmtl );
fprintf( fid, fmt, c_d_max, c_l_max, St, dp_f2 );
end
fprintf( fid, fmtl );
fmt = ['| Ref. ',repmat(fmtf,1,4+double(opt.instatbc)),'\n'];
fprintf( fid, fmt, ref_data );
fprintf( fid, fmtl );
fprintf( fid, '\n\n' );
end
% Postprocessing.
if( opt.iplot>0 )
figure
fea = parseprob( fea );
subplot(2,2,1)
postplot( fea, 'surfexpr', 'sqrt(u^2+v^2)', 'arrowexpr', {'u' 'v'}, 'arrowcolor', 'k' )
title('Velocity field at t=8')
ix = tlist > 1;
if( any(ix) )
tlist = tlist(ix);
c_d = c_d(ix);
c_l = c_l(ix);
dp = dp(ix);
end
subplot(2,2,3)
plot( tlist, c_d )
title('drag coefficient')
subplot(2,2,4)
plot( tlist, c_l )
title('lift coefficient')
subplot(2,2,2)
plot( tlist, dp )
title('pressure difference')
end
if( nargout==0 )
clear fea out
end
%