工学解析って何？！　日本オンリーワン株式会社

イメージ（右脳）

人間の脳は、右脳と左脳が脳幹でつながった構造をしている。一般の俗説では、左脳は計算とか理論的思考とかをつかさどり、右脳はイメージとか芸術性・創造性をつかさどると信じられている。

レーシングドライバーで事故により左脳を損傷した人が、絵画のすばらしい才能を示したりすることがあるのは、左脳が機能しないため右脳をフルに使っているからであろう。

しかし、科学技術ではコツコツ理論で計算したとか、地道に１つ１つ洗い出す作業をして到達した理論であるとかが評価される場合が多いよう考えられる。

では、右脳で考えられたものと言うのはいったい何を技術的に意味するのであろうか？

可能性へのイメージと私は考える。

完璧ではなく、ボヤッとしたイメージの合成や思いつきは、理論を越えて存在するもので１つの可能性という右脳の創造物であろう。

１０００億メッシュ、計算時間１０００年という膨大な計算パワーとリソースそして時間を使わないと答えが出ない命題に、一瞬でその結果のイメージを合成するのは右脳のイメージ力である。

ここで左脳で体験・経験した事柄を右脳との連動により一瞬でイメージに合成することが科学技術では要求され、さらにそれを左脳で分析・検証しなくてはならない。

今、右脳のパワーに関する書籍などがたくさん出版されているが、科学技術においては再度述べるが左脳的思考が優先されるようである。

可能性をさっとイメージする力は、脅威の力であり、このことを科学技術の予測の分野に取り入れて解答に迅速ににアプローチすることが今後の科学技術にイノベーションを起こす原動力となるであろう。

CYGWIN

WindowsでUNIX環境を使いたいときには、フリーでCygwinという手がある。
Cygwinのホームページ
ここで開発環境を紹介する。
Fortran　F77のプログラム開発環境をしめす。
ここでメインプログラム１個、サブルーチン４個の単純プログラムと仮定する。

main:値の入力（設定）とサブルーチンコール
s1 :受け取った２つの値の合計出力
s2 :受け取った２つの値の差を出力
s3 :受け取った２つの値の積の出力
s4 :受け取った２つの値の除算を出力

ここで
>make f
とすると全てのプログラムをコンパパイルリンクしfと言う実行モジュールを作るようMakefileを作った。
プログラムを改変すると改変したプログラムのみコンパイルしなおし全体とリンクし実行モジュールを作り直すものである。

Makefile:
#
# Makefile for Fortran
#
PROGRAM=f
F77=f77
FLAG=-c -g
LINK=f77
LDFLAG=-lm -g

.SUFFIXES=:.o.f

.f.o:
$(F77) $(FLAG) $<

CFILE=main.f\
s1.f s2.f s3.f\
s4.f

F77OBJ=$(CFILE:.f=.o)

$(PROGRAM):$(F77OBJ) $(LINK) $(F77OBJ) -o $(PROGRAM) $(LDFLAG)


コンパイラーf77は、開発環境に揃えてgcc,f90,f95などと置き換えてもよい。
CFILEはコンパイルするプログラムの一覧でブランク区切りで書いてゆき、継続行はバックスラッシュで継続する。
一連のシステムをダウンロードする
初心者の方使ってみてください。

科学者になりたい！

小学生時代いやもっと前から科学的事柄に興味があった。生物いや虫の観察・収集から草花の観察などを幼少のとき体験してきた。図鑑ばっかりみてる少年であった。小学生時代の成績は、後ろから３番目でまったく出来ない生徒であった。

中学に入って、初めてこういった知見を認めてくれたのが中根先生であった。草花の観察・写生、太陽の観察、理科の実験などで評価していただいた。中根先生は、私が早稲田実業１年生のときに亡くなられたが、葬式に参列し顕花の時に、「科学者になりたい！」と心の中で祈った。

現在、いろいろ紆余曲折を踏んで東大生産研やJEITAのスーパーコンピューティング部会などと関係もできて、科学技術者の端くれになりつつある。

また幼少の頃から富士神社境内で遊び、そこの祭りが大の楽しみであった。今年も６月３０日～７月２日まで富士神社の祭りである。最近の世相を反映してかDSなどの店もあるがやはり、金魚すくいや綿飴は定番である。

あの時誓ったことが現実のものとなりつつある今、私は先生に恵まれたと思う。今年もまた富士神社の祭りが行われている。

キャビティフローでの騒音発生問題についての考察

図のようなキャビティと気体の流れが発生させる騒音について以下のように考察する。

