Rhinocerosが構造解析ツールになる!?

はじめに

この記事はAEC and Related Tech Advent Calendar 2020の17日目の記事です。

adventar.org

本題は２．Rhinoceros+Openseesの導入ですので、１．自己紹介はお暇であれば目を通して頂ければと思います。

１．自己紹介

私は、現在、某アトリエ構造設計事務所で構造設計業務をしています。他のアドベントカレンダーに投稿されてる方々のようにデジタル系の仕事をしているわけではないですが、この１，２年で業務にデジタル系のテクノロジーを少しずつ組み込んできました。

具体的には・・・

①複雑な計算式をpythonで解いて業務に役立つプログラムを作成

構造計算では様々な計算式を使用しており、業務効率化のためにプログラムを書くことがあります。
例えば、木造の垂壁＋腰壁の耐力式は以下に示すように複数の連立方程式から構成されており解くのが面倒です。。。（以下の式を解いてエクセルで数値を出すことも可能ですが、やりたくない。。）

　　　　　　　　　　　　　　　　　出典：木造住宅の耐震診断と補強方法

そこで、pythonというプログラミング言語が持つ豊富な数値計算ライブラリーを利用して、以下のコードを書くことで連立方程式を解くことができます。

#7連立方程式を解く
from sympy import *

de1=Symbol("de1")
de2=Symbol("de2")
de3=Symbol("de3")
Q1=Symbol("Q1")
Q2=Symbol("Q2")
P1=Symbol("P1")
P2=Symbol("P2")

expe1=Q1-P1-P
expe2=Q2-P2-P
expe3=Q2*h2+P1*h0-Q1*(h0-h1)
expe4=de1-de2-Q1*h1/K1
expe5=de3-Q2*h2/K2
expe6=de2-(h0-h1)*de1/h0-h1**2*(h0-h1)**2/(3*E*I*h0)*Q1+h1*h2*(h0**2-h1**2-h2**2)/(6*E*I*h0)*Q2
expe7=de3-h2/h0*de1+h2**2*(h0-h2)**2/(3*E*I*h0)*Q2-h1*h2*(h0**2-h1**2-h2**2)/(6*E*I*h0)*Q1

sol=solve([expe1,expe2,expe3,expe4,expe5,expe6,expe7])
di1=sol[de1]
di2=sol[de2]
di3=sol[de3]

sympyというライブラリーを用いれば、上記のコードで7元連立方程式を解くことができます。煩雑な計算は計算機に投げて、その結果がおかしくないか設計者が判断することで人間と計算機で分業を行っています。

②Rhinoceros+Grasshopperで意匠検討（え、構造設計の仕事じゃな、、）

また、最近はお付き合いのある意匠事務所でも３DCADのRhinocerosを導入されるところが出てきており、私がGrasshopperで意匠検討用のコンポーネントを開発するということもあります。

例えば、屋根の軒天をみせるプロジェクトでは
垂木のボリュームを確認するために、スライダーをいじることで垂木のピッチやサイズをリアルタイムに感覚的に変化させることができるコンポーネントを開発して構造だけではなく、意匠の設計支援も行っています。

↓任意のサーフェイスをつくって、そのサーフェイス上に垂木を並べる

２．Rhinoceros+Openseesの導入

構造解析のソルバーにはOpenseesを採用します。
Openseesとはカリフォルニア大学バークレー校(UCB)がオープンソースとして公開している有限要素解析フレームワークのことです。
ソルバーをオープンソースとして公開するなんて太っ腹ですね！

ただひとつ問題点があるのが、入力モデルの作成が煩雑な点と結果が数値でしか出力されずビジュアル的に分かりにくい点です。

そこで最近、国内外で先ほど紹介したRhinocerosとOpenseesを組み合わせたソフトウェアが開発されております。

以下に紹介します。

①OpenSees for Grasshopper

北九州市立大学の藤田先生はRhino+GrasshopperでOpenseesが使えるようにGrasshopperのコンポーネントを開発されています。
設計実務でも使いやすいように面荷重の入力など、設計者目線の開発をされています。
ちなみに、こちらは開発中とのことですが、フリーでダウンロードできます。興味のある方は以下のnotionを参照して頂ければと思います。

Notion – The all-in-one workspace for your notes, tasks, wikis, and databases.

②Alpaca4d

Marco Pellegrino氏とDomenico Gaudioso氏が開発されたGrasshopper のコンポーネントです。food4Rhinoからダウンロードできます。

Alpaca4d

個人的に物凄く興味があるのが、これがpythonで実装されたOpenseespyをソルバーに使っている点です。このOpenseespyはpython3系でしか動かないのですが、Rhinocero及びGrasshopperはIronpython2.7というpython2系をインタプリターとして採用しており、本来であればOpenseespyはRhinoceros上では動かすことさえできないはずなのです。

