DATCOMで安定微係数の計算を行うエクセルファイルを作ったので紹介する。
(2024.02.11 随時更新中)
↓計算に使用したエクセルファイル（マクロ付き）
使用しているマクロ
線形補間
Function Interpolation3D(ByVal x As Double, ByVal y As Double, ByVal z As Double, ParamArray ranges() As Variant) As Double
'https://excel-ubara.com/excelvba1/EXCELVBA433.html
Dim z1, z2 As Double
Dim z1_index, z2_index As Long
Dim nz As Long
Dim f1, f2, f As Double
'Get reference values of z and its index
nz = (UBound(ranges, 1) + 1) / 4 - 1
For i = 0 To nz - 1
If ranges(4 * i + 2) <= z And ranges(4 * (i + 1) + 2) >= z Then
z1_index = i
z2_index = i + 1
z1 = ranges(4 * z1_index + 2)
z2 = ranges(4 * z2_index + 2)
End If
Next i
'Constant value extrapolation in out of range of known values
If z <= ranges(4 * 0 + 2) Then z1_index = 0
If z >= ranges(4 * nz + 2) Then z1_index = nz - 1
'Get interpolated values against x and y
f1 = Interpolation2D(x, y, ranges(4 * z1_index + 0), ranges(4 * z1_index + 1), ranges(4 * z1_index + 3))
f2 = Interpolation2D(x, y, ranges(4 * z2_index + 0), ranges(4 * z2_index + 1), ranges(4 * z2_index + 3))
'Get interpolated values against z
If z <= ranges(2) Then
f = f1
ElseIf z >= ranges(4 * nz + 2) Then
f = f2
Else
f = f1 + (f2 - f1) / (z2 - z1) * (z - z1)
End If
Interpolation3D = f
End FunctionFunction Interpolation2D(ByVal x As Double, ByVal y As Double, ByRef x_range As Variant, ByRef y_range As Variant, ByRef f_range As Variant) As Double
'https://codepool.hatenablog.com/entry/2020/06/03/015921
'1:lower, 2:upper
Dim y1, y2 As Double 'reference value both sides of y
Dim y1_index, y2_index As Long 'reference value index both sides of y
Dim f1, f2, f As Double 'interpolated value
Dim tmp() As Variant
Dim i As Long
' Get reference value index of y
With WorksheetFunction
If y < .Min(y_range) Then y = .Min(y_range) 'Constant value extrapolation in range under lower limit
y1_index = .Match(y, y_range, 1)
y2_index = .Min(y1_index + 1, y_range.Count) 'Constant value extrapolation in range over the upper limit
End With 'if x or y out of known value range, 1 and 2 indexes become same value
' Get reference values of y
y1 = y_range(y1_index)
y2 = y_range(y2_index)
' Get reference values of f
f1 = Interpolation1D(x, x_range, f_range(y1_index))
f2 = Interpolation1D(x, x_range, f_range(y2_index))
'Get interpolated values against y
If y2 <> y1 Then
f = f1 + (f2 - f1) * (y - y1) / (y2 - y1)
Else
f = f1
End If
Interpolation2D = f
End FunctionFunction Interpolation1D(ByVal x As Double, ByRef x_range As Variant, ByRef f_range As Variant) As Double
'https://codepool.hatenablog.com/entry/2020/06/03/015921
'1:lower, 2:upper
Dim x1, x2 As Double 'reference value both sides of x
Dim x1_index, x2_index As Long 'reference value index both sides of x
Dim f, f1, f2 As Double 'interpolated value
Dim tmp() As Variant
Dim i As Long
' Get reference value index of x
With WorksheetFunction
If x < .Min(x_range) Then x = .Min(x_range) 'Constant value extrapolation in range under lower limit
x1_index = .Match(x, x_range, 1)
x2_index = .Min(x1_index + 1, x_range.Count) 'Constant value extrapolation in range over the upper limit
End With 'if x or y out of known value range, 1 and 2 indexes become same value
' Get reference values of x
x1 = x_range(x1_index)
x2 = x_range(x2_index)
' Get reference values of f
f1 = f_range(x1_index)
f2 = f_range(x2_index)
'Get interpolated values against x
If x2 <> x1 Then
f = f1 + (f2 - f1) * (x - x1) / (x2 - x1)
Else
f = f1
End If
Interpolation1D = f
End FunctionFunction InterpolationLinear(ByVal x As Double, ByVal x1 As Double, ByVal f1 As Double, ByVal x2 As Double, ByVal f2 As Double) As Double
'https://codepool.hatenablog.com/entry/2020/06/03/015921
'1:start, 2:end
Dim f As Double
Dim i As Long
'Get interpolated values against x
If x2 <> x1 Then
f = f1 + (f2 - f1) / (x2 - x1) * (x - x1)
Else
f = f1
End If
InterpolationLinear = f
End FunctionVLM
Sub VLM()
'シートの定義
Dim sht As Worksheet
Dim sht0 As Worksheet
'VLMの変数の定義
Dim alpha As Double '迎角 [deg]
Dim chord() As Variant 'コード長 [m]
Dim chord_cp() As Variant 'コード長 [m] 'コントロールポイントでのコード長 [m]
Dim cp() As Variant 'i番目の要素のcontrol point [m]
Dim sp1() As Variant 'i番目の要素のsurface point1 [m]
Dim sp2() As Variant 'i番目の要素のsurface point2 [m]
Dim x() As Variant 'x座標 [m]
Dim y() As Variant 'y座標 [m]
Dim z() As Variant 'z座標 [m]
Dim xc, zc As Double '翼型座標
Dim dzdx() As Variant 'キャンバーラインの傾き @cp [-]
Dim Cij() As Variant 'j番目のU字渦がi番目のcpに誘導する速度を表す係数
Dim dx1 As Double 'dx1=xi-x1j
Dim dx2 As Double 'dx2=xi-x2j
Dim dy1 As Double 'dy1=yi-y1j
Dim dy2 As Double 'dy2=yi-y2j
Dim r1 As Double 'r1=sqrt(dx1^2+dy1^2)
Dim r2 As Double 'r2=sqrt(dx2^2+dy2^2)
'翼の空力計算
Dim circulation() As Double '循環 [m^2/s]．Γc．Circulation
Dim circulation_chord() As Variant 'コード方向の循環の和 [m^2/s]．Γc．Circulation
Dim Lift As Double
Dim rho As Double
'一般
Dim sum As Double
Dim pi As Double: pi = 4 * Atn(1) '円周率
Dim rad As Double: rad = (pi / 180) 'degからradへの変換
Dim deg As Double: deg = (180 / pi) 'radからdegへの変換
'ガウスの消去法に代入するための配列
Dim matrix1() As Variant
Dim matrix2() As Variant
'カウンター
Dim nx, ny As Long 'コード方向分割数、セミスパン方向分割数
Dim i As Integer
Dim j As Integer
Dim k As Integer
Set sht = Worksheets("Wing Geometry")
Set sht0 = Worksheets("VLM")
alpha = 0 '迎角 [deg]
Sref = sht.Cells(3, 4) '面積 [m^2]
b = sht.Cells(4, 4) 'スパン [m]
cMAC = sht.Cells(6, 4) '平均空力翼弦長 [m]
lambda025 = sht.Cells(11, 4) '1/4コードライン後退角
Lambda = sht.Cells(13, 4) 'テーパー比
c0 = sht.Cells(14, 4) '翼根コード長 [m]
theta_r = sht.Cells(23, 4) '翼根取付角 [deg]
theta = sht.Cells(25, 4) 'ねじり下げ [deg]
tbarc = sht.Cells(32, 4) '翼厚比 [-]
V = Worksheets("Summary").Cells(2, 4) '対気速度 [m/s]
rho = Worksheets("Summary").Cells(20, 4) '空気密度 [kg/m^3]
'VLM
'********************
'read_value_from_sheet
Dim tmp As Variant
sht0.Range(sht0.Range(sht0.Cells(5, 5), sht0.Cells(5, 5).End(xlDown)), sht0.Range(sht0.Cells(5, 5), sht0.Cells(5, 5).End(xlDown)).End(xlToRight)).ClearContents
sht0.Range(sht0.Cells(2, 5), sht0.Cells(2, 5).End(xlToRight)).ClearContents
tmp = sht0.Range(sht0.Cells(5, 2), sht0.Cells(5, 2).End(xlDown))
ny = (UBound(tmp) - 1) / 2 'セミスパン方向分割数
nx = 16
ReDim chord(2 * ny)
ReDim chord_cp(2 * ny - 1)
For i = 0 To 2 * ny
chord(i) = c0 * (1 + (Lambda - 1) * Abs(sht0.Cells(5 + i, 4)))
Next i
For i = 0 To 2 * ny - 1
chord_cp(i) = (chord(i) + chord(i + 1)) / 2
Next i
'********************
'calculation_xyz
ReDim x((2 * ny + 1) * nx + 2 * ny)
ReDim y((2 * ny + 1) * nx + 2 * ny)
ReDim z((2 * ny + 1) * nx + 2 * ny)
rad = (4 * Atn(1)) / 180
'x,y,zを計算 '桁位置が原点 '左翼端が0，右翼端が2*span_divisions
For i = 0 To 2 * ny
For j = 0 To nx
num = i * nx + i + j
xc = (j / nx)
zc = 5 * tbarc * (0.2969 * Sqr(xc) - 0.126 * xc - 0.3516 * (xc ^ 2) + 0.2843 * (xc ^ 3) - 0.1015 * (xc ^ 4))
z(num) = chord(i) * zc
y(num) = sht0.Cells(5 + i, 4) * (b / 2)
x(num) = (Abs(y(num)) * Tan(rad * lambda025) - 0.25 * chord(i)) + chord(i) * xc
Next j
Next i
'********************
'calculation_cp_sp
ReDim cp(2, 2 * ny * nx - 1)
ReDim sp1(2, 2 * ny * nx - 1)
ReDim sp2(2, 2 * ny * nx - 1)
ReDim dzdx(2 * ny * nx - 1)
'sp,cp,dz/dxを計算
For i = 0 To 2 * ny - 1
For j = 0 To nx - 1
num1 = i * nx + j 'cp
num2 = i * nx + i + j 'xyz
'平面翼なのでz=0
cp(0, num1) = ((x(num2) * 1 + x(num2 + 1) * 3) / 4 + (x(num2 + nx + 1) * 1 + x(num2 + nx + 2) * 3) / 4) / 2
cp(1, num1) = ((y(num2) * 1 + y(num2 + 1) * 3) / 4 + (y(num2 + nx + 1) * 1 + y(num2 + nx + 2) * 3) / 4) / 2
cp(2, num1) = 0
sp1(0, num1) = (x(num2) * 3 + x(num2 + 1) * 1) / 4
sp1(1, num1) = (y(num2) * 3 + y(num2 + 1) * 1) / 4
sp1(2, num1) = 0
sp2(0, num1) = (x(num2 + nx + 1) * 3 + x(num2 + nx + 2) * 1) / 4
sp2(1, num1) = (y(num2 + nx + 1) * 3 + y(num2 + nx + 2) * 1) / 4
sp2(2, num1) = 0
dzdx(num1) = ((z(num2 + 1) + z(num2 + nx + 2)) / 2 - (z(num2) + z(num2 + nx + 1)) / 2) / ((x(num2 + 1) + x(num2 + nx + 2)) / 2 - (x(num2) + x(num2 + nx + 1)) / 2)
Next j
Next i
'********************
'calculation_Cij
ReDim Cij(2 * ny * nx - 1, 2 * ny * nx - 1)
'Cijを計算
For i = 0 To 2 * ny * nx - 1
For j = 0 To 2 * ny * nx - 1
dx1 = cp(0, i) - sp1(0, j) 'dx1=xi-x1j
dx2 = cp(0, i) - sp2(0, j) 'dx2=xi-x2j
dy1 = cp(1, i) - sp1(1, j) 'dy1=yi-y1j
dy2 = cp(1, i) - sp2(1, j) 'dy2=yi-y2j
r1 = Sqr(dx1 ^ 2 + dy1 ^ 2) 'r1=sqrt(dx1^2+dy1^2)
r2 = Sqr(dx2 ^ 2 + dy2 ^ 2) 'r2=sqrt(dx2^2+dy2^2)
Cij(i, j) = (1 / (dx1 * dy2 - dy1 * dx2)) * (((dx1 - dx2) * dx1 + (dy1 - dy2) * dy1) / r1 - ((dx1 - dx2) * dx2 + (dy1 - dy2) * dy2) / r2) - (1 / dy1) * (1 + dx1 / r1) + (1 / dy2) * (1 + dx2 / r2)
Next j
Next i
For k = 0 To 10
alpha = sht0.Cells(3, 5 + k)
'********************
'calculation_BC
'境界条件の方程式を計算する
ReDim matrix1(2 * ny * nx - 1, 2 * ny * nx - 1)
ReDim matrix2(2 * ny * nx - 1)
For i = 0 To 2 * ny * nx - 1
For j = 0 To 2 * ny * nx - 1
matrix1(i, j) = Cij(i, j) ' * cMAC
Next j
Next i
For i = 0 To 2 * ny * nx - 1
matrix2(i) = -4 * pi * V * (alpha * rad - dzdx(i) + (theta_r + theta * Abs(cp(1, i) / (b / 2))) * rad) ' * (2 / V)
Next i
'********************
'境界条件の方程式を解いて循環を求める
circulation = ans(2 * ny * nx, matrix1, matrix2)
'For i = 0 To 2 * ny - 1
' If pb2V <> 0 Then circulation(i) = circulation(i) / pb2V
'Next i
ReDim circulation_chord(2 * ny - 1)
For i = 0 To 2 * ny - 1
For j = 0 To nx - 1
num1 = i * nx + j 'cp
circulation_chord(i) = circulation_chord(i) + circulation(num1)
Next j
Next i
'********************
'clの計算
For i = 1 To 2 * ny - 1
sht0.Cells(5 + i, 5 + k) = 0.5 * (circulation_chord(i - 1) + circulation_chord(i)) * (2 / chord(i) / V)
Next i
sht0.Cells(5 + 0, 5 + k) = 0
sht0.Cells(5 + 2 * ny, 5 + k) = 0
Lift = 0
For i = 0 To 2 * ny - 1
Lift = Lift + rho * V * circulation_chord(i) * Abs(y(i * nx) - y((i + 1) * nx))
Next i
sht0.Cells(2, 5 + k) = Lift / (0.5 * rho * V * V * Sref)
Next k
sht0.Range("U3").GoalSeek Goal:=0, ChangingCell:=sht0.Range("T3")
End Sub行列計算
Public Function ans(ByVal n As Integer, ByRef A() As Variant, ByRef b() As Variant) As Double()
'ガウスの消去法によって，連立方程式Ax=bを解く．
'Aはn×n行列，x，bはn×1行列
'参考文献：「解析塾秘伝！有限要素法のつくり方」 p.82-
'変数の宣言
Dim pivot As Double 'マトリックスの体格成分
Dim p As Double '計算に使用するマトリックスの成分
Dim x() As Double '解を入れる配列
ReDim x(n - 1)
'カウンター
Dim r As Integer
Dim c As Integer
Dim rr As Integer
Dim cc As Integer
'前進消去
For r = 0 To n - 1
pivot = A(r, r) '対角成分をpivotに代入
For c = r To n - 1
A(r, c) = A(r, c) / pivot
Next c
b(r) = b(r) / pivot
For rr = r + 1 To n - 1
p = A(rr, r)
For cc = r To n - 1
A(rr, cc) = A(rr, cc) - p * A(r, cc)
Next cc
b(rr) = b(rr) - p * b(r)
Next rr
Next r
'後退代入
For r = n - 1 To 0 Step -1
x(r) = b(r)
For c = r + 1 To n - 1
x(r) = x(r) - A(r, c) * x(c)
Next c
Next r
ans = x
End Function
コメント