はじめに
なかなかゴルフが上手くならないのでゴルフボールの周りの空気の流れを計算して最も飛距離を伸ばすための打ち出し角度と回転速度を求めてみました。
- バウンドの計算は着地点の角度や芝の状態などにより大きく変化すると考えられるため、着地するまでのキャリー距離を飛距離として計算していきます。
飛距離を追求してもゴルフの上達に繋がらない可能性があります。
2019.5.6 17:40修正 コメント頂いた回転行列の部分を修正しました。これにより最も飛距離が伸びるのは低角度でバックスピンという結果になりました。
方針
まず弾道を計算するにあたり以下の条件を固定します。
- ボールの初速は60m/s (ヘッドスピード45m/sで真芯に当たった場合)
- ボールの進行方向に対して左右方向への移動が無い
- ボールの回転方向もトップスピン/バックスピンの方向に固定
この状況で飛距離を計算するために以下のようにボールの弾道を計算する必要があります。
弾道計算を行う中では、空気抵抗とそれにより影響されるボールの速度や回転速度を同時に計算しつつ、その時のボールの位置を計算していくことになります。 この計算を行うには以下の方法が考えられます。
- 空気抵抗と弾道を同時に計算：
- 非定常ソルバーの中で、各時刻でボールにかかる力を求める↔境界条件を更新を繰り返す中で飛距離を計算するソルバーを作成する
- このソルバーによる飛距離計算を複数パターンで行い最も飛距離の長い条件を探索する。
- 空気抵抗と弾道を別々に計算：
- 定常ソルバーを用いて複数条件で定常状態のボールにかかる力を求める
- その結果から条件と力の関係を近似する
- その近似式による飛距離計算を複数パターンで行い最も飛距離の長い条件を探索する。
最終的に今回は2を採用しました。
今回のケースでは1の条件はあまりに計算に時間がかかり過ぎてしまい、PCのスペックの関係で断念しました。
今回のように速度に対して対象となるメッシュが小さい場合にはどうしても時間解像度が細かくなり計算負荷が大きくなってしまいます。
ただ本格的に計算する場合には非定常でのみ計算できる状況（カルマン渦など）の影響を無視することができないため、相応の環境を用意して1を計算する必要があるかと思います。
というわけで、具体的には
- OpenFOAMのsimpleFoamで回転するゴルフボール周りの定常流体解析
- そこで得られた抗力・揚力・モーメントを用いてPythonで弾道と飛距離を数値計算
という流れで飛距離を計算していきます。
計算環境
- OS : Ubuntu 18.04 LTS
- メッシャ：SALOME 8.4.0
- ソルバー：simpleFoam (OpenFOAM v1812)
- Python：3.6.6 (Anaconda 4.5.2)
OpenFOAMによる空気抵抗の計算
まずはOpenFOAMを用いてゴルフボール周りの流れを計算していきます。乱流にはk-ω SSTモデルを使用し、回転はMRFで計算しています。
メッシュ
ゴルフボールの規格を参考に以下のようなゴルフボールを想定しました。
- 直径：42.67 mm
- 重量：45.93 g
- ディンプル
- 数：362個
- 直径：1.8 mm
- 深さ：0.6 mm
SALOMEを用いて下図のようなメッシュを作成しました。
全てヘキサで作成した後にrefineMeshにて一部をポリヘドラに細分化しています。また水色と黄色のエリアは別々に作成しており、その境界はcyclicAMIにより接続されています。
$ checkMesh
Mesh stats
points: 2428023
faces: 6752074
cells: 2162148
Overall number of cells of each type:
hexahedra: 2108612
polyhedra: 53536
Checking geometry...
Overall domain bounding box (-0.2 -0.1 -0.1) (0.3 0.1 0.1)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-1.07072e-16 7.16432e-18 -1.11005e-15) OK.
Max cell openness = 3.42628e-16 OK.
Max aspect ratio = 4.95719 OK.
Minimum face area = 1.84412e-08. Maximum face area = 3.3335e-05. Face area magnitudes OK.
Min volume = 7.28758e-12. Max volume = 1.66675e-07. Total volume = 0.0199604. Cell volumes OK.
Mesh non-orthogonality Max: 55.035 average: 14.9621
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.987239 OK.
Coupled point location match (average 0) OK.
基本となるケースの作成
OpenFOAMのソルバーであるsimpleFoam用のケースディレクトリを作成します。
$FOAM_TUTORIALS/incompressible/simpleFoam/motorBikeとほぼ一緒の条件でメッシュとcontrolDictと緩和係数を変更しました。
ここでcontrolDict内ではfunctionObjectを使ってボールにかかる力を求めています。
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v1812 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application simpleFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 20000; deltaT 1; writeControl timeStep; writeInterval 1000; purgeWrite 1; writeFormat binary; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; functions { force1 { type forces; libs ( "libforces.so" ); writeControl writeTime; log yes; patches ( wall_ball ); rho rhoInf; rhoInf 1; CofR (0 0 0); } } // ************************************************************************* //
また以下のように強めな緩和係数を設定しました。
SIMPLE
{
nNonOrthogonalCorrectors 1;
consistent yes;
}
potentialFlow
{
nNonOrthogonalCorrectors 10;
}
relaxationFactors
{
fields
{
p 0.4;
}
equations
{
U 0.4;
k 0.3;
omega 0.3;
}
}
条件毎にソルバー実行
今回は速度と回転軸を固定し、複数の回転速度で計算を行います。そのために以下のようにPythonにてケーススタディを行います。
まず以下のような関数を作成します。
import os, shutil from datetime import datetime import threading, time, glob, subprocess, requests import pandas as pd import paraview.simple as pv from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile from math import pi, sin, cos def makeGolfCase(resultDir, orgCase, **kwargs): if 'newName' in kwargs.keys(): baseName = kwargs['newName'] else: baseName = os.path.basename(orgCase) n = 0 newCase = os.path.join(resultDir, '{}_{}'.format(baseName, n)) while os.path.exists(newCase): newCase = os.path.join(resultDir, '{}_{}'.format(baseName, n)) n += 1 print('Created new case named', newCase) name = os.path.basename(newCase) shutil.copytree(orgCase,newCase) if not os.path.exists(newCase+'/0'): if os.path.exists(newCase+'/0.orig'): shutil.copytree(newCase+'/0.orig', newCase+'/0') if 'U' in kwargs: u = float(kwargs['U']) with ParsedParameterFile(newCase+'/0/U') as UFile: if kwargs['isInitialCase']: UFile.content['internalField'] = 'uniform ({:.2f} 0 0)'.format(u) UFile.content['boundaryField']['inlet']['value'] = 'uniform ({:.2f} 0 0)'.format(u) UFile.writeFile() with ParsedParameterFile(newCase+'/system/controlDict') as controlDict: controlDict.content['functions']['force1']['magUInf'] = '{:.2f}'.format(u) controlDict.writeFile() if 'rpm' in kwargs: if kwargs['rpm'] > 0: omega = kwargs['rpm'] * 2 * pi / 60.0 with ParsedParameterFile(newCase+'/constant/MRFProperties') as MRF: MRF.content['MRF1']['active'] = 'true' MRF.content['MRF1']['omega'] = 'constant {:.6f}'.format(omega) MRF.content['MRF1']['axis'] = '({:.6f} {:.6f} {:.6f})'.format(kwargs['axis'][0],kwargs['axis'][1],kwargs['axis'][2]) MRF.writeFile() return newCase
この中では
- ケースのコピー（名前が重複しないように）
- 速度Uを設定
- MRF回転条件を設定
を行っています。この関数を用いて以下のように実行します。
masterDir = '../study' if not os.path.exists(masterDir): os.mkdir(masterDir) # このkwargsをmakeGolfCase関数に渡す kwargs = { 'U' : 5, 'rpm' : 100, 'axis' : [0,0,1], 'isInitialCase' : True } for U in [40,50,60]: for rpm in [60, 1000, 3000, 5000 ]: for angle in [0]: kwargs['U'] = U kwargs['rpm'] = rpm kwargs['axis'] = [sin(angle*pi/180), 0, cos(angle*pi/180)] newName = makeGolfCase(masterDir, '../makeMesh/merge/case', **kwargs) cwdDir = os.getcwd() os.chdir(newName) subprocess.check_output(['./Allrun']) os.chdir(cwdDir)
これにより速度を[40,50,60]m/s × 回転速度[60, 1000, 3000, 5000]rpmの12パターンで計算を行います。
これで以下のようなボールにかかる力の計算が12パターン行われました。
Pythonによる弾道の計算
それでは得られた結果をもとに弾道と飛距離を計算していきます。使用するのはOpenFOAMのfunctionObjectsで出力されるファイルです。これをpandasで読み込んで加工しながら使っていきます。
結果の確認
計算が完了したらまずはちゃんと収束しているかどうかを確認してみましょう。 以下のようにfunctionObjectsにより生成されたforceとmomentをtimeStepを横軸にプロットしてみます。 今回は一度pandasに放り込むことで簡単に描画しています。
caseDir = '../study' cases = [] for d in glob.glob(caseDir+'/*'): if os.path.exists(d+'/system/controlDict'): cases.append(d) cases = sorted(cases) df = pd.DataFrame() fig = plt.figure(figsize=(20,10)) ind = 0 for case in cases: name = os.path.basename(case) forceFile = case+'/postProcessing/force1/0/force.dat' momdentFile = case+'/postProcessing/force1/0/moment.dat' forceDF = pd.read_csv(forceFile, sep='\t', header=3,index_col=0) force = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in forceDF[forceDF.columns[0]].items() ]) momentDF = pd.read_csv(momdentFile, sep='\t', header=3,index_col=0) moment = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in momentDF[momentDF.columns[0]].items() ]) for i in range(3): plt.subplot(2,3,i+1) plt.plot(range(1000),force[-1000:,i], label=name) plt.subplot(2,3,i+4) plt.plot(range(1000),moment[-1000:,i], label=name) i = 0 for axis in ['_X','_Y','_Z']: plt.subplot(2,3,i+1) plt.title('force'+axis) plt.subplot(2,3,i+4) plt.title('moment'+axis) if i == 2: plt.legend() i += 1 fig.savefig('fig/1.results.png', bbox_inches="tight")
force_Xとforce_Yとmoment_Zはいい感じに収束しています。 なお今回の計算条件ではZ軸が回転軸となっているので
- X方向：抗力
- Y方向：揚力
- Z軸モーメント：回転モーメント
となります。
そのため今回の弾道計算の中では進行方向に対して左右にかかる力(force_Z)は無視します。
またモーメントについても回転を弱める方向のモーメント(moment_Z)以外については無視することとします。
というわけでよく収束されていたforce_Xとforce_Yとmoment_Zのみを用いて弾道を計算していきます。
各条件における計算結果を抽出
各条件によりかかる力/モーメントを同じくpandasのDataFrame形式にまとめていきます。 今回は定常計算の各timeStepで得られる数値の後半100step分の結果を平均した値を結果として使用します。
df = pd.DataFrame() ind = 0 for case in cases: name = os.path.basename(case) forceFile = case+'/postProcessing/forceCoeffs1/0/force.dat' momdentFile = case+'/postProcessing/forceCoeffs1/0/moment.dat' forceDF = pd.read_csv(forceFile, sep='\t', header=3,index_col=0) force = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in forceDF[forceDF.columns[0]].items() ]) momentDF = pd.read_csv(momdentFile, sep='\t', header=3,index_col=0) moment = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in momentDF[momentDF.columns[0]].items() ]) rpm = int(round(ParsedParameterFile(case+'/constant/MRFProperties').content['MRF1']['omega'][1] * 60 / 2 / pi, 0)) U = int(round(ParsedParameterFile(case+'/system/controlDict').content['functions']['forceCoeffs1']['magUInf'])) ser = {'rpm':rpm,'U':U} n = 0 for axis in ['_X','_Y','_Z']: ser['force'+axis] = force[-100:,n].mean() n += 1 n = 0 for axis in ['_X','_Y','_Z']: ser['moment'+axis] = moment[-100:,n].mean() n += 1 df[ind] = pd.Series(ser) ind += 1 # 逆回転 ser = {'rpm':-rpm,'U':U} n = 0 for axis in ['_X','_Y','_Z']: if n != 1: ser['force'+axis] = force[-int(len(force)/10):,n].mean() else: ser['force'+axis] = -force[-int(len(force)/10):,n].mean() n += 1 n = 0 for axis in ['_X','_Y','_Z']: if n != 2: ser['moment'+axis] = moment[-int(len(moment)/10):,n].mean() else: ser['moment'+axis] = -moment[-int(len(moment)/10):,n].mean() n += 1 df[ind] = pd.Series(ser) ind += 1 df = df.T df = df.set_index(['rpm','U'], drop=True) df
定常計算の中では重力加速度については計算されていないため、Y方向を反転させた結果をそのまま逆回転の結果として使用しています。 あえてこの結果をDataFrameに入れた理由は、後に出てくる近似式の無回転状態の精度を高めるためです。
データの傾向を確認
以上のようにまとめた結果を以下のように図示してみます。
fig = plt.figure(figsize=(20,10)) ind = 0 for name in ['force_X', 'force_Y', 'moment_Z']: ans = {'rpm':[],'U':[]} for key,val in df[name].items(): rpm = key[0] U = key[1] rpmKey = '{}rpm'.format(int(rpm)) UKey = '{}m/s'.format(int(U)) if not rpmKey in ans: ans[rpmKey] = [] ans['rpm'].append(rpm) if not UKey in ans: ans[UKey] = [] ans['U'].append(U) ans[rpmKey].append(val) ans[UKey].append(val) plt.subplot(2,3,ind+1) for rpm in ans['rpm']: rpmKey = '{}rpm'.format(int(rpm)) plt.scatter(ans['U'],ans[rpmKey],label=rpmKey) plt.title(name) plt.legend() plt.subplot(2,3,ind+4) for U in ans['U']: UKey = '{}m/s'.format(int(U)) plt.scatter(ans['rpm'],ans[UKey],label=UKey) plt.title(name) plt.legend() ind += 1 fig.savefig('fig/2.plot.png', bbox_inches="tight")
なんとなく綺麗な規則性がありそうですね。 少々雑ですが以下のような式で近似してみましょう。
- U : 速度 [m/s]
- rot : 回転速度 [rpm]
- a,b,c,d,e,f : 定数
近似にはscipyのcurve_fitを用います。ます以下のような関数を用意します。
from scipy.optimize import curve_fit def func(xy, a,b,c,d,e,f): x = xy[0] y = xy[1] return (a*x**2 + b*x + c)*(d*y**2+e*y+f)
この上で以下のように実行することでforce_X、force_Y、moment_Zそれぞれについてフィッティングが行われます。
fig = plt.figure(figsize=(20,20)) rpmX = list(range(-5000,5001,100)) UX = list(range(0,61,5)) betas = {} ind = 0 for name in ['force_X', 'force_Y', 'moment_Z']: ans = {'rpm':[0.0],'U':[0.0]} for key,val in df[name].items(): rpm = key[0] U = key[1] rpmKey = '{}rpm'.format(int(rpm)) UKey = '{}m/s'.format(int(U)) if not rpmKey in ans: ans[rpmKey] = [0.0] ans['rpm'].append(rpm) if not UKey in ans: ans[UKey] = [0.0] ans['U'].append(U) ans[rpmKey].append(val) ans[UKey].append(val) ans['0rpm'] = [v for v in ans[rpmKey]] ans['0m/s'] = [v for v in ans[UKey]] beta, pcov = curve_fit(func, np.array(list(df[name].index)).T, df[name].values) betas[name] = beta plt.subplot(3,3,ind+1) for rpm in ans['rpm']: rpmKey = '{}rpm'.format(int(rpm)) plt.plot(UX, [func((rpm,u),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) for u in UX], linestyle='dashed',label='fitting curve') plt.scatter(ans['U'],ans[rpmKey],label=rpmKey) plt.title(name+' vs velocity') plt.xlabel('velocity [m/s]') plt.legend() plt.subplot(3,3,ind+4) for U in ans['U']: UKey = '{}m/s'.format(int(U)) plt.plot(rpmX, [func((r,U),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) for r in rpmX], linestyle='dashed',label='fitting curve') plt.scatter(ans['rpm'],ans[UKey],label=UKey) plt.title(name+' vs rotation') plt.xlabel('rotation [rpm]') plt.legend() rpms, Us = np.mgrid[-5000:5000:100, 0:61:5] Z = func((rpms, Us),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) plt.subplot(3,3,ind+7) plt.contour(rpms, Us, Z, linewidths=0.5, colors='k') plt.contourf(rpms, Us, Z, cmap=plt.cm.Spectral) plt.title(name+' contour') plt.xlabel('rotation [rpm]') plt.ylabel('velocity [m/s]') ind += 1 fig.savefig('fig/3.curvefit.png', bbox_inches="tight")
この中の下のコンター図は近似式を用いて計算した値をプロットしたものです。 OpenFOAMで計算していない条件についてもいい感じに補完されているように見えます。
なお得られた近似式の各定数は以下のようになりました。
tmpdf = pd.DataFrame(betas).T tmpdf.columns = ['a', 'b', 'c', 'd', 'e', 'f'] tmpdf
以上のようにOpenFOAMの結果に則った近似式がそれぞれの力についてできました。
近似式を用いて弾道を計算する
速度と回転速度を与えたらボールに加わる力が得られるようになりましたのでこれを使って弾道を計算していきます。
弾道計算にはオイラー法っぽい感じで時刻を細かく刻んで進行していき
前の時刻の結果から次の時刻の条件を計算するといったことを繰り返していきます。
2019.5.6 17:40修正 コメントをいただき回転行列の箇所を修正しました。
r = 42.67e-3/2 # m m = 45.93e-3 # kg def newStatus(rpms, vecU, dt, pos): pos.append([pos[-1][0]+vecU[0]*dt, pos[-1][1]+vecU[1]*dt]) rpm = rpms[-1] magU = (vecU[0]**2 + vecU[1]**2)**0.5 cosT = vecU[0] / magU sinT = vecU[1] / magU f = {} for key,beta in betas.items(): f[key] = func((rpm, magU),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) force = [-f['force_X']*cosT - f['force_Y']*sinT, -f['force_X']*sinT + f['force_Y']*cosT] # 修正しました。 #force = [-f['force_X']*sinT - f['force_Y']*cosT, -f['force_X']*cosT + f['force_Y']*sinT] vecU = [vecU[0] + force[0]/m*dt, vecU[1] + (force[1]/m - 9.81)*dt] rpms.append(rpm + f['moment_Z']/((2/5)*m*r**2)*dt*60/2/pi) #print(vecU, pos[-1]) if pos[-1][1] > 0 and len(pos) < 1000: newStatus(rpms, vecU, dt, pos) def distance(magU, rpm, angle): pos = [[0.0, 0.0]] vecU = [magU*cos(pi*angle/180),magU*sin(pi*angle/180)] rpms =[rpm] dt = 1e-3 newStatus(rpms, vecU, dt, pos) return np.array(pos)[:,0].max() def distancePlot(magU, rpm, angle): pos = [[0.0, 0.0]] vecU = [magU*cos(pi*angle/180),magU*sin(pi*angle/180)] rpms =[rpm] dt = 1e-3 topback = 'top' if rpm > 0: topback = 'back' label = '{}m/s, {}deg, {}rpm {} spin'.format(int(magU), int(angle), int(abs(rpm)), topback) newStatus(rpms, vecU, dt, pos) p = np.array(pos) plt.subplot(2,1,1) plt.plot(p[:,0],p[:,1],label=label) plt.ylabel('height [m]');plt.xlim(0,200);plt.ylim(0,100) plt.subplot(2,1,2) plt.plot(p[:,0],rpms,label=label) plt.xlabel('length [m]');plt.ylabel('rotating speed [rpm]'); plt.xlim(0,200);plt.ylim(-5000,5000)
説明の前に実行してみましょう。以下のように速度一定=60m/s、仰角一定=45°、回転速度を[-5000,-2500,0,2500,5000]の5パターンで計算を行うと各弾道を得ることができます。
fig = plt.figure(figsize=(15,7)) for rpm in [-5000,-2500,0,2500,5000]: distancePlot(60, rpm, 45) plt.legend() fig.savefig('fig/4.lines.png', bbox_inches="tight")
2019.5.6 17:40修正 図を差し替えました。
この上の図が弾道です。横軸に距離、縦軸に高さを取っています。また下の図では回転速度が減衰していっている様子が図示されています。
newStatus関数は再帰関数になっており、posというリストに各時刻の位置がappendされていきます。 この各時刻において速度と回転速度から力が計算され、そこからさっきの近似式を用いて次の時刻の回転速度と速度ベクトルが求められます。
同じような内容のdistance関数では飛距離を返すようになっています。
複数条件で計算
以上のようにdistance関数で各条件で飛距離を求めることができるようになりましたので、多くのパターンで計算を行ってみましょう。 なお速度は60m/sに固定した条件で計算を行います。 また最終的に2変数のコンター図にしたいのでX(回転速度)、Y(仰角)、Z(飛距離)で結果をまとめられるように計算を行います。
X, Y = np.mgrid[0:8001:200, 6:61:2] Z = np.zeros_like(X) maxDist = 0.0 maxCond = [0,0] for i in range(len(X)): for j in range(len(X[i])): val = distance(60, X[i][j], Y[i][j]) Z[i][j] = val if val > maxDist: maxDist = val maxCond = [X[i][j], Y[i][j]] print('Maximum distance is {:.1f} m when {:.1f} rpm and {:.1f} degree'.format(maxDist, maxCond[0], maxCond[1])) fig = plt.figure(figsize=(15,10)) plt.contour(X, Y, Z, linewidths=0.5, colors='k') plt.pcolor(X,Y, Z, cmap=plt.cm.Spectral) plt.xlabel('rotation [rpm] (+ : back spin, - : top spin)') plt.ylabel('shot angle [degree]') plt.colorbar() fig.savefig('fig/5.pcolor.png', bbox_inches="tight")
2019.5.6 17:40修正 図を差し替えました。
Maximum distance is 180.2 m when 6000.0 rpm and 18.0 degree
という訳で60m/sで打った時に最も飛距離がでるのは仰角18°の方向に6000rpmのバックスピンであることがわかりました。
最も飛ぶ弾道を図示
最後にこの条件を図示してみましょう。
fig = plt.figure(figsize=(15,7)) distancePlot(60, maxCond[0], maxCond[1]) plt.legend() fig.savefig('fig/6.mostFar.png', bbox_inches="tight")
2019.5.6 17:40修正 図を差し替えました。
低めに高回転でホップさせる感じがいいのでしょうか。
最後に
いかがでしたか？
みなさんもドラコンなどの際には仰角18°の方向に6000rpmのバックスピンをかけて打ってみてください。
