PR

DATCOMで安定微係数の計算を行うエクセルファイル

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 Function
Function 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 Function
Function 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 Function
Function 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 Function

VLM

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

DATCOM
質問・感想・意見などあれば気軽にTwitterのDMかコメントお願いします!
スポンサーリンク
mtk_birdman's blog

コメント