First version of code

MOC_2D_steady_irrotational_IVLINE.m

@@ -0,0 +1,30 @@
+function [xsonic,xvnull,u] = MOC_2D_steady_irrotational_IVLINE ( geom , params , y )
+
+  % Origin of the coordinate system relative to the nozzle throat [m]
+  eps = -geom.yt * sqrt((params.gamma+1)*(geom.delta+1)*geom.yt/geom.rhou) * 0.5 / (3+geom.delta) ;
+
+  % Sonic speed a*
+  astar = sqrt(2*params.gamma*params.R*params.T/(params.gamma+1)) ;
+
+  % Coefficient of the linear axial perturbation velocity [1/m]
+  alpha = sqrt((1+geom.delta)/(geom.yt*geom.rhou*(params.gamma+1)));
+
+  % Equation for the sonic line x = coeff_sonic * y^2, where coeff_sonic is [1/m]
+  coeff = -(params.gamma+1)*alpha*0.5/(1+geom.delta);
+  xsonic=coeff*y.^2 - eps;
+
+  % Equation for the v=0 line: x=coeff_vnull * y^2, where coeff_vnull is [1/m]
+  coeff_vnull = -(params.gamma+1)*alpha*0.5/(3+geom.delta);
+  xvnull=coeff_vnull*y.^2 - eps;
+
+  % Perturbation velocity field on the v=0 line
+  u = astar * (1 + alpha * (xvnull+eps) + (params.gamma+1)*(alpha^2)*(y.^2)*0.5/(1+geom.delta) );
+
+endfunction
+
+function a = get_speed_sound(params,V)
+% Speed of sound a² = a0²       - 0.5*(gamma-1)*V²
+%                   = gamma*R*T - 0.5*(gamma-1)*V²
+  a = sqrt ( params.gamma * params.R * params.T - ...
+             0.5 * ( params.gamma - 1.) * V.^2 ) ;
+endfunction

MOC_2D_steady_irrotational_get_geometry.m

@@ -0,0 +1,33 @@
+function [abc,y,tangent] = MOC_2D_steady_irrotational_get_geometry(x,step,geom)
+% Get the geometry of a nozzle, see page 128, Gas Dynamics V2
+
+  radtu = geom.rhou ; %0.05   ; % Radius of curvature of upstream   circular arc
+  radtd = geom.rhod ; %0.0125 ; % Radius of curvature of downstream circular arc
+  yt    = geom.yt   ; %0.025  ; % Radius of nozzle throat
+  thetaa= geom.ta   ; %35     ; % Attachment angle between downstream arc and 2nd order poly
+  thetae= geom.te   ; %10     ; % Nozzle exit lip angle
+  xe    = geom.xe   ; %0.25   ; % Nozzle length
+
+  % Coordinates of attachment point between downstream arc and 2nd order poly
+  xa    = radtd * sind ( thetaa ) ;
+  ya    = yt + ( 1 - cosd ( thetaa ) ) * radtd ;
+
+  % 2nd order poly       y = a + bx + cx^2
+  %                tangent = b + 2cx
+  matA = [ 1 , xa , xa^2 ;
+           0 , 1  , 2*xa ;
+           0 , 1  , 2*xe  ] ;
+  RHS = [ ya ;
+          tand(thetaa) ;
+          tand(thetae) ] ;
+  abc = matA\RHS ;
+
+  y       = 0.;
+  tangent = 0.;
+  if (step==2) % get the geometry after the downstream circular arc
+    assert ( x >= xa );
+%    assert ( x <= xe );
+    y       = abc(1) + abc(2) * x + abc(3) * x.^2 ;
+    tangent = abc(2) + 2*abc(3) * x ;
+  endif
+endfunction

MOC_2D_steady_irrotational_get_thermo.m

@@ -0,0 +1,10 @@
+function [Mach,pressure,density,temperature,sound] = MOC_2D_steady_irrotational_get_thermo(u,v,params)
+% Get the thermodynamic data from the velocity components
+  CP          = params.gamma * params.R / ( params.gamma - 1 );
+  Q           = sqrt( u.^2 + v.^2 ) ;
+  temperature = params.T - 0.5 * Q.^2 / CP ;
+  sound       = sqrt ( params.gamma * params.R * temperature ); 
+  Mach        = Q ./ sound ;
+  pressure    = params.P .* ( temperature / params.T ).^(params.gamma/(params.gamma-1)) ;
+  density     = pressure ./ ( params.R * temperature ) ;
+endfunction

