翼型解析③：翼の寸法を計算する

翼の座標データからリブの面積，周長，プランクの長さ，キャンバーライン，桁位置での翼厚を計算する

目的とフローチャート

今回のプログラムの目的は，Excelに入力した翼型の座標データからリブの面積，周長，プランクの長さ，キャンバーライン，桁位置での翼厚を計算し，キャンバーラインに関しては最小二乗法を用いて多項式近似することである

リブの面積，周長，プランクの長さは翼の重量推算，キャンバーラインはVLM(渦格子法)の計算，桁位置での翼厚は桁設計の際に桁が翼からはみ出さないことを確認するために用いる

プログラムのフローチャートを以下に示す

テンプレートファイルのダウンロード

エクセルシートの配置などの参考になるようにテンプレートファイルを添付しておく

※数式やマクロはすべて削除してあり，値も適当である

プログラムの解説

それでは実際に解説していく

使用する変数

マクロに使用する変数を以下に示す

最小二乗法に用いる変数については(≫翼型解析②：解析データを最小二乗法で多項式近似する)を参考にしてほしい

'変数の宣言
Dim last As Integer
Dim num As Integer '翼型の座標点数
Dim num_up As Integer '上面側の座標点数
Dim dif As Double 'x座標の差
Dim min As Double '最小値
Dim x() As Double 'x座標
Dim z() As Double 'z座標
Dim z_lerp As Double '線形補完されたz座標
Dim zc() As Double 'キャンバーラインのz座標
Dim Length As Double '翼の外周の長さ
Dim Plank As Double '翼のプランクを張る部分の外周の長さ
Dim start As Double '上面側のプランクの位置
Dim fin As Double '下面側のプランクの位置
Dim spar As Double '桁位置
Dim thickness As Double '桁位置での翼厚
Dim Area As Double 'リブの面積
'カウンター
Dim i As Integer
Dim j As Integer
Dim k As Integer '上面側のある点に1番近いx座標を持つ下面側の点の座標番号

値の読み込み

エクセルシートから必要な値を読み込む

設計者が入力するパラメータは上下面のプランク位置と桁位置である

'値読み込み
last = Cells(Rows.Count, 31).End(xlUp).Row '座標データの最終行を取得
num = last - 12 + 1 '座標点数
start = Range("AO17") '上面側のプランク位置
fin = Range("AO18") '下面側のプランク位置
spar = Range("AO19") '桁位置

翼の座標データの読み込み

読み込むべき座標点数がわかったので座標データの配列の大きさを宣言し値を読み込んでいく

このとき翼型のデータは後縁から始まり，上面を通って前縁を回り，後縁に戻っていくことに注意する

ReDim x(num - 1)
ReDim z(num - 1)

'座標データの読み込み
For i = 0 To num - 1 '後縁側から上面を通って1周するように
    x(i) = Cells(12 + i, 31)
    z(i) = Cells(12 + i, 32)
Next i

キャンバーラインと翼厚を計算

先に述べたように，翼型のデータは上下面の点数が異なる

そのため，まず上面のある点と一番x座標が近い下面の点を探し出し，そこから上面の点と同じx座標における下面でのz座標を計算する必要がある

上面の点と同じx座標における，下面のz座標の計算には線形補完を用いる

線形補完で求めた下面のz座標と上面のz座標を用いて，中心線の座標と翼厚を計算する

'上面と下面でx座標が少し違うので線形補完を行う
For i = 0 To num_up '後縁側から前縁側への上面のループ
    'キャンバーラインの計算．
    min = 1 '上下面のx座標の差の最大値
    k = 0
    '上面側のある点に1番近いx座標を持つ下面側の点を見つける
    For j = 0 To (num - 1) - num_up '後縁側から前縁側への下面のループ
        dif = x(num - 1 - j) - x(i) '上下面のx座標の差を計算(dif>0)
        If dif < 0 Then
            Exit For 'difが負になったらループを抜ける
        ElseIf min > dif Then '上下面のx座標の差が最小値よりも小さかったら
            min = dif '上下面のx座標の差の最小値を更新
            k = num - 1 - j '上下面のx座標の差が最小値となる下面側の点の座標番号
        End If
    Next j
    '下面のx(i)におけるz座標を線形補完で求める
    z_lerp = z(k - 1) + ((z(k) - z(k - 1)) / (x(k) - x(k - 1))) * (x(i) - x(k - 1))
    '中心線の座標を計算
    zc(i, 0) = (z(i) + z_lerp) / 2
    '翼厚の計算
    If x(i) > spar Then thickness = z(i) - z_lerp 'x座標が桁位置よりも大きかったら翼厚を計算
