【FreeCAD】CSVファイルの圧力分布をもとにFEM解析を行うPythonスクリプト

XFLR5から出力した翼表面の圧力分布の.csvファイルをもとにFreeCADでFEM解析を行うPythonスクリプトについての解説する

概要
使用するデータの形式
フローチャートとソースコード全文
ソースコードの解説
使い方
おわりに

概要

最終的にはこんな感じのものを作成する

真ん中の穴は人力飛行機の翼構造をイメージしている

リブ全部付いた pic.twitter.com/e5a3jEgAUc
— いーそー (@mtk_birdman) May 10, 2019

↓FreeCADの公式ドキュメントはこちら

FreeCAD Documentation

インストールとチュートリアル

↓参考サイトと公式サイト

FreeCAD Windows でのインストール方法 - XSim

インストーラーを以下のサイトからダウンロードします……

FreeCAD: Your own 3D parametric modeler

FreeCAD, the open source 3D parametric modeler

使い方は下のサイトの「形状作成チュートリアル」を一通りやれば完全に理解できる

FreeCAD 使い方メモ - XSim

オープンソースの汎用 3D-CAD である FreeCAD の使い方のメモです……

使用するデータの形式

FreeCADにはデフォルトで2次元翼型読み込みの機能が入っている

Common Airfoil Data Import - FreeCAD Documentation

FreeCAD can import airfoil data such as that found on the UIUC Airfoil Coordinates Database or files produced by airfoil creation and annalizing software like XFLR5.
https://wiki.freecadweb.org/Common_Airfoil_Data_Import

公式が2次元翼型の読み込みにXFLR5の形式を採用しているので，この記事でも圧力分布の読み込みにXFLR5の形式を使うことにする

↓XFLR5についてはこちら

NACA4412

今回の記事で使うのはXFLR5で解析したこの翼型

この圧力分布をXFLR5で出力すると次のような.csvファイルができる

Cp_Graph.csv

X,NACA4412-Re=   1e+06-Alpha= 5.35-NCrit=  9.0-XTrTop=0.000-XtrBot=1.000, , 
            1,     0.1400708, , 
    0.9928006,     0.1306438, , 
    0.9798895,     0.1109444, , 
    0.9635184,    0.08210746, , 
    0.9445469,     0.0446714, , 
    0.9234998,  4.864211e-05, , 
    0.9007495,   -0.04960863, , 
     0.876583,    -0.1020141, , 
    0.8512328,    -0.1552933, , 
    0.8248944,    -0.2082254, , 
    0.7977362,    -0.2601881, , 
    0.7699063,     -0.310967, , 
    0.7415373,    -0.3606241, , 
    0.7127487,    -0.4093158, , 
      0.68365,    -0.4572609, , 
    0.6543419,     -0.504709, , 
    0.6249182,    -0.5518682, , 
    0.5954661,     -0.599014, , 
    0.5660677,    -0.6463763, , 
    0.5368005,    -0.6942692, , 
    0.5077378,    -0.7430246, , 
    0.4789491,    -0.7931529, , 
    0.4505011,    -0.8456015, , 
    0.4224571,    -0.9012967, , 
    0.3948782,    -0.9703361, , 
    0.3678229,     -1.032724, , 
    0.3413476,     -1.084167, , 
    0.3155065,     -1.129935, , 
    0.2903524,     -1.170667, , 
    0.2659361,     -1.206881, , 
    0.2423068,     -1.239076, , 
    0.2195125,     -1.267461, , 
    0.1975998,     -1.292415, , 
    0.1766138,     -1.314321, , 
    0.1565988,     -1.333565, , 
    0.1375977,     -1.350631, , 
    0.1196525,     -1.366144, , 
    0.1028043,     -1.380608, , 
    0.0870931,     -1.394651, , 
    0.0725582,     -1.409189, , 
    0.0592379,     -1.425967, , 
      0.04717,     -1.445794, , 
    0.0398384,      -1.46135, , 
    0.0363913,     -1.468926, , 
    0.0330915,     -1.476358, , 
    0.0299394,      -1.48357, , 
    0.0269381,     -1.490287, , 
     0.024088,     -1.496443, , 
    0.0213896,     -1.501661, , 
    0.0188458,      -1.50578, , 
    0.0164571,     -1.507067, , 
    0.0142241,     -1.505324, , 
    0.0121494,     -1.499479, , 
    0.0102336,     -1.487902, , 
    0.0084774,     -1.468712, , 
    0.0068832,      -1.43877, , 
    0.0054516,     -1.395099, , 
    0.0041835,     -1.334336, , 
    0.0030809,     -1.252873, , 
    0.0021446,     -1.146472, , 
    0.0013756,     -1.013646, , 
    0.0007756,    -0.8451148, , 
    0.0003455,    -0.6391955, , 
     8.66e-05,    -0.4096659, , 
            0,    -0.1653529, , 
     8.66e-05,    0.08356277, , 
    0.0003455,     0.3209784, , 
    0.0007756,     0.5341777, , 
    0.0013756,     0.7114345, , 
    0.0021446,     0.8456459, , 
    0.0030809,     0.9348456, , 
    0.0041835,     0.9839083, , 
    0.0054516,     0.9998923, , 
    0.0068832,     0.9911922, , 
    0.0084774,     0.9651613, , 
    0.0102336,      0.928468, , 
    0.0121494,     0.8853146, , 
    0.0142241,     0.8397434, , 
    0.0164571,     0.7934297, , 
    0.0188458,      0.748363, , 
    0.0213896,     0.7052552, , 
     0.024088,     0.6646386, , 
    0.0269381,     0.6267616, , 
    0.0299394,      0.591703, , 
    0.0330915,     0.5594039, , 
    0.0363913,     0.5297343, , 
    0.0398384,     0.5022834, , 
      0.04717,     0.4539824, , 
    0.0592379,     0.3972975, , 
    0.0725582,     0.3552181, , 
    0.0870931,     0.3236172, , 
    0.1028043,     0.3000543, , 
    0.1196525,     0.2826884, , 
    0.1375977,     0.2700834, , 
    0.1565988,     0.2611257, , 
    0.1766138,     0.2549584, , 
    0.1975998,     0.2508825, , 
    0.2195125,     0.2482936, , 
    0.2423068,     0.2467346, , 
    0.2659361,     0.2457315, , 
    0.2903524,     0.2448956, , 
    0.3155065,     0.2437871, , 
    0.3413476,     0.2419036, , 
    0.3678229,     0.2386706, , 
    0.3948782,     0.2318488, , 
    0.4224571,     0.2217901, , 
    0.4505011,     0.2161533, , 
    0.4789491,     0.2122116, , 
    0.5077378,     0.2093985, , 
    0.5368005,     0.2071908, , 
    0.5660677,     0.2052845, , 
    0.5954661,     0.2034736, , 
    0.6249182,     0.2016096, , 
    0.6543419,     0.1995884, , 
      0.68365,      0.197338, , 
    0.7127487,     0.1947957, , 
    0.7415373,      0.191957, , 
    0.7699063,     0.1887996, , 
    0.7977362,     0.1853389, , 
    0.8248944,     0.1816187, , 
    0.8512328,     0.1776779, , 
     0.876583,     0.1736001, , 
    0.9007495,     0.1694696, , 
    0.9234998,     0.1654165, , 
    0.9445469,     0.1615532, , 
    0.9635184,      0.158121, , 
    0.9798895,     0.1548783, , 
    0.9928006,     0.1534044, , 
            1,     0.1400708, ,