問題を解析で解く事と実験的に気体粒子の軌跡を考えてみると想像ではあるが以下の図のようになると考えられる。

気体は、キャビティのエッジからバックステップフローによるウエイクを生じ、流線の巻き込み現象が起こり、一部はキャビティ内で旋回する。この時に生じる圧力変動は、線形音響理論から次の2つの周波数の式で表される値に近似される。

騒音もこの圧力変動によって生じる気体の混合から生じることが予想される。

よって周波数の応答は、圧力変動の周波数応答に準じるものと成るであろう。また微小な渦の生成・崩壊による高周波成分が観測されるであろう。

東京大学生産技術研究所で行われている新幹線車両つなぎ目で生じる騒音の解析について、このことを提案したい。

OpenCascade

OpenCascadeのホームページ
Salomeプロジェクトのホームページ

OpenCascadeとSALOMEは、フランスのOpenCascade　Technology社のリリースするGNUオープンソースの構造・流体解析のツールである。オープンソースであり高機能であるため、最低限の投資でトップパフォーマンスをたたき出す。

OpenCascadeTechnology社は、多くのサポートベンダーとカスタマーを抱え、日本では矢崎総業やMITSUTOYOがサポートベンダーであり開発サポート体制が敷かれている。

商用ライセンスを含み、OpenCascadeとSALOMEは、全てオープンでフリーで誰でも入手可能であり5万ダウンロードの実績がある。

OpenCascadeは、３Dモデリングと可視化部分のPRE・POSTでC++のライブラリーに加え、多くのフリーソフトたとえばVTK、Qｔなどで構成されている。３Dモデリング機能は、IGES/STEPのデータ読み込み機能を備え、UG社のCADなどの対応している。

マルチフィジカル・マルチスケール・連成シミュレーションや統合されたシミュレーションを実現する。 またSALOMEは、OpenCascade上で動作するシミュレーション統合プラットフォームであり、IGES/STEPの形状データ読み込み後、テトラ・プリズム・ヘキサ・混合のメッシュ生成や物理データ貼り付け、境界条件の貼り付けを行い、解析計算と可視化までサポートしている。CADと解析をドラマティックに融合リンクするツールである。

OpenCascade Technology社カスタマーサポートマネージャSergey ZertchaninovにパーミッションをいただいたのでOpenCascade本体の画像を載せる。

Maxima

GNUのフリーソフトMaximaを紹介する。
Maximaのホームページ

Maximaは、MITのMacsymaSystemをWilliam　F.SchelterがCommonLispに実装したものです。

(%i3) solve(x*x*x+1=0,x);
(%o3) [x=-(sqrt(3)*%i-1)/2,x=(sqrt(3)*%i+1)/2,x=-1]

など簡単に数式を処理いたします。
(%i3)の入力で(%o3)を返します。


図形描画機能もあります。
(%i4)plot3d(x^4+y^4,[x,-1,1],[y,-1,1]);

使ってみてください。

OPENFOAM

OpenFoamホームページ

OpenFoamは、ロンドン大学インペリアルカレッジ・ゴスマン研究室に在籍していたHenry Weller(現在英国：OpenCFD Ltd.代表)らによって開発が進められ2004年にGNUのオープンソースとしてリリースされている流体解析のツールである。

Javaで駆動されるPre部分とC++ライブラリーによる各種ソルバー(PDEソルバー)とユーザーアプリの部分から構成されており、解析結果の可視化部分は、別のオープンソースのParaviewで行う様にデータ変換のツールが用意されている。Javaで駆動されるPre部分は管理画面のケースブラウザとメッシャーツール（ブロックメッシャー）とからなる。

ソルバーは、PDEソルバーといい、非圧縮性・圧縮性、乱流（多数のLESモデル）、2層流、自由表面、燃焼問題を扱える。またユーテリティーには、メッシュ操作（並列計算用分割、ダイナミックメッシュ）、LESモデル用POST、渦度、揚力、抗力算出が可能であり、ジョブ管理からデバックレベル管理などの管理画面が用意され、並列処理（LAM/MPI）も可能である。

海外では勿論のことで、オープンソースのCFDツール・高機能であるため日本でも学会等が取り上げ始めている。日本機械学会・日本建築学会・日本音響学会などで、OpenFoamを用いて解析がなされ論文が出ている。

Wikki.co.ukのHrvoje JasakにパーミッションをいただいたのでUSNavyのLarge Eddy Simulationの解析例を載せる。

