special feature 防護技術の動向
特 集

構造解析ソフトと衝撃解析コードの違いはどこにあるのか
伊東 雅晴（伊藤忠テクノソリューションズ（株））
はじめに

構造解析ソフトとの数値

最後に、本稿で扱う衝撃

ときどき「構造解析ソフ

解法の違いについて述べ

解析コードの計算要素に

トと衝撃解析コードはどてみたい。両コードの相違ついて付記しておきたい。

こが違うのでしょう」とか、点を明示することが本稿

ここではシェル要素やビ「汎用の構造解析ソフトの目的である。

ーム要素に関する数値解を持っているので、それをところで、巷には「誰で法は扱わない。空間的な拡使って衝撃問題を解くこもわかる」という惹句を売がりをもつ物体内部を伝とは可能でしょうか」といりものにしている理系教

播する応力波が、その物体

う質問を受けることがあ

科書があり、その多くは方の変形と破壊に大きく寄る。このような疑問を呈す

程式をなるべく使わずに与する現象の解析に用いる原因として考えられる説明をしていることに特るコンピュータ・プログラのは、構造解析ソフトと比徴がある。しかし、このよ

ムを本稿では衝撃解析コ

較した場合、衝撃解析コー

うな教科書の欠点は、読みードと呼び、その数値解法

ドの特徴、言い替えればそ終えるとなんとなくわかについて論じることが本の数値解法が世の中にそった気持ちにはしてくれ

稿の目的である。

れほど知られていないこるが、重要な事柄についてとが挙げられるのではな

はなんら理解を深める手いだろうか。構造解析の数

助けにはなってくれない

衝撃解析コードの数値解法

値解析手法に関しては数という点にある。そもそも（１） 保存方程式

多くの教科書があるのに方程式を使った説明がわ－連続な場

対し、衝撃の方はといえばかりづらいのは、読者に行

構造解析ソフトも衝撃解

計算コードの理論マニュ間を埋める作業を強いる析コードも、解析の対象とアルを除けば Wilkins のからである。「明らかに次する物体を連続体という数

テキスト [1] 以外には見あ式が得られる」と書いてあ

学的な概念を用いて近似したらないというのが現状るのに、実際に導出してみている点は同じである。連である。

ると数ページの計算が必

続体は無数の物質粒子の集

筆者は 1980 年代の後半要だったという経験は筆まりでできているとみなさから汎用衝撃解析コード者だけに留まるものではれるが、ここで注意しなけ

ANSYS®AUTODYN® の開ないであろう。本稿では方ればならないのは、この物発に係わり、その計算要素

程式を用いて説明を進め質粒子がニュートン力学でや接触条件の数値解法をるが、読者に数式を展開すいう質点という概念とは異開発し組み込んだ経験をる煩雑な作業を強いるこなり、体積をもち、密度やもつ。そこで本稿では衝撃とがないように配慮した応力など熱力学的な概念を

解析コードの側からみた、

つもりである。

付与することができる仮想- 46 -

special feature 防護技術の動向

ラグランジュ表示
ラグランジュ表示

オイラー表示オイラー表示

質量保存式
質量保存式

∂u
ρ&ρ&==−−ρρ ∂ ui i
∂∂xxi
i

∂∂ρρ
∂ (ρ u )==−− ∂ (ρ ui i )
∂∂t t
∂∂xxi

運動量保存式
運動量保存式1 ∂σ ij
u&u&i == 1 ∂σ ijiρρ ∂∂xxjj∂∂(ρ(ρuui )i )== ∂∂ (σ(σij −−ρρuuiuuj ) )
ij
i j
∂∂t t
∂∂xxj
j

1 ∂ (σ ij u i )e&e&== 1 ∂ (σ ij u i )ρρ ∂∂xxjj∂∂   uu2 2 ∂∂ 
u 2 2
ρρee++ == (σ(σij uuj ) )−−ρρuui ee++ u 
∂∂t t  22 ∂∂xxi  ij j
i  22

i 

 


エネルギー保存式
エネルギー保存式

表1i

ラグランジュ表示とオイラー表示による保存方程式的な粒子であるという点で特定の箇所の変形をみるの

標的に衝突して貫通する現ある．ニュートン力学では、に適していることから、こ象を、ANSYS®AUTODYN®