データの出力は以下の手順で行う

①「XFoil Direct Analysis」の「Operationg Points>Cp Graph>Export Graph」をクリック

②好きな名前で保存

以上

各データは左から順に以下の通り

名称	単位
Operationg Pointの位置（x座標）	m
圧力係数（Cp）	-

.csvファイルのデータ構成

2次元翼のデータはこんな感じ

NACA4412.dat

NACA4412
 1.0000000     0.0012600
 0.9928006     0.0032198
 0.9798895     0.0066835
 0.9635184     0.0109823
 0.9445469     0.0158354
 0.9234998     0.0210604
 0.9007495     0.0265227
 0.8765830     0.0321168
 0.8512328     0.0377572
 0.8248944     0.0433730
 0.7977362     0.0489045
 0.7699063     0.0543008
 0.7415373     0.0595184
 0.7127487     0.0645195
 0.6836500     0.0692712
 0.6543419     0.0737448
 0.6249182     0.0779148
 0.5954661     0.0817592
 0.5660677     0.0852584
 0.5368005     0.0883955
 0.5077378     0.0911557
 0.4789491     0.0935267
 0.4505011     0.0954983
 0.4224571     0.0970626
 0.3948782     0.0982103
 0.3678229     0.0988053
 0.3413476     0.0987893
 0.3155065     0.0981781
 0.2903524     0.0969911
 0.2659361     0.0952510
 0.2423068     0.0929845
 0.2195125     0.0902210
 0.1975998     0.0869933
 0.1766138     0.0833369
 0.1565988     0.0792897
 0.1375977     0.0748918
 0.1196525     0.0701852
 0.1028043     0.0652129
 0.0870931     0.0600189
 0.0725582     0.0546478
 0.0592379     0.0491435
 0.0471700     0.0435495
 0.0398384     0.0397909
 0.0363913     0.0379076
 0.0330915     0.0360235
 0.0299394     0.0341394
 0.0269381     0.0322577
 0.0240880     0.0303793
 0.0213896     0.0285049
 0.0188458     0.0266369
 0.0164571     0.0247759
 0.0142241     0.0229225
 0.0121494     0.0210790
 0.0102336     0.0192459
 0.0084774     0.0174236
 0.0068832     0.0156139
 0.0054516     0.0138171
 0.0041835     0.0120339
 0.0030809     0.0102667
 0.0021446     0.0085156
 0.0013756     0.0067792
 0.0007756     0.0050573
 0.0003455     0.0033493
 0.0000866     0.0016607
 0.0000000     0.0000000
 0.0000866    -0.0016261
 0.0003455    -0.0032112
 0.0007756    -0.0047474
 0.0013756    -0.0062299
 0.0021446    -0.0076600
 0.0030809    -0.0090391
 0.0041835    -0.0103693
 0.0054516    -0.0116513
 0.0068832    -0.0128843
 0.0084774    -0.0140686
 0.0102336    -0.0152049
 0.0121494    -0.0162931
 0.0142241    -0.0173340
 0.0164571    -0.0183285
 0.0188458    -0.0192762
 0.0213896    -0.0201778
 0.0240880    -0.0210342
 0.0269381    -0.0218453
 0.0299394    -0.0226118
 0.0330915    -0.0233344
 0.0363913    -0.0240133
 0.0398384    -0.0246491
 0.0471700    -0.0257940
 0.0592379    -0.0272029
 0.0725582    -0.0282569
 0.0870931    -0.0289743
 0.1028043    -0.0293755
 0.1196525    -0.0294825
 0.1375977    -0.0293193
 0.1565988    -0.0289118
 0.1766138    -0.0282876
 0.1975998    -0.0274762
 0.2195125    -0.0265089
 0.2423068    -0.0254180
 0.2659361    -0.0242376
 0.2903524    -0.0230023
 0.3155065    -0.0217477
 0.3413476    -0.0205094
 0.3678229    -0.0193230
 0.3948782    -0.0182234
 0.4224571    -0.0171746
 0.4505011    -0.0160650
 0.4789491    -0.0149118
 0.5077378    -0.0137351
 0.5368005    -0.0125542
 0.5660677    -0.0113870
 0.5954661    -0.0102497
 0.6249182    -0.0091567
 0.6543419    -0.0081203
 0.6836500    -0.0071506
 0.7127487    -0.0062555
 0.7415373    -0.0054401
 0.7699063    -0.0047076
 0.7977362    -0.0040587
 0.8248944    -0.0034919
 0.8512328    -0.0030041
 0.8765830    -0.0025904
 0.9007495    -0.0022449
 0.9234998    -0.0019608
 0.9445469    -0.0017313
 0.9635184    -0.0015496
 0.9798895    -0.0014106
 0.9928006    -0.0013115
 1.0000000    -0.0012600