ローレンツ方程式(その２）

前回に、余計な項が含まれていたので、解きなおした。

>function xdot=f(x,t)
xdot(1)=-10*x(1)+10*x(2)
xdot(2)=-x(1)*x(3)+28*x(1)-x(2)
xdot(3)=x(1)*x(2)-8/3*x(3)
endfunction
>x0=[1;1;1];
>t=linspace(0,500,1000);
>x=lsode("f",x0.t);
>plot(t,x);

Oh,No!

アメリカの金星探査ロケットで起こったことである。それまで順調に飛行していたのが、軌道修正時に行方不明になった。何回か軌道修正して順調に飛行を続けていたが、金星目前で軌道修正した途端、予定されている軌道をはずれ、地球から所在が把握できなくなったのである。

巨額の投資をして開発した高価なロケットが行方不明となり、徹底的な調査が行われた。

原因が判明した。なんとたった一つのコンマがピリオドに化けていたことがその重大な事態を導いたことが分かる。

DO 100 I=1.10と言う記述が見つかりプログラムがDO100I=1.10と解釈して変数DO100Iに1.10を代入しただけとなっていることが判明した。

Oh,No!

ローレンツ方程式

流体のカオス的ふるまいを表す方程式としてローレンツ方程式がある。1963年Lorenz,E.N.が”Deterministic Nonperiodic Flow",Journal of Atmospheric Sciencesに発表した決定論的非周期流れの方程式である。

dx/dt=-px+py
dy/dt=-xz+rx-y
dz/dt=xy-bz

ローレンツがこの論文で与えたのは、
p=10
r=28
b=8/3
であり、この元に解かれるX,Y,Zの軌跡（アトラクタ）が示された。

この方程式を私も解いてみた。OctaveというGNUのソフトを使い解く。

>function xdot=f(x,t)
xdot=zeros(1,3)
xdot(1)=-10*x(1)+10*x(2)
xdot(2)=-x(1)*x(3)+28*x(1)-x(2)
xdot(3)=x(1)*x(2)-8/3*x(3)
endfunction
>x0=[1;1;1];
>t=linspace(0,500,1000);
>x=lsode("f",x0.t);
>plot(t,x);

以上の設定で解けるはずである。

アトラクタは示せないが、確かにカオス的解が得られる。このことが流体のふるまいを示していると言うのいである。アトラクタを示せばよりわかり易いであろうが、たとえば小川の中の石の周りの流れがこのことを表している。同じようなパターンの流れを示したかと思えば異なるパターンにぱっと変わり、また元に戻りこれが非周期的に起こる。川の中の石の周りの流れをじっと見て余生を過ごしたレオナルド・ダ・ビンチもこういったことが連立常微分方程式の初期値問題でこんなことが表せると気づいたらさぞかし喜んだであろう。

地球シミュレータとTSUBAME

2002年3月に登場したES（Earth Simulator)は、ベクトルプロセッサを5,120個並列につなぐことによって40TFLOPSの性能を示し当時コンピュートニックショックと言うほど世界を圧倒してきた。しかしESの登場によって触発されたアメリカは、巻き返しのため躍起で高性能のスパコンを開発し始める。今やBlueJeanという米国IBMのスパコンがESの10倍ほど高速で世界一の高速演算を誇る。

現在日本で一番高速演算するスパコンは、東工大のTSUBAMEであり、世界8位である。X86のCPUを並列につなげて高速演算させるのでX86対応のソフトなら何でも走る。また省電力化にも工夫がなされ空調に工夫がある。ダウンフローとサイドフローの空調で温度調節がうまくなされ、ESよりかなり省電力で運転できる。
　　　<

職人と日本が捨ててはならぬもの

大阪下町の職人。
東京大田区の町工場の職人。
に代表される人たちは、こう思うであろう。
”今俺が作ってやるから、見てろ！！！”

時には図面を見ずに、時には地べたに蝋で図面を書いて作り始める。私は、東京電力㈱や中部電力㈱の発電所サイトで登録され働く職人を見てきた。図面を完全に頭に描いて作っている。手間が早い。段取りに時間をかけ、段取り７割りというまで入念に仕事の準備をする。仕事始めから終わりまでをきっちりこなし養生して次の仕事につなげる。

挨拶は、”ご安全に！”である。
無精ひげを生やしている人もいたり、首に手ぬぐいを巻いている人もいて、地べたに座って昼食をとったりする。新日鉄㈱設計会社の東海設計㈱の工務部の職人は、首に手ぬぐいを巻いたりしてSUN端末のCADで図面を書いたりしている。

