セプタム磁石ビームチューブの解析

2009/6/16

和気正芳

KEKDAのセプタム磁石はアウトガスを低減するために磁場中にビームチューブを入れることになった。この場合、ビームチューブには渦電流による発熱があり、大気圧と熱変形応力がかかる。設計者川久保さんの依頼で、このビームチューブの有限要素解析を行ったので報告する。

磁場、渦電流、大気圧、熱変形のすべてを含めた連携解析は単機能のプログラムでは出来ず、マルチフィジックスを謳う有限要素法の昨今の流行でもあるので、この報告には若干のANSYSでの手法解説を含めて、後日必要になる他の磁石解析の参考となるようにしておく。

1 磁場計算

ANSYSでは幾つかの種類の3次元磁場計算法が用意されているが、これらは全て2次元とは全く違った計算方法になる。今回はベクトルポテンシャル法を用い、SOLID97要素を使った。解析の対象はビームチューブであり、磁石自体は必要ないので、空間にベクトルポテンシャルを設定して一様な磁場空間を作った。こうすれば鉄を使わなくてすむので、計算は非常に早くなる。ANSYSではこのように作った磁場も時間変動させて計算することが出来る。実際の運転はパルス波形であるが、簡単のために交流磁場の半波を使うという解析にした。結果的に計算時間はFermilabの64bitPCを使えば20分以内になった。 しかしながら、安達さんが行ったOPERA3Dの結果と比較してみると、Eddy Currentの作る磁場には大きな違いがあることがわかった。結論として、今回の解析にはこれで良いが、今後Eddy Currentが作る磁場を計算する場合にはこの方法は問題があることになる。

2 要素分割

モデルを要素分割するが、収束性からも、精度からもなるべく同じような大きさの直方体に近い形状のメッシュを切るのがよい。三角メッシュや扁平な要素を作らないように気をつける。ANSYS3次元磁場計算は実用に耐えないなどと言う人もいるが、事例を調べて見ると多くはメッシングに問題があったと思われる。適正なメッシングが行われないと計算は収束せず精度も悪い。このためにはなるべく6面体の集合でモデルを作る必要があり、モデル作成段階で配慮する必要がある。近年はCADファイルを自動的にメッシュするといったことが行われているが、このような安直なメッシングでは精度と収束性が悪く、後で述べる検証に大きな問題となる。メッシングは単なるCPUの節約ではないという認識が大切だと思う。ブーリアン操作はモデル作成に有用だが、多用するとMKS単位系では分解能がなくなり失敗する。mm単位系で作るとモデルしやすいが、今回のような連携解析の場合、単位系の整合性を極めるのが難しく、結局MKSに戻ってモデルをつくる工夫をした。図のようなメッシュで、楕円断面の筒の上半と、その周りの空間を構成した。見やすくするため要素サイズを大きくしてあるが、実際にはこの五分の一のサイズにした。いくつか試して見たがこのメッシュが最も収束の早いものだった。全て6面体要素になっており、4面体は使われていない。

3 渦電流

磁場は外から一様な6.6kGをかけたが、磁場変化に対してビームチューブにはEddy Currentが流れる。その様子は図に示すように、上面から見て回転する電流であり、中央部にはほとんど流れず、周辺を回るように流れる。SOLID97はEddy Currentの自由度を持つ要素と持たない要素に分けて、空間にはEddy Currentを流さないようにする。こうしないと空間のわずかな電流も計算して計算時間が長くなる。交流解析を選択すればこの設定で自動的にEddyCurrentは算入される。交流であるから、磁場の状態を表示すると、そのままではImaginaryPart表示されるのでload stepを戻ってReal Partを表示するようにする。Currentそのものはout phaseだからImaginary partになる。

ビームチューブの中の磁場は、渦電流で遮蔽されるので、少し下がるが、結果的には大きなものではない。逆にビームチューブの少し外側では少し磁場が上がる。mid plane中央部の磁場分布を示すと、図のようになる。Real Partでの磁場のずれは最大で1.8GだがImaginary Partでは100G位の磁場がでている。Imaginaryの部分も位相が90°ずれるだけだからパルス前半後半に100Gの磁場ズレが起こると考えることになる。この結果については安達さんのOPERA3D計算結果300Gとはかなり違う。それでもジュール熱に関しては一致しているということになった。

