|
FEATool Multiphysics
v1.17.2
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.2
Finite Element Analysis Toolbox
|
EX_NONNEWTONIAN1 2D Example for non-Newtonian flow in a channel.
[ FEA, OUT ] = EX_NONNEWTONIAN1( VARARGIN ) Sets up and solves stationary non-Newtonian flow in a rectangular channel. The fluid is governed by a power-law model such that the viscosity can be expressed as mu_eff = mu0*K*D^(n-1). A constant pressure drop is prescribed resulting in an analytical solution for the velocity u(y) = n/(n+1)*(1/(mu0*K)*dp)^(1/n)*((h/2)^((n+1)/n)-abs(h/2-y)^((n+1)/n)). Accepts the following property/value pairs.
Input Value/{Default} Description
-----------------------------------------------------------------------------------
rho scalar {1} Density
mu0 scalar {1.29684} Newtonian viscosity
dp scalar {1} Pressure drop
n scalar {0.25} Power-law exponent
K scalar {1} Power-law constant
h scalar {1} Channel height
l scalar {0.2} Channel length
igrid scalar 1/{0} Cell type (0=quadrilaterals, 1=triangles)
hmax scalar {0.1} Max grid cell size
sf_u string {sflag1} Shape function for velocity
sf_p string {sflag1} Shape function for pressure
iplot scalar 0/{1} Plot solution and error (=1)
.
Output Value/(Size) Description
-----------------------------------------------------------------------------------
fea struct Problem definition struct
out struct Output struct
cOptDef = { ...
'rho', 1;
'mu0', 1.29684;
'K', 1;
'n' , 0.25;
'dp' , 1;
'h', 1;
'l', 2;
'igrid', 1;
'hmax', 0.05;
'sf_u', 'sflag2';
'sf_p', 'sflag1';
'iplot', 1;
'fid', 1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid = opt.fid;
% Geometry and grid parameters.
h = opt.h; % Height of rectangular domain.
l = opt.l; % Length of rectangular domain.
% Discretization parameters.
sf_u = opt.sf_u; % FEM shape function type for velocity.
sf_p = opt.sf_p; % FEM shape function type for pressure.
% Geometry definition.
gobj = gobj_rectangle( 0, l, 0, h );
fea.geom.objects = { gobj };
fea.sdim = { 'x' 'y' }; % Coordinate names.
% Grid generation.
if ( opt.igrid==1 )
fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
else
fea.grid = rectgrid(round(l/opt.hmax),round(h/opt.hmax),[0 l;0 h]);
if( opt.igrid<0 )
fea.grid = quad2tri( fea.grid );
end
end
n_bdr = max(fea.grid.b(3,:)); % Number of boundaries.
% Boundary conditions.
dtol = opt.hmax;
i_inflow = findbdr( fea, ['x<',num2str(dtol)] ); % Inflow boundary number.
i_outflow = findbdr( fea, ['x>',num2str(l-dtol)] ); % Outflow boundary number.
% Problem definition.
fea = addphys(fea,@nonnewtonian); % Add non-Newtonian flow physics mode.
fea.phys.nn.eqn.coef{1,end} = { opt.rho };
fea.phys.nn.eqn.coef{2,end} = { 1 }; % Select power-law viscosity model.
fea.phys.nn.eqn.coef{4,end} = { opt.mu0 };
fea.phys.nn.eqn.coef{6,end} = { opt.K };
fea.phys.nn.eqn.coef{7,end} = { opt.n };
fea.phys.nn.eqn.coef{13,end}{1} = [num2str(opt.dp),'*(1-x/',num2str(l),')']; % Initial pressure.
fea.phys.nn.sfun = { sf_u sf_u sf_p }; % Set shape functions.
fea.phys.nn.bdr.sel([i_inflow,i_outflow]) = 4;
fea.phys.nn.bdr.coef{4,end}{3,i_inflow} = l/opt.dp; % Set inflow pressure.
% Parse and solve problem.
fea = parsephys(fea); % Check and parse physics mode.
fea = parseprob(fea); % Check and parse problem struct.
[fea.bdr.n{1}{[i_inflow,i_outflow]}] = deal( '-p*nx' ); % Offset natural Neumann BC.
[fea.bdr.d{2}{[i_inflow,i_outflow]}] = deal(0); % Set v=0 at inlet and outlet.
fea.sol.u = solvestat( fea, 'relchg', false, 'fid', fid, 'maxnit', 50 ); % Call to stationary solver.
% Postprocessing.
s_velm = 'sqrt(u^2+v^2)';
s_refsol = [num2str(opt.n/(opt.n+1)*(1/(opt.mu0*opt.K)*opt.dp)^(1/opt.n)),'*(',num2str((h/2)^((opt.n+1)/opt.n)),'-abs(',num2str(h/2),'-y)^',num2str((opt.n+1)/opt.n),')'];
p = [ l/2*ones(1,25); linspace(0,h,25) ];
u = evalexpr( s_velm, p, fea );
u_ref = evalexpr( s_refsol, p, fea );
s_velm = 'sqrt(u^2+v^2)';
if ( opt.iplot>0 )
figure
subplot(3,1,1)
postplot(fea,'surfexpr',s_velm,'evaltype','exact')
title('Velocity field')
subplot(3,1,2)
postplot(fea,'surfexpr','miu_nn')
title('Effective viscosity')
subplot(3,1,3)
plot( u, p(2,:) )
hold on
plot( u_ref, p(2,:), 'k.' )
legend( 'Computed solution', 'Reference solution','Location','West')
title( 'Velocity profile at x=l/2' )
ylabel( 'y' )
end
% Error checking.
err = sqrt(sum((u-u_ref).^2)/sum(u_ref.^2));
if( ~isempty(fid) )
fprintf(fid,'\nL2 Error: %f\n',err)
fprintf(fid,'\n\n')
end
out.err = err;
out.pass = err<0.05;
if ( nargout==0 )
clear fea out
end