FEATool Multiphysics
v1.17.0
Finite Element Analysis Toolbox
|
EX_COMPRESSIBLEEULER1 SOD Shock tube example.
[ FEA, OUT ] = EX_COMPRESSIBLEEULER1( VARARGIN ) Sets up and solves an instationary SOD shock tube compressible Euler equation example on the unit line. Initial data are rhoL, uL=0, pL=1 and rhoR=0.125, uR=0, pR=0.1 with an inital discontinuity at t=0 and x=0.5. The computed solution with shocks, expansion, and compression waves is compared with the analytical Riemann solution.
[1] G.A. Sod, A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws, JCP, 1978.
Accepts the following property/value pairs.
Input Value/{Default} Description ----------------------------------------------------------------------------------- hmax scalar {0.002} Max grid cell size tmax scalar {0.2} Maximum time tstep scalar {0.0025} Time step size ischeme scalar {2} Time discretization scheme sfun string {sflag1} Shape function iplot scalar 0/{1} Plot solution and error (=1) . Output Value/(Size) Description ----------------------------------------------------------------------------------- fea struct Problem definition struct out struct Output struct
cOptDef = { 'hmax', 0.002; 'tmax', 0.2; 'tstep', 0.0025; 'sfun', 'sflag1'; 'ischeme', 2; 'iplot', 1; 'tol', 0.05; 'fid', 1 }; [got,opt] = parseopt(cOptDef,varargin{:}); fid = opt.fid; fea.sdim = { 'x' }; fea.geom.objects = { gobj_line(0,1) }; fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid ); fea = addphys(fea,@compressibleeuler); rhoL = 1; uL = 0; pL = 1; rhoR = 0.125; uR = 0; pR = 0.1; x = fea.grid.p(:); x0 = (max(x) - min(x))/2; init0 = { sprintf('%g+(x>%g)*%g',rhoL,x0,rhoR-rhoL), ... sprintf('%g+(x>%g)*%g',uL,x0,uR-uL), ... sprintf('%g+(x>%g)*%g',pL,x0,pR-pL) }; fea = parsephys(fea); fea = parseprob(fea); [fea.sol.u,fea.sol.t] = solvetime( fea, 'init', init0, 'fid', fid, ... 'ischeme', 2, 'tstep', opt.tstep, 'tmax', opt.tmax ); % Error checking. [r_ref,u_ref, p_ref] = analytic_sod( opt.tmax, x ); r = evalexprp( fea.dvar{1}, fea ); u = evalexprp( fea.dvar{2}, fea ); p = evalexprp( fea.dvar{3}, fea ); out.err = [ sum(abs(r_ref-r))/length(x), ... sum(abs(u_ref-u))/length(x), ... sum(abs(p_ref-p))/length(x) ]; out.pass = all(out.err<opt.tol); % Postprocessing. if( opt.iplot>0 ) figure subplot(2,2,1) plot( x, r, 'k-' ) hold on plot( x, r_ref, 'b--' ) legend('Computed','Analytic solution') title('Density (rho)') subplot(2,2,2) plot( x, u, 'k-' ) hold on plot( x, u_ref, 'b--' ) title('Velocity (u)') subplot(2,2,3) plot( x, p, 'k-' ) hold on plot( x, p_ref, 'b--' ) title('Pressure (p)') subplot(2,2,4) plot( x, p./r, 'k-' ) hold on plot( x, p_ref./r_ref, 'b--' ) title('Non-dimensional temperature (p/rho)') end if( nargout==0 ) clear fea out end %