NACA0012
 1.0000000     0.0012600
 0.9928006     0.0022657
 0.9798895     0.0040471
 0.9635184     0.0062660
 0.9445469     0.0087833
 0.9234998     0.0115106
 0.9007495     0.0143838
 0.8765830     0.0173536
 0.8512328     0.0203807
 0.8248944     0.0234325
 0.7977362     0.0264816
 0.7699063     0.0295042
 0.7415373     0.0324793
 0.7127487     0.0353875
 0.6836500     0.0382109
 0.6543419     0.0409325
 0.6249182     0.0435357
 0.5954661     0.0460044
 0.5660677     0.0483227
 0.5368005     0.0504749
 0.5077378     0.0524454
 0.4789491     0.0542192
 0.4505011     0.0557817
 0.4224571     0.0571186
 0.3948782     0.0582169
 0.3678229     0.0590641
 0.3413476     0.0596494
 0.3155065     0.0599629
 0.2903524     0.0599967
 0.2659361     0.0597443
 0.2423068     0.0592012
 0.2195125     0.0583649
 0.1975998     0.0572348
 0.1766138     0.0558122
 0.1565988     0.0541007
 0.1375977     0.0521056
 0.1196525     0.0498338
 0.1028043     0.0472942
 0.0870931     0.0444966
 0.0725582     0.0414523
 0.0592379     0.0381732
 0.0471700     0.0346718
 0.0363913     0.0309604
 0.0299394     0.0283756
 0.0269381     0.0270515
 0.0240880     0.0257068
 0.0213896     0.0243414
 0.0188458     0.0229565
 0.0164571     0.0215522
 0.0142241     0.0201283
 0.0121494     0.0186860
 0.0102336     0.0172254
 0.0084774     0.0157461
 0.0068832     0.0142491
 0.0054516     0.0127342
 0.0041835     0.0112016
 0.0030809     0.0096529
 0.0021446     0.0080878
 0.0013756     0.0065046
 0.0007756     0.0049023
 0.0003455     0.0032802
 0.0000866     0.0016434
 0.0000000     0.0000000
 0.0000866    -0.0016434
 0.0003455    -0.0032802
 0.0007756    -0.0049023
 0.0013756    -0.0065046
 0.0021446    -0.0080878
 0.0030809    -0.0096529
 0.0041835    -0.0112016
 0.0054516    -0.0127342
 0.0068832    -0.0142491
 0.0084774    -0.0157461
 0.0102336    -0.0172254
 0.0121494    -0.0186860
 0.0142241    -0.0201283
 0.0164571    -0.0215522
 0.0188458    -0.0229565
 0.0213896    -0.0243414
 0.0240880    -0.0257068
 0.0269381    -0.0270515
 0.0299394    -0.0283756
 0.0363913    -0.0309604
 0.0471700    -0.0346718
 0.0592379    -0.0381732
 0.0725582    -0.0414523
 0.0870931    -0.0444966
 0.1028043    -0.0472942
 0.1196525    -0.0498338
 0.1375977    -0.0521056
 0.1565988    -0.0541007
 0.1766138    -0.0558122
 0.1975998    -0.0572348
 0.2195125    -0.0583649
 0.2423068    -0.0592012
 0.2659361    -0.0597443
 0.2903524    -0.0599967
 0.3155065    -0.0599629
 0.3413476    -0.0596494
 0.3678229    -0.0590641
 0.3948782    -0.0582169
 0.4224571    -0.0571186
 0.4505011    -0.0557817
 0.4789491    -0.0542192
 0.5077378    -0.0524454
 0.5368005    -0.0504749
 0.5660677    -0.0483227
 0.5954661    -0.0460044
 0.6249182    -0.0435357
 0.6543419    -0.0409325
 0.6836500    -0.0382109
 0.7127487    -0.0353875
 0.7415373    -0.0324793
 0.7699063    -0.0295042
 0.7977362    -0.0264816
 0.8248944    -0.0234325
 0.8512328    -0.0203807
 0.8765830    -0.0173536
 0.9007495    -0.0143838
 0.9234998    -0.0115106
 0.9445469    -0.0087833
 0.9635184    -0.0062660
 0.9798895    -0.0040471
 0.9928006    -0.0022657
 1.0000000    -0.0012600