質量と速度の二つの物理量の方法は一般に固体の運動のラグランジュ表示とオイを質点に対して付与するこの記述に用いられる。要す

ラー表示を用いて解析したとができるだけである。

るにオイラー表示の計算メ

ときの現象図を示す。解析

さて、連続体の運動の記

ッシュは時間が経過しても体系は２次元円柱座標であ述は、この物質粒子の存在

空間に固定されて動かないる。飛翔体の貫入距離は両を明示的にとらえるかどうが、ラグランジュ表示の計表示でほぼ同じ結果となる。

かでオイラー表示とラグラ

算メッシュは空間内を移動

一般的に構造解析ソフトでンジュ表示という二つの方するという違いが存在する。

固体を模擬する場合にはラ法に分けられる。オイラー

表１に示したのはラグラ

グランジュ表示を適用し、表示では、物質粒子が意識ンジュ表示とオイラー表示

オイラー表示を用いることされないで、空間に固定した座標をもちいて空間内の特定の１点における物理量の変化が観測される。流体の定常流れではそもそも特定の物質粒子という考えが意味をなさないことから明らかなように、この方法は一般には流体の運動の記述に適しているといえる。一で 表 現 し た 保 存 方 程 式 [2]である。気体も液体も固体もすべての物質の力学現象
はこの一組の方程式に従う。
こ こ で ρ は 密 度 、 ui
（ i = 1, 2, 3 ）は速度ベクト
ル、 σ ij （ i , j = 1, 2, 3 ）は応力テンソル、 e は内部エネルギー密度を表す。本稿ではない。しかし、衝撃解析
コードでは解析する対象に応じてラグランジュ表示とオイラー表示のいずれかを選択して用いる。
ここから本論に入る。表１に示した保存方程式が解をもつためには、式中で使われている密度、速度ベク

方、ラグランジュ表示では、

は、場の変数はすべてデカ物質粒子に固定された座標

ルト座標系 xi （ i = 1, 2 , 3 ）物理量が連続で、かつ微分をもちいて運動が記述されを用いて表現することとす

可能でなければならない。

る。この場合は物質粒子に物理量に関するこの連続条

乗って物理量の変化を観測る。
図 1 に、鋼鉄製の飛翔体することになる。連続体のが高速でアルミニウム製の解析に用いられるのが構造- 47 -

トル、応力テンソルなどの件が満足されている現象のspecial feature 防護技術の動向

ラグランジュ表示
図1

オイラー表示二つの表記法を用いた高速衝突の数値シミュレーション（使用計算コード ANSYS®AUTODYN®）

解析ソフトである。

は、現在の計算機の処理能に代わり、衝撃波の前後ででは、その条件が満足さ力では不可能である。何ら物理量が満足しなければなれないと解を得ることがでかの条件を加えなければ衝らない一組の方程式である．

きないことになるが、物理

撃波の前後で物理量は不連

それは以下に述べる方法で量の滑らかさが常に満足さ続となるために、保存方程導くことができる。

れるという保証はない。そ

式が成り立たず、数値解は

いま粘性のない完全流体の一例が衝撃波である。衝不安定となってしまう。そで満たされた空間領域の中撃波の厚さは、気体の場合の条件が次ぎに述べる跳躍を、衝撃波が定速 U で伝播

で分子の平均自由行程（例

条件である。衝撃解析コーしている状況について考え

え ば 圧 力 100kPa 、 温 度

ドには、安定な数値解が得る。熱の発生はないものと300K の窒素で約 0.1μm）

られるように、跳躍条件にする。図 2 に示したように、の数倍から 10 数倍程度と基づくいくつかの数値解法

衝撃波の波面に垂直な方向

非常に薄い。圧力や密度のが組み込まれている。

に x 軸をとり、それに直角
な y 方向と z 方向には物理

空間勾配を計算するために（２） 跳躍条件とは、衝撃波の厚さ方向に数量が変化しないと仮定した個の計算メッシュが必要でランキン・ユゴニオ

１次元の平板体系で議論をあるが、そのサイズの計算

関係式－不連続な場進める。

メッシュで解析を行うこと

跳躍条件は、保存方程式
衝撃波yu0 = 0u1 = u

z

ρ1

U

p1

x

（ 非定常な状態）

