水ロケット2号機の設計 その2 -水ロケットの推力を求める-

だいぶ前になりますが、ペットボトルの形状を計測して数式化する準備をしました。
今回は推力を計算したいと思います。

基本的には以下参考URLの1番をなぞったうえで、ペットボトルの形状を反映させた式にして解いてみます。2つ目のリンクでは大口径ノズルで慣性も考慮した噴出の運動方程式を解かれています。今回は使わないけど参考まで。

スポンサーリンク
レクタングル(大)広告

参考URL

  1. 笹川民雄 「水ロケットと空気ロケット」(PDF)
    http://www.mars.dti.ne.jp/~stamio/rocket.pdf
    (最終閲覧日2017年5月4日)
  2. まいすた ヘタに考え休みにニタリ
    「ロケットを打ち上げろ」
    http://primzahl.seesaa.net/article/393289936.html
    「水ロケットの運動方程式」
    http://primzahl.seesaa.net/article/393289945.html
    (最終閲覧日 2017年5月4日)

基礎式の導出

変数定義

まずは変数の定義から。以下ペットボトルを横から見た図。

この絵だと右側に水を噴出して、その反作用でペットボトルは左に飛んでいきます。ペットボトル内の圧力をP、ボトル内の水の体積をV、水面の面積をS、ボトル内空気の体積をVair、温度をT、ノズル出口の面積をSt、ロケットから見た水の噴出速度をvとします。ロケット固定座標系としてボトル底面から口に向かってx軸を取って底面を0とします。またボトル全体の体積をVpet、水の密度をρ、大気圧をPA、空気の比熱比をγ(ガンマ)、添え字0は初期状態を表すものとします。

水の噴出速度

まずは水の噴出速度を求めます。水の流れを定常流としてベルヌーイの定理から1)、

$$ \frac{1}{2}0+\frac{P}{\rho} = \frac{1}{2}v^2+\frac{P_A}{\rho} \\
v=\sqrt{ \frac{2(P-P_A)}{\rho} } … ① $$

※水面の位置で流速を0としています。ここの流速をv1と置いて水の圧縮性を無視するとv1=v*St/S(V)となるので、上の式に代入して解いてもいいですが、0と置いても計算結果にほとんど差が無かったので0としてます。

水の体積

次に水の体積V。質量保存則より、

$$ \rho \frac{dV}{dt}=\rho S_t v … ② $$

ボトル内圧力

最後にボトル内圧力Pを求めます。水の噴出は短時間で行われるため断熱膨張とすると、

$$ P V_{air}^\gamma = P_0 V_{air0}^\gamma \\
\Leftrightarrow P=P_0(\frac{V_{air0}}{V_{air}})^\gamma … ③\\ $$

水の体積の変数分離型へ変形

式がそろったので変形してVの変数分離型へ変形します。
③を①へ代入して、

$$ v = \sqrt{ \frac{2}{\rho} \{P_0(\frac{V_{air0}}{V_{air}})^\gamma – P_A \} } $$

これを②に代入。またVair=Vpet-Vを代入。

$$ \frac{dV}{dt} = S_t \sqrt{ \frac{2}{\rho} \{P_0(\frac{V_{air0}}{V_{pet}-V})^\gamma – P_A \} } … ④ $$

というわけで、この④の式をルンゲクッタ法を使って数値的に解けば、水の体積が求まって、そこから圧力と水の噴出速度も求まります。

推力

ここまでで水の噴出速度が求まるので、推力を求めることができます。

推力をFとすると、Fは噴出した水の反作用なので、噴出した水の運動量変化に等しい。単位時間あたりに噴出する水の質量は噴出速度vからρ*St*vなので、

$$ F=\rho S_t v * v \Leftrightarrow F=\rho S_t v^2 $$

と、推力が求まります。

実際の計算コード

実際の計算はプログラムにやらせます。せっかくなのでなにか言語を覚えようと迷った挙句に科学技術計算に強いと評判のpythonにしてみました。

プログラミング自体は学生時代にCを少しかじった程度なのでトリッキーなことや最適化はやってません。というかできません。

以下、コードを掲載して今回は終わります。空気の噴出ステージは次の機会に追加予定。

# coding: utf_8
# water rocket thrust calcuration
# Date: 2017-05-12
# Version: 0.1
'''
   Copyright 2017 Yasubumi KANAZAWA (camelinsect@wings2fly.jp)

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
'''

import math
import sys