今や職人たちがCAD図面や官庁申請書を書いている時代であり、そうできないと工事できない時代である。

元々勘が良く、運動神経抜群なのであるから物覚えも速いし、物づくりを毎日やっているのであるからCADなんか簡単に操作する。出来ないわけがない。

職人は、”俺は学校が嫌いだった”といい、今では人生の中で職聖を見出し立派に毎日働いている。こういった人の中に超ベテランがボウシンとして存在して、勤めている。

この人たちの物づくりは、経験から来る理屈と自信に満ちており筋が通っている。時には理論解や解析解を凌駕するほど精度が高い。

日本は、第二次世界大戦の時、零式戦闘機や戦艦大和や１万メートル射程の大砲を作り出した脅威の物づくり大国である。戦争兵器だけでなく今日は職人技によって開発されたものに溢れている。ありとあらゆる建物は全部職人が作ったものである。

これらのことを職人に感謝して、日本が絶対に捨ててはならぬものとして継承してゆきたい。このような分野で解析や分析の仕事をしてゆこうと思う。

根気

解析の仕事。
何かアカデミックで花形な仕事と捕らえる人が多い。

しかし実際は、泥臭い体力勝負のところが多い仕事である。データを一つ一つ洗い出してゆくような根気が必要である。総研とかリサーチ会社に入ってみたもののその泥臭い・つらい仕事で幻滅してやめる人もいる。ツールを使って簡単に答えが出るものとそうではなくややこしい問題が多くある。

解析の仕事は、一つずつデータを洗い出し、的を絞って行き、最終ターゲットである現象の解明を目的とする根気の要る仕事であって、データ変換とかスクリプト生成とかその都度用があれば作りながら進めることも要求される。

根気こそ、解析の仕事の原動力である。

想像力＆ひらめき

解析の仕事をしていて、「パッ」と閃くことがある。

結果見て、今までの常識で読み解こうとする。力学の常識・実験での経験知により類推する。

本当にその常識どおりか？！！！

流体の現象は、もう古典の力学と呼ばれるほど研究し尽くされたが、まだ非定常・非線形の流体現象に研究の余地がある。材料の応力解析も同じで線形の静弾性応力解析は古典的と呼ばれる分野に入る。

どのような現象の解析でも理論と実験の常識から発想しなければならない。

想像力や閃きも、これら今まで培った理論や実験の常識の中で発想されるべきである。

面白い曲面（その１５）

Sphere9
X = cos(u)*cos(v)*sin(u)
Y = cos(u)*sin(v)*cos(u)*sin(v)
Z = sin(u)*sin(v)*sin(u)*sin(v)
With U[ -pi/2, 0] And V[ 0, pi]



Sphere10
X = cos(u)*cos(v)*sin(u)*sin(v)*sin(v)*sin(v)
Y = cos(u)*sin(v)*cos(u)*sin(v)
Z = sin(u)*sin(v)*sin(u)*sin(v)*cos(u)
With U[ -pi/2, 0] And V[ 0, pi]



Star
X = cos(u)*cos(v)*(abs(cos(1*u/4))^0.5 + abs(sin(1*u/4))^0.5)^(-1/0.3)*(abs(cos(5*v/4))^1.7 + abs(sin(5*v/4))^1.7)^(-1/0.1)
Y = cos(u)*sin(v)*(abs(cos(1*u/4))^0.5 + abs(sin(1*u/4))^0.5)^(-1/0.3)*(abs(cos(5*v/4))^1.7 + abs(sin(5*v/4))^1.7)^(-1/0.1)
Z = sin(u)*(abs(cos(1*u/4))^0.5 + abs(sin(1*u/4))^0.5)^(-1/0.3)
With U[ -pi/2, pi/2] And V[ 0, 2*pi]



Star7
X = cos(u)*cos(v)*(abs(cos(7*v/4))^1.7 + abs(sin(7*v/4))^1.7)^(-1/0.2)*(abs(cos(7*u/4))^1.7 + abs(sin(7*u/4))^1.7)^(-1/0.2)
Y = cos(u)*sin(v)*(abs(cos(7*v/4))^1.7 + abs(sin(7*v/4))^1.7)^(-1/0.2)*(abs(cos(7*u/4))^1.7 + abs(sin(7*u/4))^1.7)^(-1/0.2)
Z = sin(u)*(abs(cos(7*u/4))^1.7 + abs(sin(7*u/4))^1.7)^(-1/0.2)
With U[ -pi/2, pi/2] And V[ 0, 2*pi]