Description

EX_COMPRESSIBLEEULER5 3D Compressible inviscid flow over a bump.

[ FEA, OUT ] = EX_COMPRESSIBLEEULER5( VARARGIN ) Sets up and solves a stationary 3D compressible Euler equation problem for supersonic flow over a cylindrical bump.

Reference

[1] Lynn J.F., van Leer B., Lee D. (1997) Multigrid solution of the euler equations with local preconditioning. In: Kutler P., Flores J., Chattot JJ. (eds) Fifteenth International Conference on Numerical Methods in Fluid Dynamics. Lecture Notes in Physics, vol 490. Springer, Berlin, Heidelberg.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.05}          Max grid cell size
sfun        string {sflag1}        Shape function
solver      string openfoam/su2/{} Use OpenFOAM, SU2, or default solver
iplot       scalar 0/{1}           Plot solution and error (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { 'hmax',     0.08;
             'sfun',     'sflag1';
             'solver',   '';
             'iplot',    1;
             'tol',      0.11;
             'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


 fea.sdim = { 'x', 'y', 'z' };
 r = (0.5^2/0.042+0.042)/2;
 fea.geom.objects = { gobj_block(-1,4,0,0.5,0,2), gobj_cylinder([.5,0,0.042-r],r,.5,[0,1,0]) };

 fea = geom_apply_formula(fea,'B1-C1');
 fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);

 fea = addphys(fea,@compressibleeuler);

 Ma    = 1.4;
 rho0  = 1;
 p0    = 1;
 u0    = Ma*sqrt(1.4*p0/rho0);
 v0    = 0;
 w0    = 0;
 init0 = {rho0, u0, v0, w0, p0};
 fea.phys.ee.eqn.coef{6,end}{1} = rho0;
 fea.phys.ee.eqn.coef{7,end}{1} = u0;
 fea.phys.ee.eqn.coef{8,end}{1} = v0;
 fea.phys.ee.eqn.coef{9,end}{1} = w0;
 fea.phys.ee.eqn.coef{10,end}{1} = p0;

 i_in = findbdr( fea, 'x<=-1+sqrt(eps)' );
 i_out = findbdr( fea, 'x>=4-sqrt(eps)' );
 fea.phys.ee.bdr.sel(i_in) = 1;
 fea.phys.ee.bdr.sel(i_out) = 2;
 fea.phys.ee.bdr.coef{1,end}{1,i_in} = rho0;
 fea.phys.ee.bdr.coef{1,end}{2,i_in} = u0;
 fea.phys.ee.bdr.coef{1,end}{3,i_in} = v0;
 fea.phys.ee.bdr.coef{1,end}{4,i_in} = w0;
 fea.phys.ee.bdr.coef{1,end}{5,i_in} = p0;

 fea = parsephys(fea);
 fea = parseprob(fea);


 if( strcmp(opt.solver,'openfoam') )
   logfid = fid; if( ~got.fid ), fid = []; end
   fea.sol.u = openfoam( fea, 'deltaT', 0.01, 'endTime', 5, 'fid', fid, 'logfid', logfid );
   fid = logfid;
 elseif( strcmp(opt.solver,'su2') )
   logfid = fid; if( ~got.fid ), fid = []; end
   fea.sol.u = su2( fea, 'upwind', 'jst', 'fid', fid, 'logfid', logfid, 'nproc', 1 );
   fid = logfid;
 else
   fea.sol.u = solvestat( fea, 'init', init0, 'maxnit', 30, 'nlrlx', 1.0, 'fid', fid );
 end



% Postprocessing.
 s_Ma = fea.phys.ee.eqn.vars{end-1,2};
 if( opt.iplot>0 )
   postplot( fea, 'sliceexpr', s_Ma, 'slicex', [], 'slicey', 0.45, 'slicez', 0.05 )
   title(fea.phys.ee.eqn.vars{end-1,1})
 end


% Error checking.
 [Ma_min,Ma_max] = minmaxsubd( s_Ma, fea );
 out.err = [ abs(Ma_min-1.0)/1.0, ...
             abs(Ma_max-1.8)/1.8 ];
 out.pass = all(out.err<opt.tol);

 if( nargout==0 )
   clear fea out
 end