MOC_2D_steady_irrotational_internal_point.m

@@ -0,0 +1,85 @@
+function [xo,yo,uo,vo] = MOC_2D_steady_irrotational_internal_point ( xp,yp,up,vp,...
+                                                                     xm,ym,um,vm,...
+                                                                     params )
+% xp,yp,up,vp are conditions at left -running characteristic C+
+% xm,ym,um,vm are conditions at right-running characteristic C-
+   x_orig = [ xp , xm ] ; y_orig = [ yp , ym ] ;
+   u_orig = [ up , um ] ; v_orig = [ vp , vm ] ;
+
+   step_current = 0;
+   step_max = 50;
+   eps_pos = 1.e-8;
+   eps_vel = 1.e-5;
+
+% Predictor step
+   [xo,yo,uo,vo] = MOC_2D_steady_irrotational_solve_internal_point ( xp,yp,up,vp,...
+                                                                     xm,ym,um,vm,...
+                                                                     x_orig,y_orig,u_orig,v_orig,...
+                                                                     params ) ;
+   while 1
+      step_current++;
+
+% Corrector step
+      xcp = xp; ycp = 0.5*(yp+yo) ; ucp = 0.5*(up+uo) ; vcp = 0.5*(vp+vo) ;
+      xcm = xm; ycm = 0.5*(ym+yo) ; ucm = 0.5*(um+uo) ; vcm = 0.5*(vm+vo) ;
+     [xn,yn,un,vn] = MOC_2D_steady_irrotational_solve_internal_point ( xcp,ycp,ucp,vcp,...
+                                                                       xcm,ycm,ucm,vcm,...
+                                                                       x_orig,y_orig,u_orig,v_orig,...
+                                                                       params ) ;
+      error_pos = max( [ xo-xn , yo-yn ] );
+      error_vel = max( [ uo-un , vo-vn ] );
+      xo = xn ;
+      yo = yn ;
+      uo = un ;
+      vo = vn ;
+      if (abs(error_pos)<eps_pos && abs(error_vel)<eps_vel)
+        break;
+      endif
+      if (step_current>step_max)
+        error('The maximum of iterations for the predictor-corrector algorithm has been reached. Stopping the execution...')
+      endif
+   end
+
+endfunction
+
+function [xn,yn,un,vn] = MOC_2D_steady_irrotational_solve_internal_point ( xp,yp,up,vp,...
+                                                                           xm,ym,um,vm,...
+                                                                           x_orig,y_orig,u_orig,v_orig,...
+                                                                           params )
+
+  xo = [xp xm] ;
+  yo = [yp ym] ;
+  uo = [up um] ;
+  vo = [vp vm] ;
+
+  V     = sqrt  ( uo.^2 + vo.^2 ) ;
+  a     = get_speed_sound( params, V ) ;
+  theta = atand (vo./uo) ;
+  alpha = asind ( a./V ) ;
+  lambda(:,1) = tand ( theta(1) + alpha(1) );  % lambda+
+  lambda(:,2) = tand ( theta(2) - alpha(2) );  % lambda-
+  Q     = uo.^2 - a.^2 ;
+  R     = 2*uo.*vo - Q.*lambda ;
+  S     = a.^2 .* vo ./ yo ;
+
+% Solve the system matrix to get the position of the intersection of the C+ and C- characteristics
+  matA    = [ -lambda(1) 1 ; -lambda(2) 1 ] ;
+  vecB    = [ y_orig(1)-lambda(1)*x_orig(1) ; y_orig(2)-lambda(2)*x_orig(2) ] ;
+  new_pos = matA\vecB ;
+  xn   = new_pos(1);
+  yn   = new_pos(2);
+
+  T       = S(:).*(xn-x_orig(:)) + Q(:).*u_orig(:) + R(:).*v_orig(:) ;
+  matA    = [ Q(1) , R(1) ; Q(2) , R(2) ] ;
+  new_vel = matA\T ;
+  un   = new_vel(1) ;
+  vn   = new_vel(2) ;
+
+endfunction
+
+function a = get_speed_sound(params,V)
+% Speed of sound a² = a0²       - 0.5*(gamma-1)*V²
+%                   = gamma*R*T - 0.5*(gamma-1)*V²
+  a = sqrt ( params.gamma * params.R * params.T - ...
+             0.5 * ( params.gamma - 1.) * V.^2 ) ;
+endfunction