Mathematicaでスキージャンプの飛距離を計算してみる

◼
  • 執筆者：技術部・降旗
    • ◼
  • 使用したバージョン：Mathematica14.3
    • ◼
  • 使用した主な関数：NDSolve

    • はじめに

    Mathematicaを使って微分方程式の解き方を説明するときに、高校物理で馴染みのある斜方投射の運動方程式を使うことがありますが、高校物理では空気抵抗は考慮せずに運動方程式を立てます。
    ですが、2月に幕を閉じたミラノ・コルティナ冬季オリンピックのスキージャンプを見ていて、空気抵抗が飛距離に影響を与える様子を見て、空気抵抗を考慮して運動方程式を立てるとどうなるかが気になったため、今回はスキーのジャンプの飛距離をMathematicaで計算できないか試してみました。

    STEP 1：空気抵抗を考慮せずに飛ばしてみる

    まずはできる限り条件をシンプルにするため、空気抵抗が存在しない真空中のジャンプを考えます。このとき、ジャンパーに働く力は重力だけです。ニュートンの第2法則に基づき、水平方向 (x) と垂直方向 (y) の動きを考えます。水平方向は、力が働かないため、
    m
    2
    d
    x
    2
    dt
    =0
    、垂直方向は重力
    mg
    が下向きに働くため、
    m
    2
    d
    y
    2
    dt
    =-mg
    となります。これより​
    x''(t)=0​​y''(t)=-g
    ​の式が求められ、空気抵抗を考えないと一般的な投射の運動方程式であり、体重が飛距離に関係しないことが分かります。
    続いてランディングバーン(着地斜面)について考えます。​
    ​下記のページによると、ジャンプ台の一つである大倉山ジャンプ競技場ではランディングバーンの斜度は最大37度であるとのことです。​
    ​https://www.age.ne.jp/x/sas/sas2005/Jump/Info_JP/JP_Setumei_2.htm​
    ​こちらもシンプルな条件にするため、今回は35度の直線であると設定します。​
    ​また、カンテ(踏み切り)の高さは3.3mであるそうですので、ジャンパーが飛び出す位置を原点(0,0)とし、(0,-3.3)の位置から35度の斜面があると設定します。
    In[]:=
    (*定数の定義*)​​height=3.3;(*カンテの高さm*)​​phi=35Degree;(*斜面の角度*)​​(*斜面の定義*)​​slope[x_]:=-height-x*Tan[phi]​​(*斜面を表示*)​​Plot[slope[x],{x,0,100},Filling->Bottom,FillingStyle->GrayLevel[0.9],PlotStyle->Black]
    Out[]=
    最後にジャンパーの初速について考えます。ジャンパーの飛び出す時の速さは約90km/hとのことですので、初速
    v
    0
    =25m/s
    とします。そして、先ほどのページによるとカンテは斜度が-11度とのことですので、初速のx成分、y成分をそれぞれ​
    x'(0)=
    v
    0
    Cos(-11°)​​y'(0)=
    v
    0
    Sin(-11°)
    ​とします。これらを踏まえて計算を行ってみます。
    In[]:=
    (*定数の定義*)​​g=9.81;(*重力加速度m/s^2*)​​theta=-11Degree;(*カンテの角度（下向き）*)​​v0=25;(*初速m/s*)​​
    In[]:=
    Module{sol,tEnd,landingX,landingY,distance},​​(*運動方程式を解く*)​​sol=NDSolve​​x''[t]==0,(*水平方向は等速*)​​y''[t]==-g,(*垂直方向は重力加速度*)​​x[0]==0,y[0]==0,(*初期位置は原点*)​​x'[0]==v0*Cos[theta],(*初速のx成分*)​​y'[0]==v0*Sin[theta],(*初速のy成分*)​​WhenEvent[y[t]<=slope[x[t]],​​"StopIntegration"](*y[t]が斜面以下になったら計算を終了*)​​,{x,y},{t,0,40}[[1]];​​tEnd=(x/.sol)["Domain"][[1,2]];(*計算終了時の時間を抽出*)​​(*着地位置と飛距離の計算*)​​landingX=x[tEnd]/.sol;​​landingY=y[tEnd]/.sol;​​distance=Norm[{landingX,landingY+height}];​​(*結果表示*)​​Show​​(*微分方程式の解をプロット*)​​ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,tEnd},PlotStyle->{Thick,Red}],​​(*斜面を表示*)​​Plot[slope[x],{x,0,300},Filling->Bottom,FillingStyle->GrayLevel[0.9],PlotStyle->Black],​​(*着地点を表示*)​​Graphics[{Red,PointSize[0.02],Point[{landingX,landingY}]}],​​Plot[slope[x],{x,0,landingX},PlotStyle->{Thick,Blue}]​​,AxesLabel->"水平距離 (m)","高さ (m)",PlotRange->{{0,Automatic},{-200,0}},​​PlotLabel->StyleStringForm"飛距離: `` m",Round[distance,0.1],14,Bold,Red​​
    Out[]=
    この条件では飛距離が約83.1mになることが分かりました。

    STEP 2：抗力と揚力を追加する

    続いて、ジャンプの飛距離に影響を与える抗力と揚力をモデルに追加してみます。抗力は進行方向と真逆の向きに働くブレーキの力で、次の式で表されるそうです。​​
    F
    D
    =
    1
    2
    C
    D
    ρS
    2
    v
    ​​・
    C
    D
    (抗力係数)：物体の形による空気の通しにくさ。・
    ρ
    (空気密度)：空気の濃さ。・
    S
    (投影面積)：正面から見た時の面積。・
    v
    (速度)：物体のスピード。​この式から抗力は速度の2乗に比例するため、スピードが出るほどブレーキがかかることが分かります。
    揚力は進行方向に対して垂直に働く力で、次の式で表されるそうです。​​
    F
    L
    =
    1
    2
    C
    L
    ρS
    2
    v
    ​​・
    C
    L
    (揚力係数)：スキー板や体の角度によって決まる浮き上がりやすさの係数。・
    ρ
    (空気密度)：空気の濃さ。・
    S
    (投影面積)：正面から見た時の面積。・
    v
    (速度)：物体のスピード。​どちらも投影面積が式に含まれており、揚力を大きくしようと面積を大きくすると抗力も大きくなってしまうということが式から分かります。
    これらの式を運動方程式に加えたいと思いますが、
    v
    は進行方向の速度であり、かつ抗力は速度と逆向き、揚力は速度と90度向きを変えた方向に働くため、これをxとy成分に分解する必要があります。それぞれの成分に分解して整理すると以下の加速度の式が得られます。​​
    x''(t)=
    1
    m
    -
    F
    D
    v
    x
    v
    -
    F
    L
    v
    y
    v
    ​​y''(t)=-g+
    1
    m
    -
    F
    D
    v
    y
    v
    +
    F
    L
    v
    x
    v
    ​ここで
    m
    が式に出てきましたので、体重により飛距離が変わるようになりました。この式を使って再度微分方程式を解いてみたいと思います。なお、定数についてはAIの回答を参考にして、投影面積は
    s=0.5
    2
    m
    、抗力係数、揚力係数はそれぞれ
    cl=1.0
    、
    cd=0.4
    とします。
    In[]:=
    (*追加の定数の定義*)​​m=70;(*体重+装備kg*)​​rho=1.225;(*空気密度kg/m^3*)​​s=0.5;(*投影面積m^2*)​​cl=1.0;(*揚力係数*)​​cd=0.4;(*抗力係数*)​​
    先ほどよりも飛距離が大きく伸びることが確認できました。
    しかし、参考となる記録がないため、残念ながらこれが現実的な値であるか分かりません。

    STEP 3：パラメーターの影響度を探る

    最後のステップでは、体重や投影面積、抗力係数、揚力係数などが飛距離にどう影響するかを分析します。
    まずは、パラメーターの変化に伴う結果をすぐに取得できるようにするために、それぞれの変数を引数に取って飛距離を求める関数と結果を表示する関数を作成します。
    体重の変化により、飛距離がどの程度変わるのか計算してみます。
    また、今までは飛距離の計算を(0,-3.3)から着地点までの直線距離として求めていましたが、これを上記の修正に伴い変更します。
    新しい斜面の設定を使った結果を表示してみます。
    今度は地面に着地するようになりました。
    この設定を使って再度体重の変化に伴う飛距離を計算してみます。
    体重が軽いほど飛距離が長くなるという様子が分かるようになりました。
    続いて揚力と抗力の変化に伴う飛距離の変化についても可視化してみます。
    着地点が268mを超えるか超えないかの2つの領域に分かれるのではないかと想像していましたが、大きく4つの領域に分かれているようです。
    抗力係数を固定し、揚力係数のみを変化させてみます。
    抗力係数を0.1としたときには、0.912、1.107、1.295を境に飛距離が大きく変わるようです。
    何が起きているか確認するために、それぞれの領域の結果を描画してみます。
    一つ目と二つ目の領域については想像通りで、斜面に着地するか、水平な箇所に着地するかの違いのようでした。
    そして面白いことに、揚力がさらに増えると紙飛行機のように下降と上昇を繰り返しながら徐々に落ちるような動きになるようです。
    最後に、ここまで数式や等高線マップで解析してきた内容を、今度は直感的に体感できるようにしてみます。
    以下のコードを実行すると、スライダーで各パラメーターを操作し、ジャンプの軌道がどう変化するかをリアルタイムにシミュレートできます。

    終わりに

    今回はスキーのジャンプの飛距離を計算できないか試してみました。
    実際にはジャンプ中にジャンパーの姿勢が変わることにより、投影面積や揚力/抗力係数なども変わってきますので、まだ現実に即した値からは遠いと考えられますが、空気抵抗を考慮した物体の投射運動としては思っていたよりもしっかりと計算できたのではないかと思います。
    今後の課題としては、ランディングバーンをより現実的な形状にすることにより実際の飛距離との比較を行ったり、追い風や向かい風の影響を加えたりするなどの変更が考えられます。