構成ファイルを見てみると、Alpaca4d内にpython3の実行ファイルがあるので恐らく、そちらにOpenseesの計算を実行させているのだと思うのですが、どのように行っているのか皆目見当がつかず、世界の広さを感じました。。
（もし、分かる方がいたら、教えてください。。。）

③compas_fea

ETHZ（スイス連邦工科大学チューリッヒ校）で開発されているCompasというデジタルファブリケーションのためのフレームワークのひとつに構造解析を行うcompas_feaというパッケージがあります。導入方法は以下のリンク先を参照してください。

compas_fea

こちらは上記の２つとは異なり、Grasshopperのコンポーネントではなくpythonのライブラリーです。

構造解析を実行するまでの流れは、
１．Rhinocerosもしくはblenderでモデル作成（Rhinoceros or blender）
２．モデルを元にopenseesが読み込める入力ファイルを作成(compas_fea)
３．入力ファイルをopenseesの実行ファイルに読み込ませる(compas_fea)
４．構造解析の実行(opensees)
５．数値結果を読み込んでRhinoceros及びblender上に描画(compas_fea)

compas_feaの個人的お気に入りポイントはpythonで操作できる点とGrasshopperを起動させなくても使える点です。

grasshopperを使用しないrhinoベースのモデリングでrhinopythonを叩くことで構造解析できる。
こちら側でUIを作りこめば、普通に実務で使えそう。
ちなみにソルバーはopenseesです。#opensees pic.twitter.com/q1qHzlZIRn

— 🦅HiguchiYUSUKE🐍 (@8U_H8) November 22, 2020

以下にソースコードとrhinoのモデルを置いておきます。（Git hub整備中につき申し訳ありませんが。。）

下のコードのpath='   '内は任意に変えてください。

from compas_fea.cad import rhino
from compas_fea.structure import ElementProperties as Properties
from compas_fea.structure import GeneralStep
from compas_fea.structure import GravityLoad
from compas_fea.structure import PinnedDisplacement
from compas_fea.structure import PointLoad
from compas_fea.structure import Steel
from compas_fea.structure import Structure
from compas_fea.structure import TrussSection

p  = 7850
A1 = 0.0008
A2 = 0.0005
A3 = 0.0001

# Structure

mdl = Structure(name='truss_frame', path='C:/Users/User/rhinopython/')

# Elements

rhino.add_nodes_elements_from_layers(mdl, line_type='TrussElement', layers='elset_main', pL=A1*p)
rhino.add_nodes_elements_from_layers(mdl, line_type='TrussElement', layers='elset_diag', pL=A2*p)
rhino.add_nodes_elements_from_layers(mdl, line_type='TrussElement', layers='elset_stays', pL=A3*p)

# Sets

rhino.add_sets_from_layers(mdl, layers=['nset_pins', 'nset_load_v', 'nset_load_h'])

# Materials

mdl.add(Steel(name='mat_steel', fy=355, p=p))

# Sections

mdl.add([
   TrussSection(name='sec_main', A=A1),
   TrussSection(name='sec_diag', A=A2),
   TrussSection(name='sec_stays', A=A3),
])

# Properties

mdl.add([
   Properties(name='ep_main', material='mat_steel', section='sec_main', elset='elset_main'),
   Properties(name='ep_diag', material='mat_steel', section='sec_diag', elset='elset_diag'),
   Properties(name='ep_stays', material='mat_steel', section='sec_stays', elset='elset_stays'),
])

# Displacements

mdl.add(PinnedDisplacement(name='disp_pinned', nodes='nset_pins'))

# Loads

mdl.add([
   PointLoad(name='load_v', nodes='nset_load_v', z=-15500),
   PointLoad(name='load_h', nodes='nset_load_h', x=5000),
   GravityLoad(name='load_gravity', elements=['elset_diag', 'elset_main']),
])

# Steps

mdl.add([
   GeneralStep(name='step_bc', displacements=['disp_pinned']),
   GeneralStep(name='step_loads', loads=['load_v', 'load_h', 'load_gravity'], factor=1.5, increments=300),
])
mdl.steps_order = ['step_bc', 'step_loads']

# Summary

mdl.summary()

# Run

mdl.analyse_and_extract(software='opensees', fields=['u', 's', 'sf', 'cf', 'rf'], ndof=3)

rhino.plot_data(mdl, step='step_loads', field='um', radius=0.1, scale=10, cbar_size=0.3)
rhino.plot_data(mdl, step='step_loads', field='sf1', radius=0.1, cbar_size=0.3)  # abaqus:sxx opensees:sf1
rhino.plot_reaction_forces(mdl, step='step_loads', scale=0.05)
rhino.plot_concentrated_forces(mdl, step='step_loads', scale=0.2)

