概要
ラバルノズルとは先細末広タイプのノズルであり、超音速流れを加速するために使用される装置である。ノズル内の圧力分布は、以下の模式図に示すように出口圧力に依存して変化する。出口圧力が入口圧力よりもわずかに小さい場合、ノズル内は亜音速流れとなる(case-A)。出口圧力をさらに下げるとスロート部でマッハ数が１の臨界状態になり下流側では超音速流れになる。この場合、出口圧力が十分に小さければ下流側全体で超音速流れになる(case-C)が、そうでない場合、途中で垂直衝撃波が発生する流れとなる(case-B)。
ここでは、ラバルノズルの検証例題を対象として、理論解の導出とOpenFOAMによる数値解を示す。
理論解
一様エントロピー関係式
これは圧縮性流体の運動方程式から導かれる式であり、非圧縮性流体のベルヌイの式に対応している。流管上の任意の２点間には以下の関係式が成り立つ。ここで$\gamma$は比熱比、$M$はマッハ数、$p$は圧力、$T$は温度、$\rho$は密度、$u$は流速、$A$は流路断面積を表す。
\begin{align}
\frac{p_1}{p_2}&=\left[
\frac{1+\frac{\gamma-1}{2}M_1^2}
{1+\frac{\gamma-1}{2}M_2^2}
\right]^{\frac{-\gamma}{\gamma-1}} \tag{1} \\
\frac{T_1}{T_2}&=\left[
\frac{1+\frac{\gamma-1}{2}M_1^2}
{1+\frac{\gamma-1}{2}M_2^2}
\right]^{-1} \tag{2} \\
\frac{\rho_1}{\rho_2}&=\left[
\frac{1+\frac{\gamma-1}{2}M_1^2}
{1+\frac{\gamma-1}{2}M_2^2}
\right]^{\frac{-1}{\gamma-1}} \tag{3}\\
\frac{u_1}{u_2}&=\frac{M_1}{M_2}\left[
\frac{1+\frac{\gamma-1}{2}M_1^2}
{1+\frac{\gamma-1}{2}M_2^2}
\right]^{\frac{1}{2}} \tag{4}\\
\frac{A_1}{A_2}&=\frac{M_2}{M_1}\left[
\frac{1+\frac{\gamma-1}{2}M_1^2}
{1+\frac{\gamma-1}{2}M_2^2}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{5}
\end{align}
よどみ点状態と臨界状態の関係
臨界状態ではスロート部においてはマッハ数が１($M_*=1$)となる。前記の一様エントロピー関係式を使うと、よどみ点（添字0で表す)とスロート位置（添字*で表す)の間には以下の関係式が成り立つ。ここで$a$は音速を表す。
\begin{align}
\frac{p_*}{p_0}&=\left(
\frac{2}{\gamma+1}
\right)^{\frac{\gamma}{\gamma-1}} \tag{6} \\
\frac{T_*}{T_0}&=\frac{2}{\gamma+1} \tag{7} \\
\frac{\rho_*}{\rho_0}&=\left(
\frac{2}{\gamma+1}
\right)^{\frac{1}{\gamma-1}} \tag{8}\\
\frac{a_*}{a_0}&=\left(
\frac{2}{\gamma+1}
\right)^{\frac{1}{2}} \tag{9}
\end{align}
臨界状態と他の断面の物理量の関係
前記の一様エントロピー関係式を使うと、臨界状態のスロート（添字*で表す。$M_*=1$）と他の任意断面との間には以下の関係式が成り立つ。
\begin{align}
\frac{A}{A_*}&=\frac{1}{M}
\left[
\frac{1+\frac{\gamma-1}{2}M^2}{\frac{\gamma+1}{2}}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{10} \\
\frac{p}{p_*}&=\left[
\frac{\frac{\gamma+1}{2}}
{1+\frac{\gamma-1}{2}M^2}
\right]^{\frac{\gamma}{\gamma-1}} \tag{11} \\
\frac{T}{T_*}&=
\frac{\frac{\gamma+1}{2}}
{1+\frac{\gamma-1}{2}M^2} \tag{12} \\
\frac{\rho}{\rho_*}&=\left[
\frac{\frac{\gamma+1}{2}}
{1+\frac{\gamma-1}{2}M^2}
\right]^{\frac{1}{\gamma-1}} \tag{13} \\
\frac{u}{u_*}&=M \left[
\frac{\frac{\gamma+1}{2}}
{1+\frac{\gamma-1}{2}M^2}
\right]^{\frac{1}{2}} \tag{14}
\end{align}
(10)式を見ると臨界状態の流れでは、断面積が決まると自動的にマッハ数が規定されることがわかる（すなわち流路形状だけで任意箇所のマッハ数が決まるのである）。与えられた断面積$A$に対してこの関係式を満足する$M$は２つ存在する。これらはスロート部の上流側（M<1）、下流側(M>1)に対応している。
垂直衝撃波
静止している垂直衝撃波を考える。衝撃波の直前（断面1）、直後(断面2)の物理量には以下の関係式が成り立つ。
\begin{align}
\frac{p_2}{p_1}&=1+
\frac{2\gamma}{\gamma+1}(M_1^2-1) \tag{15}\\
\frac{u_2}{u_1}&=
\frac{2+(\gamma-1)M_1^2}{(\gamma+1)M_1^2} \tag{16}\\
\frac{\rho_2}{\rho_1}&=
\frac{(\gamma+1)M_1^2}{2+(\gamma-1)M_1^2} \tag{17}\\
\frac{T_2}{T_1}&=1+
\frac{2(\gamma-1)(\gamma M_1^2+1)(M_1^2-1)}{(\gamma+1)^2M_1^2} \tag{18}\\
M_2^2&=
\frac{1+\frac{\gamma-1}{2}M_1^2}
{\gamma M_1^2-\frac{\gamma-1}{2}} \tag{19}
\end{align}
理論解の導出
限界圧力値の計算
出口圧力$p_e$が、限界値$p_{c1}$よりも大きい場合、ノズル内部は全領域にわたり亜音速流れとなり、限界値$p_{c2}$よりも小さい場合、スロート下流部全体で超音速流れとなる。出口圧力$p_e$が$p_{c2}<p_e<p_{c1}$の場合、ノズル内部に衝撃波が発生することになる。これらの限界圧力値($p_{c1},p_{c2}$)は以下の式を使い求める($A_*,p_*$は既知であることに注意)。
\begin{align}
\frac{A_e}{A_*}&=\frac{1}{M_e}
\left[
\frac{1+\frac{\gamma-1}{2}M_e^2}{\frac{\gamma+1}{2}}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{10} \\
\frac{p_e}{p_*}&=\left[
\frac{\frac{\gamma+1}{2}}
{1+\frac{\gamma-1}{2}M_e^2}
\right]^{\frac{\gamma}{\gamma-1}} \tag{11}
\end{align}
(10)式を解いて得られる$M_e$は２つ存在し、それぞれ亜音速流れと超音速流れに対応している。この$M_e$を(11)式に代入することで限界圧力値が求められる。
流れ場の理論解
case-A
まず(1)式から導かれる次式を使い出口マッハ数$M_e$を計算する。
\begin{align}
\frac{p_0}{p_e}&=\left(
1+\frac{\gamma-1}{2}M_e^2
\right)^{\frac{\gamma}{\gamma-1}} \tag{1}
\end{align}
次に次式を使いノズル内任意位置(断面積$A$)でのマッハ数$M$を求め、各種の特性値を計算する。
\begin{align}
\frac{A}{A_e}&=\frac{M_e}{M}\left[
\frac{1+\frac{\gamma-1}{2}M^2}
{1+\frac{\gamma-1}{2}M_e^2}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{5}
\end{align}
case-B
衝撃波位置の物理量は、６個の未知変数($M_1,M_2,A_1(=A_2),p_1,p_2,M_e$)に対して６個の非線形方程式を連立して計算する必要がある。以下に方程式を再記しておく。なお添字$e$は出口を表している。なお、これらの解には衝撃波位置$x_*$が明示的に表現されていないが、$A(x_*)=A_1$から逆算する。
- スロートと断面1の間の等エントロピー関係式
\begin{align}
\frac{A_1}{A_*}&=\frac{1}{M_1}\left[
\frac{1+\frac{\gamma-1}{2}M_1^2}{\frac{\gamma+1}{2}}
\right]^{\frac{\gamma+1}{2(\gamma-1)}} \tag{10}\\
\frac{p_1}{p_*}&=\left[
\frac{\frac{\gamma+1}{2}}{1+\frac{\gamma-1}{2}M_1^2}
\right]^{\frac{\gamma}{\gamma-1}} \tag{11}
\end{align}
- 断面1と断面2の間の衝撃波関係式
\begin{align}
\frac{p_2}{p_1}&=1+
\frac{2\gamma}{\gamma+1}(M_1^2-1) \tag{15}\\
M_2^2 &= \frac{1+\frac{\gamma-1}{2}M_1^2}
{\gamma M_1^2 - \frac{\gamma-1}{2}} \tag{19}
\end{align}
- 断面2と出口の間の等エントロピー関係式
\begin{align}
\frac{p_2}{p_e}&=\left[
\frac{1+\frac{\gamma-1}{2}M_2^2}
{1+\frac{\gamma-1}{2}M_e^2}
\right]^{\frac{-\gamma}{\gamma-1}} \tag{1}\\
\frac{A_2}{A_e}&=\frac{M_e}{M_2}\left[
\frac{1+\frac{\gamma-1}{2}M_2^2}
{1+\frac{\gamma-1}{2}M_e^2}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{5}
\end{align}
流れ場の特性値は、扱う方程式が異なるため衝撃波前後に分けて計算を行う。
衝撃波上流側ではスロート位置での臨界条件およびよどみ点状態を境界条件として設定する。断面積-マッハ数関係式（上記 式(10)）から、マッハ数を求め一様エントロピー関係式から物性値を計算する。スロート前後において亜音速と超音速が切り替わるため、マッハ数の初期値を適切に指定しないと正しい解が得られないことに注意。
\begin{align}
\frac{A}{A_*}&=\frac{1}{M}
\left[
\frac{1+\frac{\gamma-1}{2}M^2}{\frac{\gamma+1}{2}}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{10}
\end{align}
衝撃波下流側では衝撃波直後の条件、および流出口圧力を境界条件として設定する。流出口温度（あるいは圧力）を計算するためにはノズル通過流量$\dot{m}$（これは臨界流量として計算される）を使用する。物性値の計算には、まず断面積-マッハ数関係式（上記 式(5)）からマッハ数を求め、一様エントロピー関係式を使用する。
\begin{align}
\frac{A}{A_e}&=\frac{M_e}{M}\left[
\frac{1+\frac{\gamma-1}{2}M^2}
{1+\frac{\gamma-1}{2}M_e^2}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{5}
\end{align}
\begin{align}
\dot{m}&= \rho_s a_s A_s \tag{20}\\
a_e&= \frac{\gamma p_e M_e A_e}{\dot{m}} \tag{21}\\
T_e &= \frac{a_e^2}{\gamma R} \tag{22} \\
u_e &= M_e a_e \tag{23} \\
\rho_e &= \frac{p_e}{R T_e} \tag{24}
\end{align}
case-C
次式を使いノズル内任意位置(断面積$A$)でのマッハ数$M$を求め、各種の特性値を計算する。スロート前後で亜音速と超音速が切り替わるため、非線形求解時の初期値に注意すること。
\begin{align}
\frac{A}{A_*}&=\frac{1}{M}
\left[
\frac{1+\frac{\gamma-1}{2}M^2}{\frac{\gamma+1}{2}}
\right]^{\frac{(\gamma+1)}{2(\gamma-1)}} \tag{10}
\end{align}
OpenFOAMによる数値解
解析モデル
解析モデルは以下の通りである。
メッシュ作成
OpenFOAMの軸対称モデルを使って計算を行う。ここでは blockMeshを使ってメッシュを作成する。
ノズル形状
ラバルノズルは先細-末広形状となっており形状を以下の式で表現する。
if ( x .lt. 5.0 ) then
area = 1.75 - 0.75 * cos( ( 0.2 * x - 1.0 ) * pi )
else
area = 1.25 - 0.25 * cos( ( 0.2 * x - 1.0 ) * pi )
endif
注：ここで $0 \le x \le 10 [inch]$。$A=\pi r^2$から半径$r$の分布を計算すると以下のようになる。
blockMeshDictの作成
以下のblockMeshDictを作成する。ノズルの壁面形状は codeStream で記述している。blockMeshDict内で定義されたパラメータをソースプログラム内で参照することができる（ここでは ndivを使用）。
system/blockMeshDict
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2312 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
mergeType points; // Wedge geometry - Merge points instead of topology
scale 0.0254;
L 10;
ang ${{ 2.5*3.14159265358979/180.0 }};
cs ${{ cos($ang) }};
sn ${{ sin($ang) }};
rin 0.892062058076386;
rout 0.690988298942671;
yin ${{ $rin*$cs }};
zin ${{ $rin*$sn }};
yout ${{ $rout*$cs }};
zout ${{ $rout*$sn }};
ndiv 200;
vertices
(
( 0 0 0) // 0
($L 0 0) // 1
( 0 $yin -$zin) // 2
($L $yout -$zout) // 3
( 0 $yin $zin) // 4
($L $yout $zout) // 5
);
blocks
(
hex (0 1 3 2 0 1 5 4) ($ndiv 20 1) simpleGrading (1 1 1)
);
edges #codeStream
{
codeInclude
#{
#include "pointField.H"
#include "mathematicalConstants.H"
#};
code
#{
constexpr scalar xMin = 0;
constexpr scalar xMax = 10;
constexpr label nPoints = $ndiv + 1;
constexpr scalar dx = (xMax - xMin)/scalar(nPoints - 1);
constexpr scalar pi = constant::mathematical::pi;
scalar cs = cos(2.5*pi/180.0);
scalar sn = sin(2.5*pi/180.0);
os << "(" << nl << "spline 2 3" << nl;
pointField profile(nPoints);
for (label i = 0; i < nPoints; ++i)
{
scalar x = xMin + i*dx;
scalar r(0.0);
if ( x < 5.0 ) {
scalar area ( 1.75 - 0.75 * cos( ( 0.2 * x - 1.0 ) * pi ) );
r = sqrt(area/pi);
}
else {
scalar area ( 1.25 - 0.25 * cos( ( 0.2 * x - 1.0 ) * pi ) );
r = sqrt(area/pi);
}
profile[i].x() = x;
profile[i].y() = r*cs;
profile[i].z() =-r*sn;
}
os << profile << nl;
os << "spline 4 5" << nl;
for (label i = 0; i < nPoints; ++i)
{
scalar x = xMin + i*dx;
scalar r(0.0);
if ( x < 5.0 ) {
scalar area ( 1.75 - 0.75 * cos( ( 0.2 * x - 1.0 ) * pi ) );
r = sqrt(area/pi);
}
else {
scalar area ( 1.25 - 0.25 * cos( ( 0.2 * x - 1.0 ) * pi ) );
r = sqrt(area/pi);
}
profile[i].x() = x;
profile[i].y() = r*cs;
profile[i].z() = r*sn;
}
os << profile << nl;
os << ");" << nl;
#};
};
boundary
(
inlet
{
type patch;
faces
(
(0 2 4 0)
);
}
outlet
{
type patch;
faces
(
(1 3 5 1)
);
}
wall
{
type wall;
faces
(
(2 4 5 3)
);
}
wedge1
{
type wedge;
faces
(
(0 2 3 1)
);
}
wedge2
{
type wedge;
faces
(
(0 1 5 4)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
pointの論理配置は以下のようになっている。軸対称モデルなので全体は楔形形状である(中心角=5[deg])。
blockMeshの実行
ケースディレクトリ内で blockMeshとタイプすることでメッシュが作成される。ただし、標準設定のままだと以下のような警告メッセージがたくさん出力される（メッシュはできているぽい）
--> FOAM Warning :
From virtual void Foam::wedgePolyPatch::calcGeometry(Foam::PstreamBuffers&)
in file meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C at line 73
Wedge patch 'wedge2' is not planar.
At local face at (0.253364 0.0170957 0.000746413) the normal (-7.04761e-07 -0.0436194 0.999048) differs from the average normal (-2.02391e-10 -0.0436194 0.999048) by 4.97371e-13
Either correct the patch or split it into planar parts
この警告を抑制するためには、 system/controlDict の writePrecisionの数値を大きくする。例えば
writePrecision 10; // increased 6 -> 10
最終的なメッシュは以下のようになる。
流入境界付近の拡大図を以下に示す。
解析ケース
出口圧力を変えて３ケース実施した。ここで$p_e$は出口圧力、$p_t$は入口（よどみ点）圧力である。
| case | 圧力比($p_e/p_t$) | 特徴 |
|---|---|---|
| A | 0.89 | 亜音速 |
| B | 0.75 | 亜音速->超音速->(衝撃波)->亜音速 |
| C | 0.16 | 亜音速->超音速 |
解析結果
解析には、rhoCentralFoam（非定常圧縮性流体ソルバー）を使用し解が定常状態に漸近するまで計算を行った。各ケースの結果を以下に示す。なおxy図には比較のため理論解と数値解を表示している。
case-A
流体特性値を計算するjuliaプログラムを以下に示す。
profile.jl
using NLsolve
gasconst = 8.31446262 # m2 kg s-2 K-1 mol-1
MW = 29e-3
cp = 1005
R = gasconst / MW
gamma = cp/(cp-R)
A0 = 0.0016129 # =2.5[inch^2]
p0 = 6895
T0 = 100
rho0 = p0/(R*T0)
a0 = sqrt( gamma*R*T0 )
Ts = 2/(gamma+1) * T0
ps = (2/(gamma+1))^(gamma/(gamma-1)) * p0
rhos = (2/(gamma+1))^(1/(gamma-1)) * rho0
as = sqrt( 2/(gamma+1) ) * a0
us = as
As = 0.00064516 # =1.0[inch^2]
mdot = rhos*us*As
Ae = 0.00096774 # =1.5[inch^2]
pe = 6137
f(x,gamma,p0,pe) = ( p0/pe
- ( (1+(gamma-1)/2*x[1]^2) )^(gamma/(gamma-1)) )
g(x) = f(x,gamma,p0,pe)
sol = nlsolve( g, [ 0.3 ])
result = sol.zero
Me = result[1]
Te = T0 / ( (1+(gamma-1)/2*Me^2) )
rhoe = rho0 / ( (1+(gamma-1)/2*Me^2)^(1/(gamma-1)) )
ae = sqrt( R*gamma*Te )
ue = Me*ae
n = 100
for i in 0:n
x = i*(0.254/n)
xin = i*(10/n)
area = ( x<=0.127 ? (1.75-0.75*cos((0.2*xin-1.0)*pi) )*2.54e-2^2
: (1.25-0.25*cos((0.2*xin-1.0)*pi) )*2.54e-2^2 )
f(x,area,gamma,Me,Ae) = ( area/Ae
- ( (1+(gamma-1)/2*x[1]^2)/(1+(gamma-1)/2*Me^2) )^((gamma+1)/(2*(gamma-1)))*Me/x[1] )
g(x) = f(x,area,gamma,Me,Ae)
local sol = nlsolve( g, [ 0.3 ], show_trace=false )
local result = sol.zero
M = result[1]
rho = ( (1+(gamma-1)/2*M^2)/(1+(gamma-1)/2*Me^2) )^(-1/(gamma-1)) * rhoe
T = ( (1+(gamma-1)/2*M^2)/(1+(gamma-1)/2*Me^2) )^(-1) * Te
p = ( (1+(gamma-1)/2*M^2)/(1+(gamma-1)/2*Me^2) )^(-gamma/(gamma-1)) * pe
u = M/Me*( (1+(gamma-1)/2*Me^2)/(1+(gamma-1)/2*M^2) )^(0.5) * ue
println("$x $area $M $p $T $rho $u")
end
case-B
衝撃波条件（6変数の連立非線形方程式）を計算するjuliaプログラムを以下に示す。
shock.jl
gasconst = 8.31446262 # m2 kg s-2 K-1 mol-1
MW = 29e-3
cp = 1005
R = gasconst / MW
gamma = cp/(cp-R)
A0 = 0.0016129
p0 = 6895
T0 = 100
rho0 = p0/(R*T0)
a0 = sqrt( gamma*R*T0 )
As = 0.00064516
Ts = 2/(gamma+1) * T0
ps = (2/(gamma+1))^(gamma/(gamma-1)) * p0
rhos = (2/(gamma+1))^(1/(gamma-1)) * rho0
as = sqrt( 2/(gamma+1) ) * a0
Ae = 0.00096774
pe = 5171
#####################################################
using NLsolve
function f(x,gamma,As,Ae,p0,pe,ps)
# 1 2 3 4 5 6
# M1 M2 A1 p1 p2 Me
[ x[3]/As - ( (1+(gamma-1)/2*x[1]^2)/((gamma+1)/2) )^((gamma+1)/(2*(gamma-1)))/x[1],
x[4]/ps - ( ((gamma+1)/2)/(1+(gamma-1)/2*x[1]^2) )^(gamma/(gamma-1)),
x[5]/x[4] - 1 - (2*gamma/(gamma+1)*(x[1]^2-1)),
x[2]^2 - ( (1+(gamma-1)/2*x[1]^2)/( gamma*x[1]^2-(gamma-1)/2 ) ),
x[5]/pe - ( (1+(gamma-1)/2*x[6]^2)/(1+(gamma-1)/2*x[2]^2) )^(gamma/(gamma-1)),
x[3]/Ae - ( (1+(gamma-1)/2*x[2]^2)/(1+(gamma-1)/2*x[6]^2) )^((gamma+1)/(2*(gamma-1)))*x[6]/x[2] ]
end
g(x) = f(x,gamma,As,Ae,p0,pe,ps)
# M1 M2 A1 p1 p2 Me
sol = nlsolve(g, [ 2, 0.5, 7e-4, 3e3, 6e3, 0.3 ], show_trace=true)
result = sol.zero
println("M1 = ",result[1])
println("M2 = ",result[2])
println("A1 = ",result[3])
println("p1 = ",result[4])
println("p2 = ",result[5])
println("Me = ",result[6])
流体特性値を計算するjuliaプログラムを以下に示す。
-衝撃波の上流側
profile-upstream.jl
using NLsolve
gasconst = 8.31446262
MW = 29e-3
cp = 1005
R = gasconst / MW
gamma = cp/(cp-R)
A0 = 0.0016129
p0 = 6895
T0 = 100
rho0 = p0/(R*T0)
a0 = sqrt( gamma*R*T0 )
Ts = 2/(gamma+1) * T0
ps = (2/(gamma+1))^(gamma/(gamma-1)) * p0
rhos = (2/(gamma+1))^(1/(gamma-1)) * rho0
as = sqrt( 2/(gamma+1) ) * a0
us = as
As = 0.00064516
n = 100
for i in 0:76
if i == 76 # <-- shock position
area = 0.0008129027418810686
x = 0.1921177321640935
else
x = i*(0.254/n)
xin = i*(10/n)
area = ( x<=0.127 ? (1.75-0.75*cos((0.2*xin-1)*pi) )*2.54e-2^2
: (1.25-0.25*cos((0.2*xin-1)*pi) )*2.54e-2^2 )
end
f(x,area,gamma,As) = ( area/As
- ( (1+(gamma-1)/2*x[1]^2)/((gamma+1)/2) )^((gamma+1)/(2*(gamma-1)))/x[1] )
g(x) = f(x,area,gamma,As)
ini = x<=0.127 ? 0.3 : 1.1
sol = nlsolve( g, [ ini ], show_trace=false)
result = sol.zero
M = result[1]
rho = ( ((gamma+1)/2)/(1+(gamma-1)/2*M^2) )^(1/(gamma-1)) * rhos
T = ((gamma+1)/2)/(1+(gamma-1)/2*M^2) * Ts
p = ( ((gamma+1)/2)/(1+(gamma-1)/2*M^2) )^(gamma/(gamma-1)) * ps
u = M*( ((gamma+1)/2)/(1+(gamma-1)/2*M^2) )^0.5 * us
println("$x $area $M $p $T $rho $u")
end
-衝撃波の下流側
profile-downstream.jl
using NLsolve
gasconst = 8.31446262 # m2 kg s-2 K-1 mol-1
MW = 29e-3
cp = 1005
R = gasconst / MW
gamma = cp/(cp-R)
A0 = 0.0016129
p0 = 6895
T0 = 100
rho0 = p0/(R*T0)
a0 = sqrt( gamma*R*T0 )
Ts = 2/(gamma+1) * T0
ps = (2/(gamma+1))^(gamma/(gamma-1)) * p0
rhos = (2/(gamma+1))^(1/(gamma-1)) * rho0
as = sqrt( 2/(gamma+1) ) * a0
us = as
As = 0.00064516
mdot = rhos*us*As
Ae = 0.00096774
pe = 5171
Me = 0.50200727578001
ae = gamma*pe*Me*Ae/mdot
Te = ae^2/(gamma*R)
ue = Me * ae
rhoe = pe/(R*Te)
n = 100
for i in 75:n
if i == 75 # <-- shock position
area = 0.0008129027418810686
x = 0.1921177321640935
else
x = i*(0.254/n)
xin = i*(10/n)
area = ( x<=0.127 ? (1.75-0.75*cos((0.2*xin-1)*pi) )*2.54e-2^2
: (1.25-0.25*cos((0.2*xin-1)*pi) )*2.54e-2^2 )
end
f(x,area,gamma,Me,Ae) = ( area/Ae
- ( (1+(gamma-1)/2*x[1]^2)/(1+(gamma-1)/2*Me^2) )^((gamma+1)/(2*(gamma-1)))*Me/x[1] )
g(x) = f(x,area,gamma,Me,Ae)
sol = nlsolve( g, [ 0.3 ], show_trace=false, ftol=1e-15)
result = sol.zero
M = result[1]
rho = ( (1+(gamma-1)/2*M^2)/(1+(gamma-1)/2*Me^2) )^(-1/(gamma-1)) * rhoe
T = ( (1+(gamma-1)/2*M^2)/(1+(gamma-1)/2*Me^2) )^(-1) * Te
p = ( (1+(gamma-1)/2*M^2)/(1+(gamma-1)/2*Me^2) )^(-gamma/(gamma-1)) * pe
u = M/Me*( (1+(gamma-1)/2*Me^2)/(1+(gamma-1)/2*M^2) )^(0.5) * ue
println("$x $area $M $p $T $rho $u")
end
case-C
流体特性値を計算するjuliaプログラムを以下に示す。
profile.jl
using NLsolve
gasconst = 8.31446262
MW = 29e-3
cp = 1005
R = gasconst / MW
gamma = cp/(cp-R)
A0 = 0.0016129
p0 = 6895
T0 = 100
rho0 = p0/(R*T0)
a0 = sqrt( gamma*R*T0 )
Ts = 2/(gamma+1) * T0
ps = (2/(gamma+1))^(gamma/(gamma-1)) * p0
rhos = (2/(gamma+1))^(1/(gamma-1)) * rho0
as = sqrt( 2/(gamma+1) ) * a0
us = as
As = 0.00064516
n = 100
for i in 0:n
x = i*(0.254/n)
xin = i*(10/n)
area = ( x<=0.127 ? (1.75-0.75*cos((0.2*xin-1)*pi) )*2.54e-2^2
: (1.25-0.25*cos((0.2*xin-1)*pi) )*2.54e-2^2 )
f(x,area,gamma,As) = ( area/As
- ( (1+(gamma-1)/2*x[1]^2)/((gamma+1)/2) )^((gamma+1)/(2*(gamma-1)))/x[1] )
g(x) = f(x,area,gamma,As)
ini = x<=0.127 ? 0.3 : 1.1
sol = nlsolve( g, [ ini ], show_trace=false, ftol=1e-15)
result = sol.zero
M = result[1]
rho = ( ((gamma+1)/2)/(1+(gamma-1)/2*M^2) )^(1/(gamma-1)) * rhos
T = ((gamma+1)/2)/(1+(gamma-1)/2*M^2) * Ts
p = ( ((gamma+1)/2)/(1+(gamma-1)/2*M^2) )^(gamma/(gamma-1)) * ps
u = M*( ((gamma+1)/2)/(1+(gamma-1)/2*M^2) )^0.5 * us
println("$x $area $M $p $T $rho $u")
end
OpenFOAMの計算に使用したデータはここからダウンロード可能である。
まとめ
ラバルノズル流れの理論解を導出し、OpenFOAMを使った解析結果と比較した。両者はよく一致していることが分かった。この資料が何らかの参考になれば幸いである。