#----- 設定 -----
St = (7.7/2)**2*math.pi/1000.0**2 #ノズル出口断面積[m^2]
dt = 0.05 #時間ステップ[s]
m=0.134 # 機体質量 [kg]
P0 = 686000.0 # 初期タンク内圧[Pa]
V0 = 1.000/1000 # 初期タンク内水体積(0.5L) [m^3]
T0 = 300.0 #初期タンク内空気温度[K]

#----- 定数 -----
g=9.8 #重力加速度[m/s^2]
PA = 101300.0 #大気圧[Pa]
gamma = 1.4 #比熱比
rho = 1000.0 #水密度 [kg/m^3]

#----- ペットボトル関連 -----
Xpet = [65.0,187.0,200.0,240.0,280.0,306.0] # ペットボトル形状が変化する位置(ボトル底面からの距離)
Vp = [0,0,0,0,0] # ペットボトルの体積
Vp[0] = 0.0094393/1000        # Xpet[5]まで水を入れたとき…出口部
Vp[1] = Vp[0]+0.1379391/1000  # Xpet[4]まで水を入れたとき…2次関数部
Vp[2] = Vp[1]+0.2546366/1000  # Xpet[3]まで水を入れたとき…拡大部
Vp[3] = Vp[2]+0.0864325/1000  # Xpet[2]まで水を入れたとき…縮小部
Vp[4] = Vp[3]+0.7761305/1000  # Xpet[1]まで水を入れたとき…円筒部(ここまでを水を入れる許容量とする)
Vpet = Vp[4]+0.3602150/1000   # Xpet[0]まで水を入れたとき…ペットボトル全体の体積

#ペットボトル内の水量から水面位置を求める
def calc_x(V):
    dV=0.0
    dx=0.0
    x=0.0
    if V>Vpet:
        print "ペットボトル容量オーバー"
        return -1
    elif V>Vp[4]:
        print "水量制限オーバー"
        return -1
    elif V>Vp[3]:
        #print "円筒部"
        dV=V-Vp[3]
        dx=dV/(0.045**2*math.pi)
        x=0.187-dx
    elif V>Vp[2]:
        #print "拡大部"
        dV=V-Vp[2]
        dx=0.007
        while True:
            dx1 = dx -(0.024785741*dx**3-0.022716131*dx**2+0.006939778*dx-dV)/(0.074357223*dx**2-0.045432263*dx+0.006939778)
            if abs(dx1-dx) < 0.0000001: dx=dx1 break dx=dx1 x=0.2-dx elif V>Vp[1]:
        #print "縮小部"
        dV=V-Vp[1]
        dx=0.01
        while True:
            dx1 = dx -(0.010471976*dx**3+0.013508848*dx**2+0.005808805*dx-dV)/(0.031415927*dx**2+0.027017697*dx+0.005808805)
            if abs(dx1-dx) < 0.0000001: dx=dx1 break dx=dx1 x=0.24-dx elif V>Vp[0]:
        #print "2次関数部"
        dV=V-Vp[0]
        dV=0.1379391/1000 - dV
        dx=10
        while True:
            dx1 = dx -((5.667671E-05*dx**5+9.521310E-04*dx**4-4.787537E-01*dx**3-4.868640E+00*dx**2+1.852407E+03*dx)*math.pi*10**-9-dV)/((2.833836E-04*dx**4+3.808524E-03*dx**3-1.436261E+00*dx**2-9.737279E+00*dx+1.852407E+03)*math.pi*10**-9)
            if abs(dx1-dx) < 0.0000001: dx=dx1 break dx=dx1 dx = dx/1000 x=0.24+dx elif V>0:
        #print "出口部"
        dV=V
        dx=dV/(0.01075**2*math.pi)
        x=0.306-dx
    else:
        print "水体積がマイナス"
        return -1
    return x

#水面位置から水面の面積を求める
def calc_S(x):
    S=0.0
    r=0.0
    dx=0.0
    if x>0.306:
        print "水面位置がボトルの外"
        return -1
    elif x<0.0:
        print "水面位置がマイナス"
        return -1
    elif x<0.065:
        print "水容量オーバー"
        return -1
    elif x<0.187:
        #print "円筒部"
        r=0.045
    elif x<0.200:
        #print "拡大部"
        dx=x-0.187
        r=0.045+dx*2.0/13.0
    elif x<0.240:
        #print "縮小部"
        dx=x-0.200
        r=0.047-dx*1.0/10.0
    elif x<0.280: #print "2次関数部" dx=(x-0.240) r=(-0.016834*(x*1000.0)**2+7.9672*(x*1000.0)-899.45)/1000 elif x>=0.280:
        #print "出口部"
        r=0.01075
    S=r**2*math.pi
    return S