4. ジュール熱

渦電流が計算できれば、抵抗率でジュール熱が計算できるがANSYSのSolid97の場合ジュール熱についてはPOWERHというマクロコマンドが用意されており、これだけで計算が済む。しかし、このコマンドはMKSA単位系以外には動かないようになっているようだ。

ジュール熱の分布は両端では電流の分布にしたがい、上下面中央部ではmid planeに集中している。2乗だから電流分布よりもさらに極端になる。この図は見やすくするため実寸よりも短くしたモデルの結果を示してある。

5 検証

これまでの計算の精度については実は何も確証がない。幾つかの計算が重なっている上に、それぞれの計算にはパラメータの受け渡しもあり、この解析でも最終結果を得るまでに何度も間違いを見つけた。間違いの主なものは

などであるが、実験の場合、こういった間違いがあれば結果が出ないのだが、計算ではともかくも結果が出てしまう。二重にまちがった結果、もっともらしい値を得るということも珍しくない。だから論文となっている報告でも、学生に計算させた場合などはあやしいものが多い。いつでも計算精度のせいにして、不可解な点をごまかせるからである。だからメッシュの取り方が悪く、計算精度が出ず、計算時間が長い場合は、こうした間違いを助長してしまう。議論の結果、この計算では、解析解との比較をするために同じ幅、長さの平板モデルを作った。さらに平板を少しづずつ曲げてビームチューブの楕円断面にして行くようにして計算をやって見た。こうすれば、連続性で見てどこかで変な操作をしていれば見つけられると思われるからである。図のなかでhitとあるのは楕円の短半径である。pavgは上半面での発熱量だが、楕円の形状変化で平板に近づいていることがわかる。平板については川久保さんが式(3)として提示されているが、同じ式を交流磁場にあてはめて

P=

を用いて計算して、duty factorを掛けて求めると84.9wであるから、82.124Wは、まあかなり良い一致が得られたと言える。周波数は150Hzで、10Hzの半波繰り返しがduty cycleである。hitはスクリプト上でパラメータとして入れ込み、変化で他の設定に影響がないようにしてあるから、変な間違い要素は入り込まないと思う。

平板から徐々に楕円ビームチューブに近づいて行き、値も変化していることがわかる。これで解析解との連続性が保証されるのであるていど確実な結果と言える。短直径24.5mmの楕円ビームチューブの場合、全部で240Wの発熱になった。安達さんがOPERA3Dで計算した結果と比較するために長さを50cmした場合は2*1920W*10/300=128Wである。安達さんの計算結果は対称性を使っているので、半分の長さで885w*2*2* 10/300=2*1770W*10/300=118Wだから、まあ一致していると言える。

6. 熱計算

のようにして求めたジュール熱を使って今度は熱計算をするが、それにはSolid5を用いる。Solid5は汎用の要素で実は磁場計算も出来るのだが、静磁場にしか対応していない。この操作はこれまでのノードをそのまま使って要素をそっくり入れ替えることで行う。ジュール熱はマトリックスに退避させておいて新しい要素に荷重として与える。空間の要素は必要がないのでSolid0のヌル要素にしてしまう。水冷を正しく評価するためには水の要素も入れなければならないが、磁石は短く高圧の水があるので十分な流量が取れるので水冷配管部分の温度を固定することで十分である。セプタムの内側は導体を固定するアンカーが入るのでスペースが少ないので、外側のmid planeだけに一本の水冷管を配置した場合の結果を示すと図のようになる。

温度の表示は水温からの温度上昇である。当然のことながら内側が温度上昇するが、これくらいの温度ならば材料的には問題がないと思われるが、271度は、やはり少し高い温度ではある。実際には空冷効果もあり、伝導で周りにも拡散するので材料的には多分大丈夫だと思われる。

7. 大気圧による変形