NACA0009
 1.0000000     0.0009450
 0.9901282     0.0019773
 0.9701682     0.0040276
 0.9469539     0.0063503
 0.9229830     0.0086826
 0.8988830     0.0109621
 0.8747544     0.0131812
 0.8506129     0.0153402
 0.8264621     0.0174401
 0.8023039     0.0194817
 0.7781400     0.0214655
 0.7539718     0.0233914
 0.7298010     0.0252589
 0.7056291     0.0270669
 0.6814577     0.0288141
 0.6572885     0.0304984
 0.6331231     0.0321174
 0.6089632     0.0336680
 0.5848106     0.0351466
 0.5606672     0.0365489
 0.5365351     0.0378702
 0.5124163     0.0391050
 0.4883132     0.0402470
 0.4642285     0.0412893
 0.4401651     0.0422242
 0.4161262     0.0430429
 0.3921157     0.0437356
 0.3681380     0.0442917
 0.3441984     0.0446990
 0.3203032     0.0449439
 0.2964603     0.0450110
 0.2726795     0.0448829
 0.2489733     0.0445397
 0.2253582     0.0439580
 0.2018562     0.0431109
 0.1784983     0.0419661
 0.1553293     0.0404845
 0.1324180     0.0386185
 0.1098794     0.0363101
 0.0879308     0.0334933
 0.0670466     0.0301197
 0.0482645     0.0262609
 0.0330984     0.0222613
 0.0253681     0.0197423
 0.0221517     0.0185528
 0.0193096     0.0174133
 0.0167968     0.0163198
 0.0145701     0.0152677
 0.0125870     0.0142513
 0.0108176     0.0132648
 0.0092361     0.0123031
 0.0078178     0.0113606
 0.0065463     0.0104322
 0.0054085     0.0095141
 0.0043917     0.0086017
 0.0034888     0.0076918
 0.0026950     0.0067823
 0.0020056     0.0058707
 0.0014197     0.0049561
 0.0009380     0.0040392
 0.0005596     0.0031211
 0.0002811     0.0022089
 0.0000981     0.0013119
 0.0000069     0.0004347
 0.0000069    -0.0004347
 0.0000981    -0.0013119
 0.0002811    -0.0022089
 0.0005596    -0.0031211
 0.0009380    -0.0040392
 0.0014197    -0.0049561
 0.0020056    -0.0058707
 0.0026950    -0.0067823
 0.0034888    -0.0076918
 0.0043917    -0.0086017
 0.0054085    -0.0095141
 0.0065463    -0.0104322
 0.0078178    -0.0113606
 0.0092361    -0.0123031
 0.0108176    -0.0132648
 0.0125870    -0.0142513
 0.0145701    -0.0152677
 0.0167968    -0.0163198
 0.0193096    -0.0174133
 0.0221517    -0.0185528
 0.0253681    -0.0197423
 0.0330984    -0.0222613
 0.0482645    -0.0262609
 0.0670466    -0.0301197
 0.0879308    -0.0334933
 0.1098794    -0.0363101
 0.1324180    -0.0386185
 0.1553293    -0.0404845
 0.1784983    -0.0419661
 0.2018562    -0.0431109
 0.2253582    -0.0439580
 0.2489733    -0.0445397
 0.2726795    -0.0448829
 0.2964603    -0.0450110
 0.3203032    -0.0449439
 0.3441984    -0.0446990
 0.3681380    -0.0442917
 0.3921157    -0.0437356
 0.4161262    -0.0430429
 0.4401651    -0.0422242
 0.4642285    -0.0412893
 0.4883132    -0.0402470
 0.5124163    -0.0391050
 0.5365351    -0.0378702
 0.5606672    -0.0365489
 0.5848106    -0.0351466
 0.6089632    -0.0336680
 0.6331231    -0.0321174
 0.6572885    -0.0304984
 0.6814577    -0.0288141
 0.7056291    -0.0270669
 0.7298010    -0.0252589
 0.7539718    -0.0233914
 0.7781400    -0.0214655
 0.8023039    -0.0194817
 0.8264621    -0.0174401
 0.8506129    -0.0153402
 0.8747544    -0.0131812
 0.8988830    -0.0109621
 0.9229830    -0.0086826
 0.9469539    -0.0063503
 0.9701682    -0.0040276
 0.9901282    -0.0019773
 1.0000000    -0.0009450

これら2つのテキストファイル（.csvと.dat）から冒頭に示したようなFEM解析を行う

フローチャートとソースコード全文

今回のプログラムのフローチャートはこんな感じ

ソースコードはこちら

airfoil_Cp.py

import os
import PySide2
import csv
import numpy as np

import FreeCAD as App
import FreeCADGui as Gui
import importAirfoilDAT
import ObjectsFem
import FemGui
import Draft
import Part

### Read .dat
DirPath = os.path.dirname(App.ActiveDocument.FileName)
ReadName, Filter = PySide2.QtWidgets.QFileDialog.getOpenFileName(None, "Read a file", DirPath, "All files (*.*);;") # Open dialog box
with open(ReadName) as f: # Open .dat
    dat = f.readlines()
    dat = dat[1:]
dat = [line.split() for line in dat]
dat = [[float(str) for str in line] for line in dat]

# Calcurate max thickness and its position
chord = 1000
thickness = 0
x = np.array(dat)[:,0]
y = np.array(dat)[:,1]
LE_idx = np.argmin(np.abs(x-0.0))
for i in range(int(len(dat)/2)):
    x1 = x[i]
    y1 = y[i]
    idx = LE_idx+np.argmin(np.abs(x[LE_idx:]-x1))
    y2 = y[idx]
    if abs(y1-y2)*chord>thickness:
        thickness = abs(y1-y2)*chord # Max thickness
        x0 = x1*chord # x coordinate of the center of the pocket
        y0 = y2*chord+thickness/2.0 # y coordinate of the center of the pocket

### Import airfoil sections
FoilName = os.path.basename(ReadName).strip('.dat')
importAirfoilDAT.insert(DirPath+"/"+FoilName+".dat",App.ActiveDocument.Name) # Import .dat file
AirfoilObject = App.ActiveDocument.getObject(FoilName).Group[0] # Get the imported airfoil object
AirfoilObject.Label = AirfoilObject.Name+"_"+FoilName # Set label for the object
CloneObject2D = Draft.clone(AirfoilObject) # Create clones for each aifoil sections
CloneObject2D.Scale = (chord,chord,chord) # Move and Scale the cloned object

### Begin command PartDesign_Clone
App.ActiveDocument.addObject('PartDesign::Body','Body')
App.ActiveDocument.addObject('PartDesign::FeatureBase','Clone')
BodyObject = App.ActiveDocument.getObject('Body')
CloneObject = App.ActiveDocument.getObject('Clone')
CloneObject.BaseFeature = CloneObject2D
CloneObject.Placement = CloneObject2D.Placement
CloneObject.setEditorMode('Placement',0)
BodyObject.Group = [CloneObject]
BodyObject.Tip = CloneObject
CloneObject.ViewObject.Transparency=getattr(CloneObject2D.getLinkedObject(True).ViewObject,'Transparency',CloneObject.ViewObject.Transparency)
CloneObject.ViewObject.DisplayMode=getattr(CloneObject2D.getLinkedObject(True).ViewObject,'DisplayMode',CloneObject.ViewObject.DisplayMode)
### End command PartDesign_Clone