一次元衝撃波の模式図- 48 -

p0（ 静止した状態）

x1 x0
図2

ρ0

special feature 防護技術の動向
垂直衝撃波に対するランキン・ユゴニオ関係式
ρ 0U = ρ1 (U − u )p1 − p0 = ρ 0 uU1
( p1 + p0 ) 1 − 1 2
 ρ 0 ρ1 

e1 − e0 =この式からは衝撃問題の解析に適した衝撃ユゴニオ状態方程式や、衝撃波近傍の数値解の安定化に寄与する人工粘性を導くことができる

１次元の平板体系のオイ

っ て い く 極 限 操 作 x1 → x0

ラー表示の保存方程式は、を行う [3] 。このとき積分の表１に示した３次元の保存

それぞれの項について

方程式を１次元に変換する
ことで容易に得られる。速

x

1
⌠ ∂
( ) dx → 0lim 
x1 → x0 ⌡ ∂ t
x0

度 ベ ク ト ル ui を ス カ ラ ー
量の流速 u で，そして応力
テン ソ ル σ ij も スカ ラ ー 量の圧力 p で置き換え、
質量保存式

・・・(1)

理量が満足しなければならない状態を定めるものである。しかし、注意して見直すと、これらの式には衝撃

・・(4)

x1

⌠ ∂
( ) dx = ( )x1 − ( )x0lim x1 → x0 ⌡ ∂ xx0・・・(5)
が成り立つことに注意する。
式（4）から式（1）~（3）∂ρ∂= − (ρ u )∂t∂x

静止状態と非定常状態の物の左辺を積分した項はすべてゼロに等しい。そして、

波の情報すなわち衝撃波速度 U が含まれていないことに気づく。そこで、座標系を空間に固定されたEuler 座標系から、一定速度 U で進む衝撃波に貼りつけた座標系に変換して、衝撃波の速度を跳躍条件に取り込むことを考える。

それぞれの式の右辺の積分

衝撃波からみると、衝撃)項に式（5）を適用するこ波前面の静止状態の物質粒とにより、以下の表現を得・・・(2)

子は u 0 = −U の速度で飛びる。
(ρ1u1 ) − (ρ 0u0 ) = 0 ・・・(6)22
p1 + ρ1u1 − p0 + ρ 0u0 = 0

運動量保存式∂(ρ u ) = − ∂ p + ρ u 2
∂t
∂x