print(mdl.get_nodal_results(step='step_loads', field='um', nodes='nset_load_v'))
print(mdl.get_nodal_results(step='step_loads', field='rfm', nodes='nset_pins'))
print(mdl.get_element_results(step='step_loads', field='sxx', elements='elset_main'))

ライノのモデルはこちら

truss__frame.3dm

509183 Bytes

ファイルダウンロードについて

ダウンロード

3．今後の展望

今まで、人の手作業では実現できなかった複雑な形状も
最近では3DCADやCGソフトでモデリングできるようになりました。
ただ、現実世界にその造形物を出力させるためには重力やその他の外乱に対して問題ないか検証する必要があります。そのためには、CADやCGモデルに連動した構造解析機能が必須だと考えており、そのための開発を引き続き行っていきたいと考えています。使っていて楽しいUIを考えるのも楽しそうです。
また、同時に実際にどのようにつくるかも、工法の点でも考えていく必要があります。

そこで先ほど紹介したデジタルファブリケーションのフレームワークであるCompasには注目しています。

このCompasの存在は石津さんや堀川さんらが主催されている勉強会TokyoAECIndustryDevGroupで教えていただきました。構造解析機能のcompas_feaだけでなくロボットと連動する機能も公開しております。

COMPAS - Robots

こちらもRhinoceros上で動かすことができ、ROSとも連動しているため実機があれば動かすことも可能なようです。（実機で動かしてみたい。。）

Grasshopper上でロボットアームを任意の曲線パスに合わせて動かす。まだ、謎な部分があるがとりあえず動きました。ROSをつかっているため、KUKAに関わらず汎用的にどのロボットにもつかえる（はず）https://t.co/RY1vreJ9by #TokyoAECIndustryDevGroup　#grasshopper pic.twitter.com/PpNRrOO3OC

— 🦅HiguchiYUSUKE🐍 (@8U_H8) July 31, 2020

近い将来、CADで作ったモデルを構造解析し、ロボットで製作できる未来がくるかもしれません！（それができる人になりたい）

そのために取り合えず、年末はGit hubの整備とホームページの製作頑張ります（笑）

最後ですが、もしよければ私のtwitterも登録してどんなことしているのか見て頂ければと。。

蛇足

プログラミング学習について

元々、私はfortranというプログラミング言語を大学時代に研究でつかっていたため、プログラミングには苦手意識はありませんでした。ただ、fortranは研究分野では、まだ現役で使用されているようですが、実務で使うには少し不便（計算速度は速いが高級関数がなくて、すべて自分で実装する必要があります）と考えていたので、ライブラリーが豊富でユーザー数の多いpythonというプログラミング言語を就職後に少しずつ学習し始めました。

pythonの利点は、無料で使えるライブラリーが豊富で本当に有料じゃないの！？と思われるものまで無料で使える点です。

例えばmathmaticaと呼ばれる計算ソフトは個人用で購入しても２～１０万円近くしますが、pythonだと同様の機能を持つライブラリーが無料で使用できます。
そのライブラリーの導入方法を以前、noteに書いたので興味のある方はご覧ください。

私自身、学生時代にプログラミングを学習して感じていたことが、
全てのプログラムを０から実装するのは大変で、車輪の再発明ばかりをしている気がしていました。
そのため、私にとっては既に実装されている豊富なライブラリーを持つプログラミング言語（開発者が多い言語は今後さらに、そのライブラリーが増える）は魅力的で、自分のやるべきことだけに集中できます。（ただでさえ、建築設計は長時間労働のため全てを実装する時間を確保するのは難しいですよね。。）

また、pythonはCADやCGソフトと親和性が高く、最近、建築の世界でも使われて始めているRhinoceros(３DCAD)やRevit(BIMソフト)、Houdini(最近話題？のCGソフト)などでも使えるので、お勧めの言語です。

独学で学習しているが、刺激が欲しいという方は、
東京大学の松尾研究室が主催されているGCIのオンライン講座(東京大学グローバル消費インテリジェンス寄付講座)をお勧めします。私は昨年受講しましたが、知らない関数やライブラリーだらけで大変勉強になりました。最終的に機械学習まで学ぶことになります。

東京大学グローバル消費インテリジェンス寄附講座

こちらの講座は東京大学の学生でなくても、簡単なプログラミングの試験をパスすれば無料で受講できるようになります。内容は教材が週一で送られてきて、それを自学自習しながら与えられた課題をクリアしていくというものです。今年は既に募集を終えているようなので、興味のある方は来年是非参加してみてください。

プログラミングは、自学自習できるようになっておくとよいと思います。

それでは、よき計算生活を