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

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

スポンサーリンク

概要

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

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

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

FreeCAD Documentation

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

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

FreeCAD Windows でのインストール方法 - XSim
インストーラーを以下のサイトからダウンロードします……
FreeCAD: Select your platform
FreeCAD, the open source 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ファイルができる

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
 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解析を行う

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

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

ソースコードはこちら

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

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

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

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

左下に表示されているのが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

FreeCADFreeCADGuiはPython Consoleに合わせてそれぞれAppGuiとしてインポートしている

.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座標を持つ下面の点を見つけ,最大翼厚を計算している

x0y0は後で翼にあける穴の中心座標である

参考
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オブジェクトを取得している

左上のTree Viewに注目

この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オブジェクトが作成される

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

Tree Viewに注目

参考
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)で圧力の向き(正圧/負圧)を設定している

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

見えにくいがリブ側面に圧力条件が設定されている

参考
FEM ConstraintPressure – FreeCAD Documentation

これでプログラムの説明は完了である

使い方

それではこのプログラムを使って実際に翼型のFEM解析を行ってみる

本当はMaterialを発泡スチロールにして機対速度から計算した圧力を作用させたかったが,FreeCADだと柔らかすぎる材料の計算はできない(?)らしい

Various CalculiX errors - FreeCAD Forum

そこで今回は適当な材料と適当な圧力の単位を使って定性的な解析を行うことにする

ソースコードの置き場所

上で説明したソースコードを置くべきディレクトリを調べる

「Macro>Macros …」をクリックする

「Execute macro」のウィンドウが開くので,下の「User macros location」のパスを調べる

筆者の場合は「C:/Users/XXX/AppData/Roaming/FreeCAD/Macro」だったので,ここに作成した.pyファイルを置けばいい

マクロの実行は右上の「Execute」を実行すればいい

プログラムの実行

「Macro>macros …」でウィンドウを開き,.pyを選択して「Execute」をクリックする

ウィンドウが開くので「NACA4412.dat」を選択して「開く」をクリックする

続いてウィンドウが開くので「Cp_Graph.csv」を選択して「開く」をクリックする

翼型がインポートされ,固定境界条件と圧力条件が設定される

それ以外の条件の設定

解析を行うために必要な条件の設定を行う

FEMワークベンチに切り替える

左上のアイコンから「Creates a FEM material for solid」を選ぶ(黄色い球体のアイコン)

「Material card」から「ABS-Generic」を選択して「OK」をクリックする

「Body」を選択した状態で中央右上の「Create a FEM mesh from a shape by Gmesh mesher」をクリックする

「Max element size/Min element size」を適当に入力して「Apply」をクリックする

メッシュが作成されるので,納得したら「OK」をクリックする 

解析して結果を見る

Tree Viewにある「SolverCcxTools」をダブルクリックする

「Write .inp file」をクリックして「Run CalculiX」をクリックする

計算が終わったら「Close」をクリックする

Tree Viewに「CCX_Results」が追加されるのでダブルクリックする

「Max Shear Stress」を選択して「Show」にチェックを入れて「Factor」を1.0に設定すれば下の画を表示することができる

おわりに

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

結果があっているかどうかの判断には正直あまり自信がないが,それっぽい結果は得られたと思う

より正確に解析したい人は発泡スチロールの特性を調べてみてもいいかもしれない

EPS断熱建材の力学的特性
曲げ強さ、圧縮強さ、ひずみ率と荷重、圧縮クリープ、圧縮弾性率、ポアソン比、線膨張係数など力学的性質解説


↓関連記事

FreeCAD
「FreeCAD」の記事一覧です。

↓おすすめ記事

コメント

タイトルとURLをコピーしました