Next i

翼の周長とプランクの長さを計算

翼の周長と長さは，座標転換の距離を計算し，それを足し合わせていくことで計算する

'翼の外周の長さとプランクの長さを計算する
Length = 0
Plank = 0
For i = 0 To num - 2 '後縁側から上面を通って1周するように
    Length = Length + Sqr((x(i + 1) - x(i)) ^ 2 + (z(i + 1) - z(i)) ^ 2) '周長を計算
    If i < nup And x(i) < start Then '上面の点かつx座標がプランク開始位置より小さかったら
        Plank = Plank + Sqr((x(i + 1) - x(i)) ^ 2 + (z(i + 1) - z(i)) ^ 2)
    ElseIf x(i) < fin Then '下面の点かつx座標がプランク終了位置より小さかったら
        Plank = Plank + Sqr((x(i + 1) - x(i)) ^ 2 + (z(i + 1) - z(i)) ^ 2)
    End If
Next i

リブ面積を計算する

リブの面積は，リブを後縁側から小さな台形に分割したとして計算していく

実際は上底と下底は平行ではないので台形ではないが，小さく分割すれば平行に見えなくもないし，しょせん重量推算に使う推算値なので細かいことは気にしない

'リブの面積を計算する
Area = 0
For i = 0 To (num - 2) / 2 - 1 '後縁側から前縁側に
    'リブを小さな台形に分けて面積を計算する
    Area = Area + 0.5 * ((z(i + 1) - z(num - 1 - (i + 1))) + (z(i) - z(num - 1 - i))) * Abs(x(i + 1) - x(i))
Next i

キャンバーラインを最小二乗法で多項式近似する

詳しくは≫翼型解析②：解析データを最小二乗法で多項式近似するで述べているのでここではコードだけ紹介しておく

'多変数の最小二乗法
data1 = num_up

'配列の大きさの宣言
ReDim x(data1)

'値の読み込み
For i = 0 To num_up
    x(i) = Cells(12 + i, 31)
Next i

data2 = 8

ReDim matrix1(data1, data2)
ReDim matrix6(data1, 0)

'多変数関数の配列の作成
For i = 0 To data1
    For j = 0 To data2
        matrix1(i, j) = x(i) ^ j
    Next j
    matrix6(i, 0) = Cells(12 + i, 34)
Next i
matrix2 = Transpose(data1 + 1, data2 + 1, matrix1)
matrix3 = mult(matrix2, matrix1)
matrix4 = inverse(data2 + 1, matrix3)
matrix5 = mult(matrix4, matrix2)
matrix7 = mult(matrix5, matrix6)
For i = 0 To data2
    output(0, i) = matrix7(i, 0)
Next i

結果を出力

結果をシートに貼り付ける
翼厚の単位は[%]なので100倍しておく

'結果の出力
 Range("AH12:AH" & 12 + num_up & "") = zc()
 Range("AO12") = Length
 Range("AO13") = Plank
 Range("AO14") = Length - Plank
 Range("AO15") = Area
 Range("AO16") = thickness * 100
 Range(Cells(3, 29), Cells(3, 29 + data2)) = output()

ちなみに今回のプログラムで計算した長さや面積は，コード長が1である翼型データから計算されたものである

そのため任意のコード長の周長や面積は，今回計算した値にコード長，コード長の2乗をかけることで計算できる

プログラムの実行

今回作成したプログラムを実際に実行してみる

DAE31で計算してみたらこんな感じになる

まとめ

ようやくプログラミングらしくなってきた

↓次の記事

重量推算①：設計シートの下準備

エクセルシートを用いて機体の重量を推算するための下準備を行う

mtkbirdman.com

2019.11.04

↓記事一覧

設計シート作りの記事一覧

設計シートを作りについて説明した記事の一覧

mtkbirdman.com

2021.06.13