【22卒】航空宇宙民の脱航空就活記
はじめに
みなさん初めまして。黄緑ペンギンと申します。本記事は東京大学航空宇宙工学科／専攻 Advent Calendar 2020の19日目の記事となっております。本当はペンギンの流体力学特性について執筆する予定だったのですが、面白い論文が見つからなかったのと勉強時間が取れなかたので就活に関するポエムを書くことにしました。
最初に申し上げておきますが、本記事は「コロナの影響で航空産業が壊滅した結果就活に苦しんでいる航空の学生」という内容ではありません。「コロナとは無関係に文系就職を志して悩んでいる理系大学院生」という内容です。ただのポエムになってしまいますが、それでもよければ読んでいってください。今後就活する後輩や同期の参考になれば幸いです。
なぜ脱航空するのか
私は中学生のころから「宇宙に携わる仕事がしたい」と考え続けてここまで走って来ました。事実、大学2年生の頃は「博士かJAXAか三菱重工か...」みたいなことを漠然と考えていたと思います。ではなぜ昔からの夢を捨てて今就活しているのでしょうか。いくつか理由を挙げたいと思います。
そもそも大学に残れない・残りたくない
小さい頃は本気で大学教授になりたいと思っていましたが、現実はそう簡単ではありません。学部の卒論では、実験しても全然うまくいかない、考察しようにもそもそも何が起こっているのか見当もつかないなど、座学と研究の違いを垣間見ました。また、この前の研究室輪講では助教に「センスないね」と直球で言われるなど、所詮自分はペーパーテストでそこそこの点を取れるだけの人間であって、研究には少しも向いていないんだなと実感しました。
また、この道を選べば、来年のポストもわからないまま安月給で働き続ける人生を送るかもしれません。私の知っている博士の方々はとんでもなく優秀な方ばかりです。（正直就活で出会う自信満々の外コン外銀志望の学生より遥かに頭いいです。）そんな方々と将来ポスト争いしたところで勝てる未来が全く見えません。私にとって「やりがい」とは安定した生活があればこそ追求できるものなのです。
エンジニアに向いてない
これは致命的な問題で、私は工学部なのに"モノ"に全く興味がない人間でした。もう少し具体的に説明すると、例えば「ロケットがどのような原理で飛行しているか」「ロケットエンジンの作動原理は何か」といったことを学ぶのは面白く感じる一方で、「○○トンの推力を出すためにはプロペラ径やシャフト径をいくつにすべきか」や「破壊が起こらないためにはこの部材の厚みはいくつにすべきか」といったことを考えるのは苦痛の極みでした。
実際、私は街中を走っている車の名前はおろか、いまだに飛行機やロケットの名前もロクに言えません。さすがに論文で頻繁にみるものは覚えましたが、”ロケットの名前クイズ”みたいなのがあれば宇宙好きの高校生に負けちゃいますね。
一方周りはと言うと、私とは対極の"モノ"に対する興味も知識も豊富な人達ばかりでした。（飛行機に"顔がたれ目"という表現を使う人々を後にも先にも見たことがありません。）重工に行って飛行機や車を作るのは彼らの方が絶対に適任です。そういうことは得意な人に任せて、自分は他のところで活躍しようと考えたわけです。
宇宙を仕事にする必要はない
色々と書きましたが、私はやっぱり宇宙が好きです。はやぶさ２の帰還には感涙しましたし、SpaceXの一段垂直着陸はいつみても心躍ります。静かな夜にきれいな星空を眺めるのはいつだって素晴らしいことだし、その遥かなる空へ挑み続ける人は皆かっこいいと心から思っています。
しかし、それと働くということはまた別物ではないかとも思うのです。宇宙についてはニュースをみたりGoogle Scholarで最新の論文を漁ることで追いかけることができます。「好きなことを仕事にしたら嫌いになった」とはよく聞く話でもあります。自分が納得して働ける・自分の能力をフルに発揮できる環境というのは宇宙業界の外に沢山あるはずです。そう考えると、航空宇宙にこだわる必要なんて1ミリもないという結論に至りました。
何のために・どういう方向で就活するのか
そもそも就活とは何を目的にするのでしょうか？「内定を得る」のはもちろんですが、ではその企業を選ぶ理由はなんでしょうか？
就活を始めた頃にリ〇ルートの社員さんがこの話をしてくださったのですが、その一つの答えは「幸せな人生を実現するため」というものでした。なんとも抽象的な話ですし、"浅い"という感想を持たれる方もいるかもしれませんが、間違ってはいないはずです。では「幸せ」とは何か。それが個々人の価値観によって異なり、どんな企業に行きたいかの指標になるというわけです。なお、"就活の軸"についてはモチベーショングラフからの自己分析なども活用しながら決めるのが一般的で、私もそれに倣いました。
なんだかんだで私はまず社会をもっと知るべきだという結論に達し、まずは色んな業界の説明会・インターンに参加することにしました。そこで気に入る企業があれば就職するもよし、もし就職というものが想像以上に自分に向いていない場合は大学に残ればいいやくらいの気持ちで気楽に選考を受けていました。
これまでの就活状況
春
就活にむけて動き出したのは2020年の3月、M1の春でした。とはいってもこの時期は企業イベントはほとんどないので、選抜コミュニティーという就活エージェント(?)の選考をいくつか受けていました。選抜コミュニティーは現在3つ所属しているのですが、「入った意味ないな」程度のものもあった上に、とある企業の採用担当の人事の方が「選抜コミュニティーに入っているかとかは全く見ていない」と仰っていたので無理に受ける必要はないと感じました。また、就活仲間のグループを友達同士でいくつか作ってGDの練習や勉強会、ケース面接の練習をしました。
夏
夏は日経を中心に色々な企業の説明会にいったり、インターンに応募したりしました。この時期は自分の中で働くことのイメージを全く形成できておらず、志望している業界もなかったので、ちょっとでも面白そうだと思ったらとりあえず説明会を聞きに行くようにしていました。ですが今振り返ると無意識に切っていた企業がいくつもあったように感じるので、もっと広く見るべきだったかなと思います。
インターンは5社参加しました。選考の通過率は56%、思ったより落ちるなとヘコんでいましたが、今考えるといいほうだったと思います。業界としてはコンサル・金融・インフラを中心に参加しました。一番面白かったのは某新幹線の強い鉄道会社で、インターン参加前は志望度はかなり低かったのですが、「ここで働くの悪くないじゃん」と思うくらいまで志望度が上がりました。逆にコンサル系は肌に合わない企業が多かったので、あまり自分に合ってない業界なのかもしれないと思いました。インターンに参加すると自分の向き不向きが鮮明に見えてくるのでいいですね。
秋
秋は就活のモチベが低下したこと、研究室の活動が忙しかったことなどが重なり、ほとんど就活はしていませんでした。秋インターンは面白そうだと思ったとある政府系金融機関のみ選考を受け、運よく通過することができました。そのインターンはとても学びの多い有意義なものだった上に、業務もやりがいのある面白いものだと感じたため、その企業は現在の第一志望群になっています。
冬（現在）
冬も就活モチベは低下気味。とりあえず出したものはマッキ〇ゼーと某日経コンサルのITソリューション部門で、後者のみインターンへの参加権を得ました。あとは某財閥系の総合商社2社に応募する予定で、それ以外は特に考えてはいないです（もう疲れた）。年内の内定獲得は不可能、たぶん来年6月まで就活は続きます。
就活で失ったもの
ネガティブなサブタイトルですが、就活を通して色々メンタルが擦り減っているので、愚痴も兼ねて何を失ったか記していこうと思います。
時間
就活は何より時間をとられます。長期休み期間ならまだいいですが、最近は選考や面談のために講義を欠席せざるを得ないこともしばしば。奨学金を受給しているので講義には出席したいのですがどうしようもありません。もう今学期の成績を見るのが怖いです。
研究室での信頼
私の所属する研究室はとにかくミーティングが多く、毎日進捗を生むことを求めらる（たとえば毎朝今日何やるかを説明するミーティングがある）環境です。そんななか私だけ頻繁にスーツを着て「今日は○○社の人事と面談します」と宣言しなければならない辛さ。最近はボスの目線も態度も冷ややかになっているような気がしていて、居心地がよろしくないです。ただ、研究室によっては就活禁止のところもあるらしく、それに比べたら全然恵まれているので感謝しなければいけませんね。
人間性
これは失ったというより、自分の人間としての未熟さが露呈したというほうが正しいと思うのですが、周りと自分を比較して妬んだり卑屈になってしまうということがしばしばあります。例えば誰かが内定をもらった時、仲の良い友達ならば素直に祝福できるのですが、顔見知り程度の人の場合は「あの人が内定貰ってるのに俺は...」という気持ちが最初に来てしまいます。これは自分の不徳の致すところではありますが、メンタルを擦り減らす一因になってしまっています。
最後に
サンタさん内定ください。
電気推進ロケットの数値解析手法
本記事は東京大学航空宇宙工学科／専攻 Advent Calendar 201９の24日目の記事です。
はじめに
みなさんはじめまして。学部4年の黄緑ペンギンです。クリスマス・イヴいかがお過ごしでしょうか。私は毎年クリスマスキャロルの頃にはを聴きながらダラダラしています。
さて、私は人工衛星のエンジンに関する研究室に所属していますが、卒論は実験で書いたので実は数値解析は専門外。「そうだ！Advent Calenderで数値解析の話をするって言えば勉強せざるを得ないじゃん！」ということで私自身も勉強しながらこの記事を執筆しました。
電気推進ロケットとは
弊学科の民でも電気推進に馴染みがある人は極一部だと思うので、本題に入る前に「電気推進ロケットとはなんぞや」という話から始めたいと思います。
イオンエンジンという言葉を聞いたことがあるでしょうか？最近有名な「はやぶさ２」にくっついている光ってるヤツですね(図１)。実はこれが人工衛星の推力を生み出す機械で、電気推進ロケットの一種です。
イオンエンジンの他にもホールスラスタやMPDスラスタなど、様々なタイプの電気推進ロケットがありますが、どれも「電気エネルギーをプラズマ粒子の運動エネルギーに変換する」ことで推力を生成しています。その方法は静電加速だったり誘導加速だったりするのですが、長くなるのでここで深追いするのはやめておきましょう。
プラズマとは
突然"プラズマ"というワードがでてきましたが、コレが今回の主役となります。イメージとしては電離した気体だと思ってください。厳密ではありませんが、電気推進に出てくるプラズマ程度ならばこの理解で大丈夫です。細かい定義は教科書に譲ることにしましょう。
"電離した気体"なのでプラズマの中には電子やイオンがふわふわ漂っているわけですね。そこに電場や磁場をかけるとどうなるでしょうか？きっと電場から力を受けて粒子は飛んでいくでしょうし、しかもその運動は磁力線に巻きつくようなラーマ運動を伴っているはずです。このような様々な効果が絡み合うプラズマ物理は非線形でとても複雑なんです。
ここではプラズマ物理を代表する現象を一つとりあげましょう。ある空間に点電荷を置いた状況を考えます。この電荷が真空中に作る電場
は真空中の誘電率を
として
\begin{align}
\phi_0=\frac{q}{4\pi \varepsilon_0 r}
\end{align}
と表せるのでした。では、プラズマ中に点電荷をおいた場合はどうなるでしょうか。
\begin{align}
\nabla^2\phi=\frac{e}{\varepsilon_0}(n_e-n_\infty)-\frac{q}{\varepsilon_0}\delta(r)
\end{align}
とかけます。ただし、は無限遠における電子密度です。真空中では右辺第二項のみが残って
は
に一致します。
電子密度としてボルツマン分布
\begin{align}
n_e=n_\infty \exp\left(\frac{e\phi}{kT_e}\right)
\end{align}
を仮定すると、電子の運動エネルギーが静電的なポテンシャルエネルギーより十分大きい場合は、で
となる解として
\begin{align}
&\phi(r)=\frac{q}{4\pi\epsilon_0 r}\exp\left(-\frac{r}{\lambda_D}\right)\\
&\lambda_D^2=\frac{\varepsilon_0kT_e}{e^2n_e}
\end{align}
を得ます。ただし、は電子温度、
はボルツマン定数です。
と
を見比べると、プラズマ中では電場がほぼ
のところで断ち切られ、それより遠方に作用が及ばないようになっていることがわかります。何が起こっているのかというと、点電荷のまわりにプラズマ中の粒子が集まってきて電荷を打ち消そうとしているのですね(図２)。この効果をデバイ遮蔽と呼びます。
デバイ遮蔽の距離をデバイ長と呼びます。デバイ長はプラズマが準中性(巨視的に電荷が釣り合っている状態)を保つのに必要な長さで、数値解析における格子サイズもデバイ長をとることが多いです。
PIC法
本題である数値解析の話に移りましょう。プラズマの数値解析にはいくつかの手法があるのですが、そのなかで最も代表的だと思われるParticl-In-Cell (PIC) 法についてとりあげます。
PIC法は、粒子を有限の広がりのもった超粒子として考えて電場や磁場と結びつける手法であり、計算コストが重いのが特徴です。例えばイオンエンジンの放電室の２D解析に1~2ヶ月かかる例も存在します。
また、PIC法はプラズマ粒子の衝突を考慮していないので万能ではありません。とあるプラズマの不安定性など、再現できない現象も存在します。しかし、PIC法は有用な解析手法として知られていて、プラズマ数値解析の代表的な手法といえるでしょう。
そこで、本記事はPIC法の仕組みについて、PICの一種であるCHPIC[3]という手法を参考にして解説していこうと思います。
支配方程式
支配方程式はマクスウェル方程式とニュートン-ローレンツ方程式です。書き下すと以下のようになります。
\begin{align}
&\nabla\cdot\boldsymbol{B}\left(t,\boldsymbol{x}\right)=0\tag{1}\\
& \nabla\times\boldsymbol{E}\left(t,\boldsymbol{x}\right)+\frac{\partial \boldsymbol{B}\left(t,\boldsymbol{x}\right)}{\partial t}=0\tag{2}\\
&\nabla\cdot\boldsymbol{D}\left(t,\boldsymbol{x}\right)=\rho\left(t,\boldsymbol{x}\right)\tag{3}\\
&\nabla\times\boldsymbol{H}\left(t,\boldsymbol{x}\right)-\frac{\partial \boldsymbol{D}\left(t,\boldsymbol{x}\right)}{\partial t}=\boldsymbol{j}\left(t,\boldsymbol{x}\right)\tag{4}\\
&\frac{d}{dt}\gamma m \boldsymbol{v}=q\left(\boldsymbol{E}+\boldsymbol{v}\times\boldsymbol{B}\right)\tag{5}
\end{align}
ただし、光速をとして
\begin{align}
&\frac{d}{dt}\boldsymbol{x}=\boldsymbol{v}\tag{6}\\
&\gamma=\sqrt{\frac{1}{1-\left(v/c\right)^2}}\tag{7}
\end{align}
とおいています。はローレンツ因子と呼ばれる量です。
解析フローのスケマティックを図3, 4に示します。電場・磁場から運動方程式により粒子位置を更新し、粒子の移動によって生じた電流と位置更新後の電荷密度から電場・磁場を更新するというのが大まかな流れとなっています。
Yee格子
電磁場の解析では計算領域をYee格子というメッシュで区切るのが一般的です。Yee格子を図5に示します。Yee格子では、格子の辺に沿って線分の中点に電場を、面の中心に垂直方向に磁場をとります。このようにとる理由はマクスウェル方程式を見ると理解できます。
例えば、磁場と電場を同じ格子点上に定義した場合
\begin{align}
\nabla\times\boldsymbol{E}\left(t,\boldsymbol{x}\right)+\frac{\partial \boldsymbol{B}\left(t,\boldsymbol{x}\right)}{\partial t}=0
\end{align}
を計算するのが面倒になってしまいます。Yee格子を用いればこれらの計算に余計な手間がかからないわけですね。また、メッシュサイズは前述したデバイ長をとることにしましょう。
磁場・電場の更新
磁場と電場はそれぞれ(2), (4)式によって更新されます。ここでは単純な差分法を用いることにします。軸対称の場合はOOPIC[4]などの手法が参考になると思います。また、千葉大学の松本先生によるpCANSなど、共役勾配法を用いた陰的解法を利用したものも存在します。
簡易な式変形なので、ここでは結果だけ示すと、次のようになります。
\begin{align}
B_{i}^{t+\Delta t / 2}&=B_{i}^{t-\Delta t / 2}-\Delta t\\
&\quad\times\left\{\frac{E_{k, x_{j}+\Delta x_{j} / 2}^{t}-E_{k, x_{j}-\Delta x_{j} / 2}^{t}}{\Delta x_{j}}\right.\\
&\quad\left.-\frac{E_{j, x_{k}+\Delta x_{k} / 2}^{t}-E_{j, x_{k}-\Delta x_{k} / 2}^{t}}{\Delta x_{k}}\right\}\\
E_{i}^{t}&=E_{i}^{t-\Delta t}+\frac{\Delta t}{\varepsilon \mu}\\
&\quad\times\left\{\frac{B_{k, x_{j}+\Delta x_{j} / 2}^{t-\Delta t / 2}-B_{k, x_{j}-\Delta x_{j} / 2}^{t-\Delta t / 2}}{\Delta x_{j}}\right.\\
&\quad\left.-\frac{B_{j, x_{k}+\Delta x_{k} / 2}^{t-\Delta t / 2}-B_{j, x_{k}-\Delta x_{k} / 2}^{t-\Delta t / 2}}{\Delta x_{k}}\right\}\\
&\quad-\frac{\Delta t J_{i}^{t-\Delta t / 2}}{\varepsilon}
\end{align}
しかし、これだけではマクスウェル方程式のうち2本の式しか使っていないことになります。残りの2本について考慮しましょう。については、さっき導いた式を代入することで
\begin{align}
\nabla\cdot\boldsymbol{B}^{t+\Delta t / 2}=\nabla\cdot\boldsymbol{B}^{t-\Delta t / 2}
\end{align}
を得るので、時間発展の間にも差分化した計算の丸め誤差の範囲におさまります。
一方(3)式は必ずしも満たされているとは限らないため何らかの補正が必要となります。そこで、さっき求めたを
とおきなおして次のようなポアソン方程式を考えます。
\begin{align}
\nabla^2\phi=\nabla\cdot\boldsymbol{E}^{*}-\frac{\rho}{\varepsilon_0}
\end{align}
これを共役勾配法などによって解いてを求め
\begin{align}
\boldsymbol{E}^t=\boldsymbol{E}^*-\nabla\phi
\end{align}
によって補正を行います。
粒子速度および位置の更新
粒子速度は(5)式によって更新します。その式はなる量
を用いて次のように表せます。
\begin{align}
&\frac{\boldsymbol{u}^{t+\Delta t / 2}-\boldsymbol{u}^{t-\Delta t / 2}}{\Delta t}=\frac{q}{m}\biggl(\boldsymbol{E}^{t}\\
&\quad\quad+\frac{\boldsymbol{u}^{t+\Delta t / 2}+\boldsymbol{u}^{t-\Delta t / 2}}{ 2\gamma^{t}} \times \boldsymbol{B}^{t}\biggr)\tag{8}\\
&\gamma^{t}=\frac{\gamma^{t-\Delta t / 2}+\gamma^{t+\Delta t / 2}}{2}\tag{9}
\end{align}
(8)式右辺に未知数が含まれているので、このままでは計算することができません。そこで新たに次の２つの物理量を定義しましょう。
\begin{align}
&\boldsymbol{u}^{-}=\boldsymbol{u}^{t-\Delta t / 2}+\frac{q \Delta t \boldsymbol{E}^{t}}{2 m}\\
&\boldsymbol{u}^{+}=\boldsymbol{u}^{t+\Delta t / 2}-\frac{q \Delta t \boldsymbol{E}^{t}}{2 m}
\end{align}
これらを用いて(5)式を書き直すと
\begin{align}
\frac{\boldsymbol{u}^{+}-\boldsymbol{u}^{-}}{\Delta t}=(\boldsymbol{u}^{+}+\boldsymbol{u}^{-})\times\frac{q\boldsymbol{B}^t }{2\gamma m}\tag{10}
\end{align}
となります。これは、から
に変化する間に
のローレンツ力のみが作用しているとみなせることを意味しています。
また、(10)式の両辺との内積をとると、右辺が０になるので
が成立することがわかります。つまり、
は
をある角度
だけ回転させたベクトルになります。
さて、(10)式をについて解くと次のような式が得られます。
\begin{align}
&\boldsymbol{u}^{+}=\boldsymbol{u}^{-}+\boldsymbol{s}\left(\boldsymbol{u}^{-}+\boldsymbol{u}^{-}\times\boldsymbol{t}^t\right)\\
&\boldsymbol{t}^t=\hat{\boldsymbol{B}}\tan{\left(\frac{q\Delta t}{2\gamma^t m}B^t\right)}\\
&\boldsymbol{s}=\frac{2 \boldsymbol{t}^{t}}{1+\boldsymbol{t}^{t} \cdot \boldsymbol{t}^{t}}
\end{align}
これらを整理すると、粒子の速度は次のステップで更新されることがわかります。
1. 電場による加速
\begin{align}
\boldsymbol{u}^{-}=\boldsymbol{u}^{t-\Delta t / 2}+\frac{q \Delta t \boldsymbol{E}^{t}}{2 m}
\end{align}
2. 磁場による回転
\begin{align}
\boldsymbol{u}^{\prime}&=\boldsymbol{u}^{-}+\boldsymbol{u}^{-} \times \boldsymbol{t}^{t}\\
\boldsymbol{u}^{+}&=\boldsymbol{u}^{-}+\boldsymbol{u}^{\prime} \times \boldsymbol{s}
\end{align}
3. 電場による加速
\begin{align}
\boldsymbol{u}^{t+\Delta t / 2}=\boldsymbol{u}^{+}+\frac{q \Delta t \boldsymbol{E}^{t}}{2 m}
\end{align}
そして最後に
\begin{align}
\frac{\boldsymbol{x}^{t+\Delta t / 2}-\boldsymbol{x}^{t}}{\Delta t}=\frac{\boldsymbol{u}^{t+\Delta t / 2}}{\gamma^{t+\Delta t / 2}}
\end{align}
によって位置の更新をおこないます。このような電場による移流と磁場による回転を分けて考える手法をBuneman-Boris法と言います。Buneman-Boris法における各ベクトルの関係を図６に示します。
これで粒子位置の更新ができた！と言いたいところですが、(10)式の電場・磁場はどうすればいいのでしょうか？先に示したように、電場・磁場はYee格子上でのみ定義されているため、粒子のいる場所での電場・磁場を新たに与えてあげる必要があります。例えば、粒子のまわりにあるYee格子点上の値を線形に重み付けして定義する方法があります。
粒子位置に対して図７のように,
,
を定め
\begin{align}
w_{j}=\frac{\Delta l_{1}}{\Delta x_{1}}\\
w_{k}=\frac{\Delta l_{2}}{\Delta x_{2}}\\
w_{m}=\frac{\Delta l_{3}}{\Delta x_{3}}
\end{align}
とおきます。このとき粒子の位置における電場は、線形で重み付けすると次のように書けます。
\begin{align}
\boldsymbol{E}_{i}=&\boldsymbol{E}_{j, k, m}\left(1-w_{j}\right)\left(1-w_{k}\right)\left(1-w_{m}\right)\\
&+\boldsymbol{E}_{j+1, k, m} w_{j}\left(1-w_{k}\right)\left(1-w_{m}\right)\\
&+\boldsymbol{E}_{j, k+1, m}\left(1-w_{j}\right) w_{k}\left(1-w_{m}\right)\\
&+\boldsymbol{E}_{j, k, m+1}\left(1-w_{j}\right)\left(1-w_{k}\right) w_{m}\\
&+\boldsymbol{E}_{j+1, k+1, m} w_{j} w_{k}\left(1-w_{m}\right)\\
&+\boldsymbol{E}_{j+1, k, m+1} w_{j}\left(1-w_{k}\right) w_{m}\\
&+\boldsymbol{E}_{j, k+1, m+1}\left(1-w_{j}\right) w_{k} w_{m}\\
&+\boldsymbol{E}_{j+1, k+1, m+1} w_{j} w_{k} w_{m}
\end{align}
磁場についても同様にと
を計算し
\begin{align}
\boldsymbol{B}^{t}=\frac{\boldsymbol{B}^{t-\Delta t / 2}+\boldsymbol{B}^{t+\Delta t / 2}}{2}
\end{align}
によって粒子位置における値を定義します。
電荷・電流の更新
こんどは逆に、ある点に存在する荷電粒子を、格子状に定義された電荷と電流に置き換える作業をします。今回は先ほどと同様に線形で重み付けする方法を用いましょう。
電荷はYee格子の格子点上（図６における直方体の頂点）において定義します。このとき、電荷は格子内部に存在するすべての荷電粒子の効果を足し合わせて
\begin{align}
Q_{j, k, m}&=\sum_{i} q_{i}\left(1-w_{j}\right)\\
&\quad\quad\quad\times\left(1-w_{k}\right)\left(1-w_{m}\right)\\
Q_{j+1, k, m}&=\sum_{i} q_{i} w_{j}\left(1-w_{k}\right)\left(1-w_{m}\right)\\
&\quad\quad\quad\quad\vdots\\
Q_{j+1, k+1, m+1}&=\sum_{i} q_{i} w_{j} w_{k} w_{m}
\end{align}
のように書けます。
また、電荷密度は格子の体積
をもちいて次のように表します。
\begin{align}
\rho_{j, k, m}=\frac{Q_{j,k,m}}{V_{j,k,m}}
\end{align}
さて、図８のように電流の向きを定義します。粒子が格子間を移動することでによって電流が流れます。これは
\begin{align}
\Delta \boldsymbol{w}&=\boldsymbol{w}^{t+\Delta t}-\boldsymbol{w}^{t}\\
\bar{\boldsymbol{w}}&=\frac{\boldsymbol{w}^{t+\Delta t}+\boldsymbol{w}^{t}}{2}
\end{align}
を用いて
\begin{align}
I_{1,j,k,m}&=\sum_{i}\frac{q_i}{\Delta t}\Delta w_1(1-\bar{w}_2)(1-\bar{w}_3)\\
I_{1,j,k+1,m}&=\sum_{i}\frac{q_i}{\Delta t}\Delta w_1\bar{w}_2(1-\bar{w}_3)\\
&\quad\quad\quad\vdots\\
I_{3,j+1,k+1,m+1}&=\sum_{i}\frac{q_i}{\Delta t}\bar{w}_1\bar{w}_2\Delta w_3
\end{align}
のように線形に重み付けして書きます。これらの式によって電荷・電流を更新することができました。
タイムステップ
これでついに図３における１サイクルを回すことができました！最後にの取り方をみていきましょう。
CHPICではCourant-Levyの安定条件である
\begin{align}
\Delta t=\frac{0.95}{c}\left(\sum_{i}\frac{1}{\left(\Delta x_i\right)}\right)^{-1/2}
\end{align}
を提示しています。デバイ長が0.1 mm程度であるとすると、オーダーは sになります。
また、pCANSのように電磁場解析に陰的解法を用いた場合にはこの条件は必要なく、プラズマ振動を記述できるように
\begin{align}
\omega_{pe} \Delta t \leq 0.1
\end{align}
とします。ただしは電子プラズマ周波数で一般にはGHz帯になります。したがって
は
程度になります。
おわりに
今回はPIC法の概要についてざっと説明しました。粒子を広がりのある超粒子として扱うのがPIC法の肝でした。ほかの分野でもこういう考え方はあるんじゃないかな...？
本当は簡単なモデルをシミュレーションしてみたかったのですが、卒論と卒業設計で忙しくて時間が...いや、言い訳はだめですよねスミマセン。春休みにでも実装しようと思っています。
知ってることしか書いてなかったよって人も、よくわからなかったよって人も、最後まで読んでくれてありがとうございました。よいお年を。
参考文献
[1] JAXA, n.d. "小惑星探査機「はやぶさ２」Asteroid Explorer Hayabusa2".
[2] F. F. Chen: Introduction to Plasma Physics and Controlled Fusion, Second edi. Plenum Press. 1974.
[3] Jun Zhou, Dagang Liu, Chen Liao, and Zhenghao Li, "CHIPIC: An Efficient Code for Electromagnetic PIC Modeling and Simulation", 2009.
[4] J.P. Verboncoeur, A.B. Langdon and N.T. Gladd, "An object-oriented electromagnetic PIC code", 1995.
[5] 松本洋介, n.d. "電磁プラズマ粒子シミュレーション", http://www.astro.phys.s.chiba-u.ac.jp/pcans/algorithm.html