### Begin command PartDesign_Pad
Gui.ActiveDocument.ActiveView.setActiveObject('pdbody',BodyObject,'')
BodyObject.newObject('PartDesign::Pad','Pad')
PadObject = App.ActiveDocument.getObject('Pad')
PadObject.Profile = (CloneObject, ['Face1',])
PadObject.Length = 10.0
PadObject.ViewObject.ShapeColor=getattr(CloneObject.getLinkedObject(True).ViewObject,'ShapeColor',PadObject.ViewObject.ShapeColor)
PadObject.ViewObject.LineColor=getattr(CloneObject.getLinkedObject(True).ViewObject,'LineColor',PadObject.ViewObject.LineColor)
PadObject.ViewObject.PointColor=getattr(CloneObject.getLinkedObject(True).ViewObject,'PointColor',PadObject.ViewObject.PointColor)
PadObject.ViewObject.Transparency=getattr(CloneObject.getLinkedObject(True).ViewObject,'Transparency',PadObject.ViewObject.Transparency)
PadObject.ViewObject.DisplayMode=getattr(CloneObject.getLinkedObject(True).ViewObject,'DisplayMode',PadObject.ViewObject.DisplayMode)
Gui.ActiveDocument.setEdit(BodyObject,0,'Pad.')
Gui.Selection.clearSelection()
### End command PartDesign_Pad
App.ActiveDocument.recompute()

### Begin command PartDesign_NewSketch
BodyObject.newObject('Sketcher::SketchObject','Sketch')
SketchObject = App.ActiveDocument.getObject('Sketch')
SketchObject.Support = (PadObject,['Face'+str(len(PadObject.Shape.Faces)),])
SketchObject.MapMode = 'FlatFace'
### End command PartDesign_NewSketch
Gui.runCommand('Sketcher_CompCreateCircle',0)
SketchObject.addGeometry(Part.Circle(App.Vector(x0,y0,0),App.Vector(0,0,1),thickness*0.9/2),False)

### Begin command PartDesign_Pocket
BodyObject.newObject('PartDesign::Pocket','Pocket')
PocketObject =App.ActiveDocument.getObject('Pocket')
PocketObject.Profile = SketchObject
PocketObject.Length = 100.000000
PocketObject.Length2 = 100.000000
PocketObject.Type = 0
PocketObject.UpToFace = None
PocketObject.Reversed = 0
PocketObject.Midplane = 0
PocketObject.Offset = 0
PocketObject.ViewObject.ShapeColor=getattr(PadObject.getLinkedObject(True).ViewObject,'ShapeColor',PocketObject.ViewObject.ShapeColor)
PocketObject.ViewObject.LineColor=getattr(PadObject.getLinkedObject(True).ViewObject,'LineColor',PocketObject.ViewObject.LineColor)
PocketObject.ViewObject.PointColor=getattr(PadObject.getLinkedObject(True).ViewObject,'PointColor',PocketObject.ViewObject.PointColor)
PocketObject.ViewObject.Transparency=getattr(PadObject.getLinkedObject(True).ViewObject,'Transparency',PocketObject.ViewObject.Transparency)
PocketObject.ViewObject.DisplayMode=getattr(PadObject.getLinkedObject(True).ViewObject,'DisplayMode',PocketObject.ViewObject.DisplayMode)
### End command PartDesign_Pocket
App.ActiveDocument.recompute()
Gui.ActiveDocument.resetEdit()

AirfoilObject.Visibility = False
CloneObject2D.Visibility = False
CloneObject.Visibility = False
PadObject.Visibility = False
SketchObject.Visibility = False

### Read .csv
DirPath = os.path.dirname(App.ActiveDocument.FileName)
ReadName, Filter = PySide2.QtWidgets.QFileDialog.getOpenFileName(None, "Read a file", DirPath, "All files (*.*);;") # Open dialog box
with open(ReadName) as f: # Open .csv
    csv = csv.reader(f)
    csv = [line for line in csv]
    csv = csv[1:-2]
csv = [line for line in csv]
csv = [[float(str) for str in line[:2]] for line in csv]
xCp = np.array(csv)[:,0]
Cp = np.array(csv)[:,1]
LE_idx = np.argmin(np.abs(xCp-0.0))

### Begin command FEM_Analysis
ObjectsFem.makeAnalysis(App.ActiveDocument, 'Analysis')
FemGui.setActiveAnalysis(App.ActiveDocument.ActiveObject)
ObjectsFem.makeSolverCalculixCcxTools(App.ActiveDocument)
FemGui.getActiveAnalysis().addObject(App.ActiveDocument.ActiveObject)
### End command FEM_Analysis

### Begin command FEM_ConstraintFixed
App.ActiveDocument.addObject("Fem::ConstraintFixed","ConstraintFixed")
App.ActiveDocument.ConstraintFixed.Scale = 1
App.ActiveDocument.ConstraintFixed.References = [(App.ActiveDocument.Pocket,"Face"+str(len(PocketObject.Shape.Faces)))]
App.ActiveDocument.Analysis.addObject(App.ActiveDocument.ConstraintFixed)
for amesh in App.ActiveDocument.Objects:
    if "ConstraintFixed" == amesh.Name:
        amesh.ViewObject.Visibility = True
    elif "Mesh" in amesh.TypeId:
        aparttoshow = amesh.Name.replace("_Mesh","")
        for apart in App.ActiveDocument.Objects:
            if aparttoshow == apart.Name:
                apart.ViewObject.Visibility = True
        amesh.ViewObject.Visibility = False
Gui.ActiveDocument.setEdit('ConstraintFixed')
### End command FEM_ConstraintFixed