Solid5では変形計算も出来て、大気圧101300Paによる変形が計算される。拘束条件として両端は大きなフランジで動かないとしているので、中央が扁平になる方向に変形する。図では中央部での断面が見えるように半分に切った状態を示してある。長さ910mmの副セプタムの場合、肉厚は1mmであるが変形量は0.12mmなので、大気圧による変形もたしかにあるが、ストレスとしては1Gpa以下なので、力学的に破断するようなことはない。

8. 熱膨張による変形

実際には熱膨張による変形が大気圧による変形に重なる。片側が水冷配管で冷やされた場合、ビームチューブはいびつになり、熱膨張の小さい側に曲がりが生じる。総計で、0.5GPa程度の応力が働くが破断限界には程遠い。図で両端に現われている大きな応力は拘束条件で動きを固定したためで、無限に硬い熱伝導の無いフランジを仮定したことになる。実際にはフランジにも熱が伝わり、フランジも変形はするのでこのような応力が出ることはない。

変形は横方向に0.4mm。これは熱膨張による曲がりが大きい。温度の高い内側は楕円が膨らみ断面は少しいびつになる。こうした非対称な変形は冷却が非対称であるために必然的に起こる。

縦方向の変形は主として大気圧によるつぶれが最大0.094mm、これは熱膨張の広がりで相殺されて数値としては熱膨張を考慮しない場合よりちいさくなる。両端部は曲がりのしわ寄せで0.1mm持ち上がることになるが、実際は完全固定ではなくフランジなので軽減される。

9. 内側への水冷管増設

以上のように、厚さ1mmで片側水冷のデザインでも製作は可能だと考えられる。しかし、内側の上下にも水冷配管が入れられないこともないので、これを入れた場合も計算する。こちらは太い配管を入れられないので水温は20度高くなるとしている。この水冷管の配置効果は大きく熱分布は大幅に改善される。熱膨張による曲がりもこの場合非常に小さくなる。半分に切った状態で示し、歪は拡大表示されているが、温度上昇はわずかに42度に留まる。熱応力も格段に少ない。曲がりを防ぐ対称性と言うことからすれば4本の水冷管のほうが望ましい。

渦電流磁場

外部磁場をかけるためにこの解析ではベクトルポテンシャルを与える方法を使い、早い計算を行った。しかし、OPERA3Dの結果とはEddy Currentの作る磁場を一致させることが出来なかった。これは計算コードの問題ではなく、外部磁場の掛け方によるものである。安達さんの解析では、20時間もかかる計算になったが、導体を設定して電流を流し、鉄で境界条件を与える計算を行った。ベクトルポテンシャルを与える方法では10分で終わる速さがあったが、Eddy Currentの磁場は正しくない。フラックスラインはベクトルポテンシャルの等高線のようなものだから、両端のベクトルポテンシャルを固定すれば、一様磁場が現れ。計算のフラックス状態はRealPartで

のようになり、Eddy Currentはこの磁場に基づいて計算される。しかしながら、Eddy CurrentそのものはImaginaryPartであり、外部磁場のImaginary Partは無い。適当な電流を想定してフラックスを表示すると

のようになる。Eddy Currentによる磁場はかなり広がっていて境界でのベクトルポテンシャルは一様ではない。交流磁場でベクトルポテンシャルを振動させた場合実は、Imaginary Partが「無い」のではなく、「ゼロ磁場を与える」ということになる。フラックスは

のようになり、宇宙の彼方までの磁束をすべて境界内に濃縮してしまうため。全体のオフセットがついてしまうことになる。このためよほど大きな空間をメッシュしないとImaginary磁場は正しくなくなる。Eddy Currentそのものは正しい。

まとめ

KEKDAのセプタムビームチューブについて、渦電流、大気圧、熱膨張、熱伝導を考慮した解析を行った。結論としては肉厚1mmで製作は可能であり、外側の一本水冷でも持ちこたえることが出来るが、内側にも水冷管を入れた場合その効果は非常に大きいことがわかった。この解析の中で、有限要素計算をどのようにして実際のものづくりに役立てられるかを考えることができた。計算過程に検証を入れることが大切である。