(

エネルギー保存式∂ 
u2  ρ e + ρ  =
∂t 
2

∂ 
u3
 ρ eu + ρ+ pu −∂x2
・・・(3)
を得る。

(

) (

)

・・・(7)2p u 
 e1 + 1 + 1  −ρ1 2 


係を式（6）と式（7）に代入して整理すると
ρ 0U = ρ1 (U − u ) ・・・(9)・・・(10)を得る。両式を U と U − u

を衝撃波の一方の端面（ x = x 0 ）から他方の端面・・・(8)（ x = x1 ）まで積分し、同式（6）～（8）が跳躍条

時に衝撃波の幅をゼロにも件であり、衝撃波を挟んで- 49 -

速度で去っていく。この関p1 − p0 = ρ 0 uU

2

p
u 
 e0 + 0 + 0  = 0

ρ0
2 


次ぎに、これらの方程式

込んできて、u1 = −(U − u ) のについて解くと1  p0 − p1 U=
ρ 0  (1 ρ 0 ) − (1 ρ1 ) 

12・・・(11)

special feature 防護技術の動向
固体に力が加えられて結圧力と同様に線形的な和と晶格子がひずむと、結晶内して表すことができるもの

部にはひずみに応じた体積とする。

・・・(12)

変化により圧力が生じる。

が得られるので，これらを状態方程式はこの圧力の大e = e L + eT

式（8）に代入して整理す

きさを定める方程式である。

上式において eL と eＴは、ることにより、

保存方程式や跳躍条件がす式（14）の圧力の各項に対

 1
1
1 
− 
e1 − e0 = ( p1 + p 0 )
2
 ρ 0 ρ1 ・・・(13)

べての物体について平等に応する内部エネルギー密度

成り立たなければならないを表す。

のに対し、状態方程式は物

ここで Grüneisen が提を得る。式（9）、式（10）、質に固有で、例えば硬いと示した関係式 [5] 、すなわち式（13）が垂直衝撃波に対か軟らかいというような性

pT と eＴの間に線形関係がするランキン・ユゴニオ関

質はこの方程式により表現成り立つとした。

係式である。

される。

1  p0 − p1 U −u =
ρ1  (1 ρ 0 ) − (1 ρ1 ) 

12

pT =

Γ
eT
V

・・・(15)・・・(16)

このランキン・ユゴニオ

結晶格子がひずむとき体の関係式は重要な式である。

積が変化するだけでなく、を用いると、圧力は次式の以下で具体的に示すように、

一般的に温度も上昇するのように表現することができ

この式からは衝撃問題の解で、体積変化からの寄与 pLる。

析に適した衝撃ユゴニオ状と温度上昇からの寄与 p Ｔ態方程式や、衝撃波近傍のの和として、圧力 p は次式

数値解の安定化に寄与するのように表される[4]。

Γ(e − eL )V・・・(17)
ただし、V (=1/ρ) は比容
積、 Γ はグルナイゼン（Grüneisen）係数と呼ばれ、物質に固有な定数である。この式の意味するとこ

人工粘性を導くことができる。

p = p L + pT・・・(14)

また、圧力が行う仕事によ（３） 衝撃ユゴニオ

り物体には内部エネルギー

状態方程式

e が蓄積されるが、それはp = pL +Mie-Grüneisen 型衝撃ユゴニオ状態方程式
Γp = pH + (e − eH )V

ρ 0 c0 µ (1 + µ )
2

pH =

eH =

[1 − (s − 1)µ ]2

1 pH  µ 
2 ρ 0  1 + µ 

Γ はグルナイゼン（Grüneisen）係数と呼ばれ、物質に固有な定数
状態方程式は圧力の大きさを定める方程式である
保存方程式や跳躍条件がすべての物体について平等に成り立たなければならないのに対し、状態方程式は物質に固有で、例えば硬いとか軟らかいというような性質はこの方程式により表現される- 50 -

special feature 防護技術の動向
ろは、圧力が体積変化だけに依存する項（第 1 項）と温度変化に依存する項（第2 項）の和として表すことができるということである。
式（17）は Mie-Grüneisen型状態方程式と呼ばれる。
式（17）で表される状態
方程式の中でもっとも単純な形は、温度依存性のない（Γ=0）線形等方弾性体に対する次の状態方程式である。

p = Kµ・・・(18)図3

この式は構造解析では広く

2024-Alminum の衝撃波速度 Usと粒子速度 Up の関係[6]用いられているが、これを高速衝突問題の解析に適用することには精度の点で問題がある。衝突速度が大きくなると、体積変化 µ に対

e1 − e0 =1( p1 + p0 )  1 − 1 2
 ρ 0 ρ1 ・・・(21)

関係式は、いま唐突な形でここに登場させたが、多くの物質に認められる性質であり、文献[6]には 600 種以する圧力の変化が線形関係であった。いま静止状態のから外れていき、さらに温物理量（ ρ 0 ， p0 ， e0 ）の 3

度上昇の影響も無視できなつは既知とすると、値が未れている。一例として、図くなると予想されるからで知な変数は（ ρ1 ， p1 ， e1 ，3 に 2024 アルミニウムに

u ， U ）の 5 個である。し

ついて得られた実験データ

たがって次式を示す。

ある。
そこで、前に導いておいた、垂直衝撃波の前後の物
理量に関するランキン・ユ
ゴニオの関係式を利用して、衝撃波の存在下でも圧力がU = c0 + s u・・・(22)を加えて方程式の個数を 4上の物質のデータが掲載さ

式（19）～（21）において，衝撃波が伝播する前の静止状態についてつとすると、p1 と e1 を ρ1 の

p0 = e0 = 0 を仮定しても一関数として表すことができ般性は失われない。そうしる。それがここで求めようて得られたランキン・ユゴとしている衝撃ユゴニオ状

ニオ関係式と式（22）を連係式は、態方程式に他ならない。式

立させると、次式を導くこρ 0U = ρ1 (U − u ) ・・・(19)（22）は、衝撃波の速度 Uとができる。

精度良く評価できる状態方
程式を以下の手順に従い求める。
ランキン・ユゴニオの関

p1 − p0 = ρ 0 uU・・・(20)と運動する物質粒子の速度

u の線形関係を示している。
c0 と s は試験により決められる材料定数である。この- 51 -

ρ 0 c0 µ (1 + µ )2

p1 =[1 − (s − 1)µ ]2 ・・・(23)special feature 防護技術の動向
e1 =1 p1  µ 
2 ρ 0  1 + µ ・・・(24)

れている。

数の衝撃波が互いに干渉し（４） 衝撃波を緩和する

それぞれの衝撃波を追跡し、つつ伝播している状況で、Mie-Grüneisen 型 状 態 方

人工粘性力

その位置と形状をとらえる

程式（17）の（ PL ， e L ）上でも述べたように、衝には非常に複雑な数値計算に、求めた（ p1 ， e1 ）を代

撃波の厚さは実用上可能なアルゴリズムを必要とする

計算メッシュのサイズに比からである。

入して、

p = pH +Γ
(e − eH )・・・(25)Vρ 0 c0 µ (1 + µ )2pH =[1 − (s − 1)µ ]2 ・・・(26)

1 pH  µ eH =2 ρ 0  1 + µ 

較すると非常に薄いので、他の一つは、Von Neu-

その厚さ方向に数個のメッ

mann と Richtmyer が考案

シュを切ることは不可能でした人工粘性力 [8] を用いるある。

方法である。この手法は衝問題を解決する一つの方

撃波を追跡する必要がなく、法は、衝撃波を相互作用境

また複雑な計算アルゴリズ

界面として扱うものである

ムも使わず３次元への拡張[7] 。ランキン・ユゴニオ関

も比較的容易なため、発表

係式（19）～（21）と状態されてから半世紀が経過し方程式（25）から、圧力・た現在でも多くの汎用衝撃

速度・密度・内部エネルギ

解析コードに採用されていーなどの跳躍量を求め、次る。なお、このことは当然にこの跳躍量を衝撃波の両であるが、衝撃波の伝播を側の領域に境界値として負解析対象としない構造解析

る多くの物質に対して、こ荷するという手法である。

ソフトに人工粘性力は不要の状態方程式が用いられる。

この手法の欠点は、１次元である。

ANSYS®AUTODYN® の 材

よりも高次の次元に拡張す人工粘性力の基本的な機ることが難しい点にある。

能は、衝撃波の厚さを計算

すなわち３次元空間内で複メッシュ幅の 3～5 倍にぼ・・・(27)を得る。ただし、（ p1 ， e1 ）を（ p H ，eH ）と表記した。
これが Mie-Grüneisen 型
衝撃ユゴニオ状態方程式である。衝撃解析で模擬され料ライブラリには 50 種類
ほどの材料の定数が記載さ

衝撃波を緩和する人工粘性力2

q1 = ρ CQ l 2 ε&kk2

q2 = − ρ C L l c0ε&kk
q = q1 + q2
衝撃波の厚さは実用上可能な計算メッシュのサイズに比較すると非常に薄いので、その厚さ方向に数個のメッシュを切ることは不可能
問題を解決する一つの方法として Von Neumann と Richtmyer が考案した人
工粘性力を用いる方法である
この手法は衝撃波を追跡する必要がなく、また複雑な計算アルゴリズムも使わず３次元への拡張も比較的容易なため、発表されてから半世紀が経過した現在でも多くの汎用衝撃解析コードに採用されている- 52 -

special feature 防護技術の動向
かし、基礎方程式の解を空の形で人工粘性力を評価す性項を導く方向で議論を進間内のすべての場所で連続る手法を提案した。多くのめる。

にすることにある。ただし、汎用衝撃解析コードが現在

人工粘性力の影響は衝撃波

この式を採用している。

式（33）と式（35）から得られる次式が存在しない領域では、本

ここで人工粘性力を実際

来の解に影響を与えないよに 導 い て み る 。 von Neu-

うに、無視できる程度に小mann と Richtmyer の原著・・・(36)[8] に は そ の 導 出 過 程の右辺の ρ1 に、式（32）か

さくなければならない、

論文

ρ 0 uU = K (ρ1 ρ 0 − 1)

Von Neumann と Richt-が示されておらず、またmyer はそのような条件を

Wilkins [10]は理想気体の場

満足する人工粘性力の関数

合について人工粘性力を導

ら導いた
U
ρ1 =
ρ0
U −u

として、次式を提示した。

いたが、論文にはその結果を代入して整理すると

2 2

q1 = ρ CQ l ε&kk

2・・・(28)

ここで l はメッシュ幅である。CQ は定数を表し、この値を大きくすると数値的な衝撃波の厚さは厚くなる傾向を示す。ε&kk は変形速度テンソルの対角成分の和、すなわち体積変化率である。

ε&kk は衝撃波の近傍では大きく、衝撃波から離れるのにしたがい小さくなるので、
人工粘性力に用いるのに望ましい性質を有している。
しかし、式（28）で表さが記されているだけである。
このように導出過程を記述した文献の入手が難しいことを考慮して、以下では線形弾性体の場合について、
ランキン・ユゴニオ関係式と状態方程式を用いて、人工粘性力の q1 と q2 の関数形を具体的に導き、それが式（28）や式（29）と同じ形になることを確認する。
いま衝撃波が到達する前の状態量について・・・(31)p0 = e0 = 0

れる人工粘性力は、衝撃波と仮定し、この式を式（19）前面に対しては有効である～（21）に代入するとが、衝撃波が通過した後のρ 0U = ρ1 (U − u )・・・(32)

数値的振動を抑制するのに

p1 = ρ 0uU

・・・(33)は有効に機能しない欠点のあることがその後指摘された。この欠点を解決するた
めに、Landsoff [9] は ε&kk に関して線形な関数 q 2

q2 = − ρ C L l c0ε&kk・・・(29)e1 =p1  11 −  ・・・(34)2  ρ 0 ρ1 

2・・・(38)を得る。ただし c0 = K ρ 0を用いた。 c0 は物体中を伝播する音速を表す。
式（38）は U に関する２次式であり、容易に解くことができて、その解は
12
2u  u 2
U = +   + c0  ，
2  2 

(QU > 0)・・・(39)である。これを式（33）に代入すると次式を得る。
12 u  2
u2
2
p1 = ρ 0+ ρ 0 u   + c0 2 2 ・・・(40)
ここで u と c0 の大きさが極端に異なる２つの場合について考えると、式（40）を得る。
以下では、式（32）～（34）と次の線形弾性体の状態方
程式

を q1 と組み合わせ

p1 = K (ρ1 ρ 0 − 1) ・・・(35)q = q1 + q2を連立させて U について・・・(30)U (U − u ) = c0

・・・(37)

解き、それを用いて人工粘- 53 -

から次の２つの式を導くことができる。
p1 = ρ 0u 2 ， (Q u >> c0 )・・・(41)p1 = ρ 0 c0u ， (Q u << c0 )・・・(42)
この２つの式はそれぞれ２

special feature 防護技術の動向
次と１次の人工粘性、すな

スライドライン法に大きく

ドライン法については近似

わち分けられる。

リーマン・スライドライン

2

q1 = ρ 0CQ u
q2 = ρ 0C L c0u

・・・(43)・・・(44)に対応する。
以上の議論では１次元を仮定して式（43）と式（44）を導出したが、これを３次元に拡張する場合は、速度u を (lε&kk ) に 置 き 換 え る ことが必要である。その結果
は式（28）と式（29）に一致する。

接触境界条件
ラグランジュ表示を採用
した計算コードで物体の衝突問題を解析するには、接触面での計算メッシュの重なりを防ぐために接触境界条件が必要である。図 1 に示した例題では、飛翔体と標的の間に接触境界条件が設定されている。この接触境界条件の計算手法は、その開発時期と適用分野の違
いからいわゆる節点力法と節点力法は、1980 年代以法[13]をとりあげ、両者の計降に開発された非線形構造

算精度の比較をする。詳し

解析コードで一般的に採用

い説明は省くが、近似リーされている手法である。そマン・スライドライン法は、の特長は、アルゴリズムが接触面の速度を上述のラン

スライドライン法に比べる

キン・ユゴニオ関係式からと簡潔で計算時間が短くて求める手法であり、接触面

すむ点である。自動車衝突において透過あるいは反射や塑性加工など衝突速度がする圧力波を精度良くとら

比較的低い現象の解析に適えることができる。

している。一方のスライド

解析対象は無限長さの２ライン法は、開発された時本の同軸円管である。図 4

期が節点力法よりも早く、に示したように、解析体系

1950 年代の HEMP コード

は対称性を考慮して、実体[11]まで遡る。その特長は接系の 1/4 とし、z 方向の変

触面を通過する応力波を正位は拘束する。内側の円管しく模擬できることである。

には解析終了時まで、そのこの計算手法は、接触面に降伏応力より大きな一定のおいて計算要素と計算要素

圧力を負荷する。円管は弾の接触面積を求め、接触面

完全塑性体であり、以下のでの応力波の透過・反射が

材料定数を用いる:密度

精度良く模擬できるように

2700 kg/m3、体積弾性係数

考慮されている。

70GPa、横弾性係数 30GPa、以下において、節点力法

降伏応力 300MPa。

の代表的な計算手法である

図 5 に、ペナルティ法とペナルティ法[12]と、スライ近似リーマン･スライドラ

15

z

y
9. 0

10

7. 0
5. 0

Penalt y

1. 0

Appr ox.
Rieman nslide line

5

x
0. 0

図4

(10 -2 m)-15

同軸円管の衝突問題-10

図5- 54 --50

5円管のメッシュ変形形状

1015(10 -2 m)

special feature 防護技術の動向
イン法で得られた t = 30 µs

妥当な解析結果となった。

境界条件の解法の３点に注における変形図を示す。解ペナルティ法の場合、この目して論じた。保存方程式

析条件に対称性を崩すよう解析例のようにメッシュ依の極限操作から跳躍条件をな条件は設定していないの存性が発生する可能性があ経てランキン・ユゴニオ関で、円管断面はつねに円形るので問題によっては適用

係式へと進み、ランキン・

形状を保ちつつ変形するはする際に注意が必要である。

ユゴニオ関係式からは衝撃
ユゴニオ状態方程式と人工ずである。しかし、ペナル
ティ法（左側）で得られた粘性力を具体的に導き、そ変形形状は周方向の対称性

まとめが崩れ、折れ曲がる形状と構造解析ソフトと衝撃解なった。それに対し、近似

析コードの違いについて、

リーマン･スライドライン

本稿では、圧力を規定する

法（右側）の計算メッシュ

状態方程式、衝撃波を緩和

は変形後も円形が保たれ、する人工粘性の有無、接触の導出過程を明示した。

参考文献[1] Wilkins, M.L., Computer Simulation of Dynamic Phenomena, Springer, 1999.
[2] Eliezer, S. et al., Fundamentals of Equations of State, 2002, World Scientific, pp.165-168.
[3] Eliezer, S. et al., Fundamentals of Equations of State, 2002, World Scientific, p.179.
[4] Eliezer, S. et al., Fundamentals of Equations of State, 2002, World Scientific, p.153.
[5] Grüneisen, E., Handbuch der Physik, ed. H. Greiger & K. Scheel, 10, 1-59, Springer.
[6] Marsh, S.P., LASL Shock Hugoniot Data, University of California Press, 1980.
[7] Hiermaier, S., ed., Predictive Modeling of Dynamic Processes: A Tribute to Professor Klaus Thoma, 2009, Springer, p.349.
[8] Von Neumann, J. and Richtmyer, R.D. A method for the calculation of hydrodynamic shocks, Journal of Applied Physics, 21(1950), pp.232-237.
[9] Landsoff, R., A numerical method for treating fluid flow in the presence of shocks, LA-1930, Los Alamos National Laboratory.
[10] Wilkins, M.L., Use of artificial viscosity in multidimensional fluid dynamic calculations, Journal of Computational Physics, 36(1980), pp.281-303.
[11] Wilkins, M., Calculation of elastic-plastic flow, UCRL-7322, Rev.1.
[12] Hallquist, J.O., Goudreau, G.L. and Benson, D.J., Sliding interfaces with contact-impact in large-scale lagrangian computations,
Computer Method for Applied Mechanics Engineering, 51 (1985), pp.107-137.
[13]伊東雅晴，村上澄男，近似リーマン法に基づく三次元動的接触境界アルゴリズム，日本機械学会論文集 A 編，66-649(2000)，pp.130-137.

- 55 -