face_idx = 1
for face in PocketObject.Shape.Faces:
    if face.Surface.Axis.z == 0 and not face.Surface.Axis.y == 0:
        if face.Surface.Axis.y<0: # Upper surface
            idx = np.argmin(np.abs(xCp[:LE_idx]*chord-face.CenterOfMass.x))
        else: # Lower surface
            idx = LE_idx+np.argmin(np.abs(xCp[LE_idx:]*chord-face.CenterOfMass.x))
        
        ### Begin command FEM_ConstraintPressure
        ObjectName = "ConstraintPressure"+str(face_idx).zfill(3)
        App.ActiveDocument.addObject("Fem::ConstraintPressure",ObjectName)
        ConstraintPressureObject = App.ActiveDocument.getObject(ObjectName)
        ConstraintPressureObject.Pressure = abs(Cp[idx]) # MPa
        ConstraintPressureObject.Reversed = bool(Cp[idx]<0) # Positive/Negative pressure
        ConstraintPressureObject.Scale = 10
        ConstraintPressureObject.References = [(App.ActiveDocument.Pocket,"Face"+str(face_idx))]
        App.ActiveDocument.Analysis.addObject(ConstraintPressureObject)
        for amesh in App.ActiveDocument.Objects:
            if ObjectName == amesh.Name:
                amesh.ViewObject.Visibility = True
            elif "Mesh" in amesh.TypeId:
                aparttoshow = amesh.Name.replace("_Mesh","")
                for apart in App.ActiveDocument.Objects:
                    if aparttoshow == apart.Name:
                        apart.ViewObject.Visibility = True
                amesh.ViewObject.Visibility = False

        print("Face"+str(face_idx),Cp[idx],face.Surface.Axis.y > 0)
        ### End command FEM_ConstraintPressure
    
    face_idx += 1

Gui.ActiveDocument.resetEdit()
App.ActiveDocument.recompute()

これだけ見てもなんのこっちゃ分からないので，1つずつ解説していく

ソースコードの解説

FreeCADにおけるプログラムの作成方法

FreeCADで行うすべての操作はPythonのスクリプトに変換されて実行される

Python Consoleを開けば自分が行った操作がPythonでどのように記述されるかをリアルタイムで知ることができる

Python console - FreeCAD Documentation

自分で何かしらのマクロを作成したい場合は

マクロにしたい操作を自分で実行する
Python Consoleに出てくるスクリプトをコピペする
上記のスクリプトをもとに汎化性を持たせたコードを作成する

といった手順をとることで，比較的簡単に欲しいマクロを作成することができる

作成したプログラムのデバッグを行う際はReport Viewを使うといい

Report view - FreeCAD Documentation

print()の出力やエラー文などが表示される

以上のことを踏まえて，フローチャートに沿って説明していく

モジュールのインポート

今回のプログラムでは以下のモジュールを使用する

import os
import PySide2
import csv
import numpy as np

import FreeCAD as App
import FreeCADGui as Gui
import importAirfoilDAT
import ObjectsFem
import FemGui
import Draft
import Part

PySide2は「ファイルを開く」ダイアログボックスを使うのに使用する

QFileDialog — Qt for Python

FreeCADとFreeCADGuiはPython Consoleに合わせてそれぞれApp，Guiとしてインポートしている

.datの読み込み

.datファイルを読み込む

### Read .dat
DirPath = os.path.dirname(App.ActiveDocument.FileName)
ReadName, Filter = PySide2.QtWidgets.QFileDialog.getOpenFileName(None, "Read a file", DirPath, "All files (*.*);;") # Open dialog box
with open(ReadName) as f: # Open .dat
    dat = f.readlines()
    dat = dat[1:]
dat = [line.split() for line in dat]
dat = [[float(str) for str in line] for line in dat]

App.ActiveDocument.FileNameで現在編集しているFreeCADのファイル名を取得し，os.path.dirname()でFreeCADのファイルが保存されているディレクトリのPathを取得する

PySide2.QtWidgets.QFileDialog.getOpenFileName()で，「ファイルを開く」のダイアログボックスを開くことができる

最初に開くのは上で取得したFreeCADのファイルが保存されているディレクトリである

選択したファイルはwith open()以降で読み込まれている

.datファイルはスペース区切りのテキストファイルなので，split()を使って空白を削除してリストに変換している

参考
≫FreeCAD API - FreeCAD Documentation
≫Pythonで実行中のファイルの場所（パス）を取得する
≫QFileDialog — Qt for Python
≫Pythonでファイルの読み込み、書き込み（作成・追記）
≫Python, splitでカンマ区切り文字列を分割、空白を削除しリスト化

翼厚と最大翼厚位置の計算

翼厚と最大翼厚位置を計算する

# Calcurate max thickness and its position

chord = 1000
thickness = 0
x = np.array(dat)[:,0]
y = np.array(dat)[:,1]
LE_idx = np.argmin(np.abs(x-0.0))
for i in range(int(len(dat)/2)):
    x1 = x[i]
    y1 = y[i]
    idx = LE_idx+np.argmin(np.abs(x[LE_idx:]-x1))
    y2 = y[idx]
    if abs(y1-y2)*chord>thickness:
        thickness = abs(y1-y2)*chord # Max thickness
        x0 = x1*chord # x coordinate of the center of the pocket
        y0 = y2*chord+thickness/2.0 # y coordinate of the center of the pocket

np.array()を使って読み込んだ.datをnumpy配列に変換し，np.argmin()を使って前縁点のインデックス（LE_idx）を取得している

後縁から上面を通るループを回し，上面の点と最も近いx座標を持つ下面の点を見つけ，最大翼厚を計算している

x0とy0は後で翼にあける穴の中心座標である

参考
≫NumPy配列ndarrayとPython標準のリストを相互に変換 _ note.nkmk.me
≫NumPy配列ndarrayの最大値・最小値のインデックス（位置）を取得 _ note.nkmk.me

翼型のインポート

.datファイルをインポートし，Cloneオブジェクトを作成する

### Import airfoil sections
FoilName = os.path.basename(ReadName).strip('.dat')
importAirfoilDAT.insert(DirPath+"/"+FoilName+".dat",App.ActiveDocument.Name) # Import .dat file
AirfoilObject = App.ActiveDocument.getObject(FoilName).Group[0] # Get the imported airfoil object
AirfoilObject.Label = AirfoilObject.Name+"_"+FoilName # Set label for the object
CloneObject2D = Draft.clone(AirfoilObject) # Create clones for each aifoil sections
CloneObject2D.Scale = (chord,chord,chord) # Move and Scale the cloned object

importAirfoilDAT.insert()を使って翼型の.datファイルをインポートしている

