FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
ex_fluidstructure4.m File Reference

## Description

EX_FLUIDSTRUCTURE4 Fluid-structure interaction - elastic beam.

[ FEA, OUT ] = EX_FLUIDSTRUCTURE1( VARARGIN ) 3D example for fluid-structure interaction flow around an elastic beam.

Reference

[1] Wall W., Ramm E. Fluid-Interaktion mit stabilisierten Finiten Elementen. Phd Thesis, Institut fur Baustatik, Universitat Stuttgart, 1999.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
iplot       scalar 0/{1}           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { 'iplot',    1;
'igrid',    1;
'tstep',    0.01;
'tmax',     1;
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});

% Geometry.
fea.sdim = { 'x', 'y', 'z' };
gobj1 = gobj_block( 0, 19.5, 0, 12, 0, 0.2, 'B1' );
gobj2 = gobj_block( 3.5, 13.5, 3.5, 8.5, 0, 0.2, 'B2' );
gobj3 = gobj_block( 4.5, 5.5, 5.5, 6.5, 0, 0.2, 'B3' );
gobj4 = gobj_block( 5.5, 9.5, 6-0.03, 6+0.03, 0, 0.2, 'B4' );
fea.geom.objects = { gobj1, gobj2, gobj3, gobj4 };
fea.geom = copy_geometry_object( 'B2', fea.geom );
fea.geom = geom_apply_formula( fea.geom, 'B1-B5' );
fea.geom = copy_geometry_object( 'B4', fea.geom );
fea.geom = geom_apply_formula( fea.geom, 'B2-B6-B3' );

% Grid generation.
if( opt.igrid )
hmax = [ 0.06 0.8 0.4 ]*0.95;
fea.grid = gridgen( fea, 'hmax', hmax, 'gridgen', 'gmsh', 'fid', opt.fid );
else
fea.grid = impexp_gmsh( 'ex_fluidstructure4.msh', 'import' );
end

% Equation and boundary settings.
fea = addphys( fea, @fluidstructure );

rho_f = 1.18e-3;
miu_f = 1.82e-4;
rho_s = 0.1;
nu_s  = 0.25;
E_s   = 2.5e6;

dtol  = num2str(sqrt(eps)*1e3);
ind_s = unique(fea.grid.s(selcells(fea.grid,['(y>=6-0.03-',dtol,').*(y<=6+0.03+',dtol,').*(x>=5.5-',dtol,').*(x<=9.5+',dtol,')'])));
ind_f = setdiff( unique(fea.grid.s), ind_s );

[fea.phys.fsi.eqn.coef{1,end}{ind_f}] = deal(rho_f);
[fea.phys.fsi.eqn.coef{1,end}{ind_s}] = deal(rho_s);
[fea.phys.fsi.eqn.coef{2,end}{ind_f}] = deal(miu_f);
[fea.phys.fsi.eqn.coef{3,end}{ind_s}] = deal(nu_s);
[fea.phys.fsi.eqn.coef{4,end}{ind_s}] = deal(E_s);
fea.phys.fsi.prop.active = zeros(7,length(ind_s)+length(ind_f));
fea.phys.fsi.prop.active(1:4,ind_f) = 1;
fea.phys.fsi.prop.active(5:7,ind_s) = 1;

n_bdr  = max(fea.grid.b(3,:));
i_in   = findbdr( fea, ['x<=',dtol] );
i_out  = findbdr( fea, ['x>=19.5-',dtol] );
i_slip = findbdr( fea, ['y>=12-',dtol,'|y<=',dtol] );
i_per  = [ findbdr( fea, ['z>=0.2-',dtol] ), findbdr( fea, ['z<=',dtol] ) ];
i_fix  = findbdr( fea, ['y>=6-0.03-',dtol,'&y<=6+0.03+',dtol,'&x<=5.5+',dtol] );
i_beam = setdiff( findbdr( fea, ['y>=6-0.03-',dtol,'&y<=6+0.03+',dtol] ), [i_fix,i_per] );
i_ns   = setdiff( findbdr( fea, ['x>=4.5-',dtol,'&x<=5.5+',dtol] ), i_fix );
i_cont = setdiff( 1:max(n_bdr), [i_in,i_out,i_slip,i_per,i_beam,i_ns,i_fix] );

sel = ones(1,n_bdr);
sel( i_in )   =  2;
sel( i_out )  =  3;
sel( i_slip ) =  5;
sel( i_per )  =  5;
sel( i_fix )  =  6;
sel( i_beam ) = -2;
sel( i_cont ) = -1;

fea.phys.fsi.bdr.sel = sel;
fea.phys.fsi.bdr.coef{2,end}{1,i_in} = 51.3;

% Solver.
fea = parsephys(fea);
fea = parseprob(fea);

[fea.sol.u,fea.sol.t,fea.sol.grid.p] = fsisolve( fea, 'tstep', opt.tstep, 'tmax', opt.tmax, 'fid', opt.fid );

% Postprocessing.
if( opt.iplot )
postplot(fea,'slicex',[],'slicey',[],'sliceexpr','sqrt(u^2+v^2+w^2)')
rotate3d on
end

% Error checking.
v_min = inf; v_max = -inf;
for i=80:size(fea.sol.u,2)
[vm,vx] = minmaxsubd( 'dy', fea, [], [], [], i );
if( vm<v_min )
v_min = vm;
i_min = i;
end
if( vx>v_max )
v_max = vx;
i_max = i;
end
end

out.err  = [ abs((v_min + 0.5)/0.5), abs((v_max - 0.45)/0.45) ];
out.pass = all( out.err < 0.5 );

if( nargout==0 )
clear fea out
end