#----- 初期化 -----
x0 = calc_x(V0) #初期水面位置[m]
if x0 == -1:
    print "初期水面位置エラー"
    sys.exit(-1)
S0 = calc_S(x0) # タンク断面積[m^2]
if S0 == -1:
    print "初期水面面積エラー"
    sys.exit(-1)
Vair0 = Vpet-V0 # 初期タンク内空気体積[m^3]
t=0.0 # 時間
P=P0 # タンク内圧
V=V0 # タンク内水体積[m^3]
S=S0 # 水面位置のタンク断面積
Vair = Vair0 # タンク内空気体積[m^3]
T=T0 # タンク内空気温度[K]
v=0.0 # 噴出速度
x=x0 # 水面位置
F=0.0 # 推力
M=m+V0*rho # 全備質量

#----- 関数 -----
# Vの変数分離型 微分方程式
def fV(V):
    x=calc_x(V)
    S=calc_S(x)
    dVdt = St* math.sqrt(2/rho * (P0*(Vair0/(Vpet-V))**gamma -PA))
    return dVdt

# ルンゲクッタでVを求める
def calc_V(V):
    k=[0.0,0.0,0.0,0.0]
    k[0] = fV(V)
    k[1] = fV(V+dt*0.5*k[0])
    k[2] = fV(V+dt*0.5*k[1])
    k[3] = fV(V+dt*k[2])
    V1 = V-dt/6.0 *(k[0]+2.0*k[1]+2.0*k[2]+k[3])
    return V1


#----- メインプログラム -----
# 初期値確認
print "初期設定値"
print "時間ステップ = "+str(dt)+"[sec]"
print "ノズル断面積 = "+str(St*10**6)+"[mm^2]"
print "初期水量 = "+str(V0*1000)+"[L]"
print "機体質量 = "+str(m)+"[kg]"
print "初期タンク内圧 = "+str(P/100)+"[hPa]"
print "初期タンク内温度 = "+str(T0-273.15)+"[℃]"
print "\n"

# 水噴射ステージ
print "時間(t) 水面位置(x) 圧力(P) 噴出速度(v) 水面面積(S) ロケット質量(M) 推力(F) 水体積(V)"
#初期(t=0)
v=math.sqrt(2*(P-PA)/rho)
F=rho*St*v**2
print str(t),
print str(x),
print str(P),
print str(v),
print str(S),
print str(M),
print str(F),
print str(V)
t+=dt
while True:
    V1=calc_V(V)
    if V1<0: # ここで残りの水噴射の計算を実施 if V>Vp[0]:
            # 水面が出口部に達していなければエラーで終了
            print "時間当たりの水体積変化が大きすぎます。時間の刻みを小さくしてください"
            break
        else:
            # 残りは噴出速度vを一定として排出時間とタンク内圧力変化を計算。
            x=Xpet[5]/1000
            Vair=Vpet
            P=P0*(Vair0/Vair)**gamma
            if P-PA<0.1:
                print "タンク内圧が大気圧まで下がりました。初期圧力を上げてください"
                break
            S=calc_S(x)
            M=m
            v=math.sqrt(2*(P-PA)/rho)
            h=V/(v*St)
            t+=h
            F=rho*St*v**2
            V=0
            print str(t),
            print str(x),
            print str(P),
            print str(v),
            print str(S),
            print str(M),
            print str(F),
            print str(V)
            print "水噴出ステージ正常終了"
            break
    V=V1
    x=calc_x(V)
    Vair=Vpet-V
    P=P0*(Vair0/Vair)**gamma
    if P-PA<0.1:
        print "タンク内圧が途中で大気圧まで下がりました。初期圧を上げてください"
        break
    v=math.sqrt(2*(P-PA)/rho)
    S=calc_S(x)
    M=m+V*rho
    F=rho*St*v**2
    print str(t),
    print str(x),
    print str(P),
    print str(v),
    print str(S),
    print str(M),
    print str(F),
    print str(V)
    t+=dt
スポンサーリンク
レクタングル(大)広告

シェアする

  • このエントリーをはてなブックマークに追加

フォローする

スポンサーリンク
レクタングル(大)広告