インポートされた翼型は「翼型名のGroupオブジェクト＞Wireオブジェクト」という入れ子構造になっている

これから使うのはWireオブジェクトの方なので，App.ActiveDocument.getObject(FoilName).Group[0]でGroupオブジェクトに入っているWireオブジェクトを取得している

このWireオブジェクトに対して拡大縮小を行って翼断面を作成する

Draft.clone()でWireオブジェクトのCloneを作成し，Scaleで拡大縮小を行う

ここまで実行すると次のようになる

参考
≫FreeCAD 翼形状の作成 - XSim
≫Object name - FreeCAD Documentation
≫Tree view - FreeCAD Documentation
≫FreeCAD Scripting Basics - FreeCAD Documentation
≫Std Group - FreeCAD Documentation
≫Draft Wire - FreeCAD Documentation
≫オブジェクトのコピー、移動、回転
≫FreeCAD 3次元オブジェクトの拡大縮小 - XSim
≫Draft Clone - FreeCAD Documentation
≫Draft Scale - FreeCAD Documentation
≫Placement - FreeCAD Documentation

Cloneの作成

上で作成したCloneをもとに，Part DesignワークベンチのCloneを作成する

### Begin command PartDesign_Clone
App.ActiveDocument.addObject('PartDesign::Body','Body')
App.ActiveDocument.addObject('PartDesign::FeatureBase','Clone')
BodyObject = App.ActiveDocument.getObject('Body')
CloneObject = App.ActiveDocument.getObject('Clone')
CloneObject.BaseFeature = CloneObject2D
CloneObject.Placement = CloneObject2D.Placement
CloneObject.setEditorMode('Placement',0)
BodyObject.Group = [CloneObject]
BodyObject.Tip = CloneObject
CloneObject.ViewObject.Transparency=getattr(CloneObject2D.getLinkedObject(True).ViewObject,'Transparency',CloneObject.ViewObject.Transparency)
CloneObject.ViewObject.DisplayMode=getattr(CloneObject2D.getLinkedObject(True).ViewObject,'DisplayMode',CloneObject.ViewObject.DisplayMode)
### End command PartDesign_Clone

Part Designワークベンチを使うことで，自動的にBodyオブジェクトが作成される

ここまで実行すると次のようになる

参考
≫PartDesign Body - FreeCAD Documentation

Padの作成

上で作成したCloneをPadで10㎜押し出してリブを作成する

### Begin command PartDesign_Pad
Gui.ActiveDocument.ActiveView.setActiveObject('pdbody',BodyObject,'')
BodyObject.newObject('PartDesign::Pad','Pad')
PadObject = App.ActiveDocument.getObject('Pad')
PadObject.Profile = (CloneObject, ['Face1',])
PadObject.Length = 10.0
PadObject.ViewObject.ShapeColor=getattr(CloneObject.getLinkedObject(True).ViewObject,'ShapeColor',PadObject.ViewObject.ShapeColor)
PadObject.ViewObject.LineColor=getattr(CloneObject.getLinkedObject(True).ViewObject,'LineColor',PadObject.ViewObject.LineColor)
PadObject.ViewObject.PointColor=getattr(CloneObject.getLinkedObject(True).ViewObject,'PointColor',PadObject.ViewObject.PointColor)
PadObject.ViewObject.Transparency=getattr(CloneObject.getLinkedObject(True).ViewObject,'Transparency',PadObject.ViewObject.Transparency)
PadObject.ViewObject.DisplayMode=getattr(CloneObject.getLinkedObject(True).ViewObject,'DisplayMode',PadObject.ViewObject.DisplayMode)
Gui.ActiveDocument.setEdit(BodyObject,0,'Pad.')
Gui.Selection.clearSelection()
### End command PartDesign_Pad
App.ActiveDocument.recompute()

PadObject.Length = 10.0で押し出し距離を指定している

ここまで実行すると次のようになる

参考
≫PartDesign Pad - FreeCAD Documentation

Pocketの作成

上で作成したPadオブジェクトにPocketを使って穴をあける

### Begin command PartDesign_NewSketch
BodyObject.newObject('Sketcher::SketchObject','Sketch')
SketchObject = App.ActiveDocument.getObject('Sketch')
SketchObject.Support = (PadObject,['Face'+str(len(PadObject.Shape.Faces)),])
SketchObject.MapMode = 'FlatFace'
### End command PartDesign_NewSketch
Gui.runCommand('Sketcher_CompCreateCircle',0)
SketchObject.addGeometry(Part.Circle(App.Vector(x0,y0,0),App.Vector(0,0,1),thickness*0.9/2),False)

### Begin command PartDesign_Pocket
BodyObject.newObject('PartDesign::Pocket','Pocket')
PocketObject =App.ActiveDocument.getObject('Pocket')
PocketObject.Profile = SketchObject
PocketObject.Length = 100.000000
PocketObject.Length2 = 100.000000
PocketObject.Type = 0
PocketObject.UpToFace = None
PocketObject.Reversed = 0
PocketObject.Midplane = 0
PocketObject.Offset = 0
PocketObject.ViewObject.ShapeColor=getattr(PadObject.getLinkedObject(True).ViewObject,'ShapeColor',PocketObject.ViewObject.ShapeColor)
PocketObject.ViewObject.LineColor=getattr(PadObject.getLinkedObject(True).ViewObject,'LineColor',PocketObject.ViewObject.LineColor)
PocketObject.ViewObject.PointColor=getattr(PadObject.getLinkedObject(True).ViewObject,'PointColor',PocketObject.ViewObject.PointColor)
PocketObject.ViewObject.Transparency=getattr(PadObject.getLinkedObject(True).ViewObject,'Transparency',PocketObject.ViewObject.Transparency)
PocketObject.ViewObject.DisplayMode=getattr(PadObject.getLinkedObject(True).ViewObject,'DisplayMode',PocketObject.ViewObject.DisplayMode)
### End command PartDesign_Pocket
App.ActiveDocument.recompute()
Gui.ActiveDocument.resetEdit()

SketchObject.addGeometry(Part.Circle(),False)で穴をあけるためのスケッチを追加している

PocketObject.Length = 100で押し出しの深さでは，十分に貫通する長さを指定している

ここまで実行すると次のようになる

参考
≫Sketcher scripting - FreeCAD Documentation
≫PartDesign Pad - FreeCAD Documentation

.csvの読み込み

圧力分布の.csvファイルを読み込む

### Read .csv
DirPath = os.path.dirname(App.ActiveDocument.FileName)
ReadName, Filter = PySide2.QtWidgets.QFileDialog.getOpenFileName(None, "Read a file", DirPath, "All files (*.*);;") # Open dialog box
with open(ReadName) as f: # Open .csv
    csv = csv.reader(f)
    csv = [line for line in csv]
    csv = csv[1:-2]
csv = [line for line in csv]
csv = [[float(str) for str in line[:2]] for line in csv]
xCp = np.array(csv)[:,0]
Cp = np.array(csv)[:,1]
LE_idx = np.argmin(np.abs(xCp-0.0))

csvファイルの読み込みはcsv.reader()を使う

.datと同じように前縁のインデックス（LE_idx）を取得している

ここまで実行すると次のようになる

参考
≫PythonでCSVファイルを読み込み・書き込み（入力・出力） _ note.nkmk.me

FEMの立ち上げ

FEMワークベンチを立ち上げる（?）

### Begin command FEM_Analysis
ObjectsFem.makeAnalysis(App.ActiveDocument, 'Analysis')
FemGui.setActiveAnalysis(App.ActiveDocument.ActiveObject)
ObjectsFem.makeSolverCalculixCcxTools(App.ActiveDocument)
FemGui.getActiveAnalysis().addObject(App.ActiveDocument.ActiveObject)
### End command FEM_Analysis

ここら辺のコードはあんまり踏み込んで理解していないので無駄があるかもしれない

参考
≫FEM Tutorial Python - FreeCAD Documentation

固定境界条件の設定

先ほどPocketで作成した穴に固定境界条件を適用する

### Begin command FEM_ConstraintFixed
App.ActiveDocument.addObject("Fem::ConstraintFixed","ConstraintFixed")
App.ActiveDocument.ConstraintFixed.Scale = 1
App.ActiveDocument.ConstraintFixed.References = [(App.ActiveDocument.Pocket,"Face"+str(len(PocketObject.Shape.Faces)))]
App.ActiveDocument.Analysis.addObject(App.ActiveDocument.ConstraintFixed)
for amesh in App.ActiveDocument.Objects:
    if "ConstraintFixed" == amesh.Name:
        amesh.ViewObject.Visibility = True
    elif "Mesh" in amesh.TypeId:
        aparttoshow = amesh.Name.replace("_Mesh","")
        for apart in App.ActiveDocument.Objects:
            if aparttoshow == apart.Name:
                apart.ViewObject.Visibility = True
        amesh.ViewObject.Visibility = False
Gui.ActiveDocument.setEdit('ConstraintFixed')
### End command FEM_ConstraintFixed

PocketObject.Shape.FacesでPocketオブジェクトが持つすべての面のリストを入手できる

穴をあけたときにできた面（穴の側面）は最後に追加された面なので，リストの1番最後のインデックスでアクセスできる

scaleが何を設定しているかは分からない，forの中も分からない

ここまで実行すると次のようになる

参考
≫FEM ConstraintFixed - FreeCAD Documentation

圧力条件の設定

リブの側面に圧力条件を設定する

face_idx = 1
for face in PocketObject.Shape.Faces:
    if face.Surface.Axis.z == 0 and not face.Surface.Axis.y == 0:
        if face.Surface.Axis.y<0: # Upper surface
            idx = np.argmin(np.abs(xCp[:LE_idx]*chord-face.CenterOfMass.x))
        else: # Lower surface
            idx = LE_idx+np.argmin(np.abs(xCp[LE_idx:]*chord-face.CenterOfMass.x))
        
        ### Begin command FEM_ConstraintPressure
        ObjectName = "ConstraintPressure"+str(face_idx).zfill(3)
        App.ActiveDocument.addObject("Fem::ConstraintPressure",ObjectName)
        ConstraintPressureObject = App.ActiveDocument.getObject(ObjectName)
        ConstraintPressureObject.Pressure = abs(Cp[idx]) # MPa
        ConstraintPressureObject.Reversed = bool(Cp[idx]<0) # Positive/Negative pressure
        ConstraintPressureObject.Scale = 1
        ConstraintPressureObject.References = [(App.ActiveDocument.Pocket,"Face"+str(face_idx))]
        App.ActiveDocument.Analysis.addObject(ConstraintPressureObject)
        for amesh in App.ActiveDocument.Objects:
            if ObjectName == amesh.Name:
                amesh.ViewObject.Visibility = True
            elif "Mesh" in amesh.TypeId:
                aparttoshow = amesh.Name.replace("_Mesh","")
                for apart in App.ActiveDocument.Objects:
                    if aparttoshow == apart.Name:
                        apart.ViewObject.Visibility = True
                amesh.ViewObject.Visibility = False

        print("Face"+str(face_idx),Cp[idx],face.Surface.Axis.y > 0)
        ### End command FEM_ConstraintPressure
    
    face_idx += 1

PocketObject.Shape.FacesでPocketオブジェクトが持つすべての面のリストを入手できるが，リブの側面のみのリストが欲しいので面の法線ベクトルで区別する

face.Surface.Axisで各面の法線ベクトルが取得できる

face.Surface.Axis.z != 0の面がリブの表裏面で，face.Surface.Axis.y == 0の面が後縁の面である

面の重心位置の座標はface.CenterOfMassで取得できる

.csvファイルのデータからface.CenterOfMass.xと最も近い位置のデータのインデックスを取得する

Pressure = abs(Cp[idx])で圧力の絶対値（単位はおそらく[MPa]）を設定する

Reversed = bool(Cp[idx]<0)で圧力の向き（正圧/負圧）を設定している

ここまで実行すると次のようになる