動的緩和法による膜構造物のしわ解析

柏 宗孝*1
小野田 淳次郎*2

梗

概近年、動的緩和法を膜面のしわ解析に適用した研究が行われている。しかし、従来の動的緩和法を用いた方法では収束解を得るまでの計算コストが大きいという問題がある。本研究では従来の動的緩和法に
静的な反復解法を組み合わせることによって動的緩和法の持つ安定性を維持しながら解析の収束性を向上させる方法を提案し、その有効性を数値解析によって従来までの手法と比較することで確認する。

1． はじめに
近年、膜面構造物は、収納の容易さや軽量性から、宇宙用の構造物として注目を浴びており、膜面構造物を使用した宇宙ミッションが数多く提案されている。しかし、膜面は非常に小さな曲
げ剛性しか持たないため、わずかな圧縮荷重が作用しただけ
で面外方向にしわ状の変形を生じる。膜面に生じたしわは、膜面の形状精度を低下させるだけでなく、膜面の振動モードを変化させることが知られている 1)。膜面の形状精度や振動特性は、
膜面アンテナや膜面鏡などのミッションにおいては、ミッションの実現性を左右するほど重要な設計要素である。そのため、膜面に生じたしわが膜面の挙動に与える影響を精度良く予測するための解析手法の確立が求められている。
これまでに膜面のしわ解析のための手法として数多くの手法が提案されている。解析的に平衡解を求める手法もいくつか提案されているものの、その多くは有限要素法に基づいた数値解析手法である。また、膜面のしわ解析といった場合に、その解析手法の安定性と収束性は、手法が背景とするしわの表現理論の安定性と収束性と、その理論に基づいて平衡解を得るための数値解析法の安定性と収束性から決定される。そこで、本稿ではしわの表現理論と平衡解を得るための数値解析法の 2 つを分けて考えて、それぞれについての既往の研究についてまとめる。
まず、これまでに提案されているしわの表現理論は、大きく分
けて張力場理論と分岐座屈理論の 2 つに分けられる。張力場理論は、膜面内の応力場を１軸引っ張り状態と仮定し、しわと直交した方向の圧縮剛性を完全に 0 とする。解析は膜要素を使用して行われるため、分岐座屈理論に比べて解析コストが少なくて済むという利点がある一方、解析結果として得られるのはリンク*1*2

東京大学大学院航空宇宙工学専攻
宇宙航空研究開発機構宇宙科学研究本部 教授

ルの発生領域や発生方向等の面内に関する情報のみに限られるという欠点を持つ。これまでに、Miller ら 2) が Stein、Hedgepeth らの修正則 3)を実装することによって、部分的にしわが生じる膜面についての解析に張力場理論を適用している。
また、Roddeman ら 4), 5)は、変形勾配テンソルを修正することによって膜面内の応力場が張力場と一致するように修正する方法を提案している。さらに、Lu ら 6)や中篠ら 7)は Roddeman らによる手法を修正し、有限要素法の解析プログラムに実装を容易とするような手法を提案している。
分岐座屈理論に基づく手法では、しわを分岐座屈現象と捉えて分岐座屈解析を行うことによって、しわの波長や振幅などの膜面の局所的な変形についての詳細な結果を得ることが可能である。しかし、しわが生じた後の解析を行うためには分岐点の計算と座屈モードの計算が必要となり、さらに解析にシェル要素を使用するため、一般的に張力場理論に基づく手法と比べて計算コストが高いという問題がある。そのため、座屈解析を行わずに分岐経路上の平衡解を得るために、初期不整を加えた不完全系
に対して解析を行う擬似分岐座屈解析が、計算コストを減らす目
的で膜面のしわ解析に適用されている 8), 9)。
一方、しわの表現理論に基づいて定式化された方程式を求解するための数値解析法としては、主に静的な反復解法である
Newton-Raphson 法や修正 Newton 法などが使用されている。
特に Newton-Raphson 法はその局所的な収束性の高さから、非線形有限要素法の解法として広く利用されているが、解析が十分平衡解の近傍にない場合には、解析の収束性が保証されないという欠点を持つ。また、しわが生じる分岐座屈点近傍において、解析が不安定になり易いという問題がある。そこで、近年、

分岐座屈点近傍でも解析が不安定とならず、安定に解析が行える動的緩和法が、膜面のしわ解析の数値解析法として利用されるようになっている 10), 11), 12)。しかし、動的緩和法は平衡解を得るために動的解析が必要なため、静的反復法に比べて解を得る
までに非常に長く解析時間がかかるという問題がある。これまでに、動的緩和法の収束性を向上させる目的で、いくつかの手法が提案されているが 13), 14)、いずれの手法も動的な運動に収束性を依存しているため、劇的な収束性の向上は難しいのが現状である。そのため、動的緩和法の収束性を改善する手法として、動的な運動に解析の収束性が依存しないような手法が求められている。
そこで、本研究では動的な運動に解析の収束性が依存せずに、動的緩和法の収束性を向上させる新たな手法についての提案を行う。そして、提案する手法の有効性を、数値解析によって従来までの手法と比較することで確認する。

2． 膜面のしわ解析のための非線形有限要素法解析
本節では膜面のしわの数値解析を行うために必要となる非線形有限要素法の定式化についてまとめる。
2.1 シェル要素の定式化
本研究では 4 節点の MITC (Mixed Interpolation of
Tensorial Components) シェル要素 15), 16)を数値解析に使用する。MITC シェル要素において、図1 に示す初期配置における要素内の任意の点に対する位置ベクトルは
0

4

X = ∑ N k ( r1 , r2 ) 0 X k +k =1

1 3 4 kr h∑ N ( r1 , r2 ) 0 V3k
2
k =1

k

座標、N (r1，r2)は節点 k における 2 次元の形状関数、h は膜厚、そして 0 V3k は初期配置における節点 k のディレクターベクトルを表している。同様に現配置における位置ベクトルも4

k =1

1 3 4 kr h∑ N ( r1 , r2 ) t V3k2
k =1

それぞれ以下のように計算される。
G ×G
G ×G
G ×G
G1 = 2 0 3 , G 2 = 3 0 1 , G 3 = 1 0 2
V
V
V
g ×g
g ×g
g ×g
g1 = 2 t 3 , g 2 = 3 t 1 , g 3 = 1 t 2
V
V
V

(5)
(6)

ここで、 0V , tV は初期配置、現配置におけるスカラー3 重積である。
0V = ( G1 × G 2 ) ⋅ G 3(7)

V = ( g1 × g 2 ) ⋅ g 3(8)

そして、初期配置から現配置への変形勾配テンソルはる節点 k の値であることを示している。また、 ri ( i = 1, 2 ) は自然

x = ∑ N k ( r1 , r2 ) t x k +加えて、初期配置、現配置における反変基底ベクトル， G i , g i もt(1)のように表せる。ここで、左上付き添え字における 0 は初期配置の値であることを示しており、右上付き添え字k は要素内におけt

図 1 MITC シェル要素の構成t03

F = ∑ gi ⊗ G i

(9)

i =1

と表され、この時 Green-Lagrange 歪テンソルは
1 t T tt
(10)0E =0F 0F − I
2と求められる。ここで、右上付き添え字 T は転置を表している。

()

第2 Piola-Kirchhoff応力テンソルは構成則テンソル 0t C を用い(2)のように表せる。また、埋め込み座標系における初期配置と玄配置の共変基底ベクトル G i と g i が、それぞれ以下のように求
められる。
Gi =

∂ 0X
∂r i

(3)

gi =

∂tx
∂r i

(4)

て構成方程式より
t0

S = 0t C : 0t E

(11)と求められる。
MITC シェル要素では、膜厚が薄い場合に生じる Shear
Locking を回避するために、面外せん断歪テンソルの成分についての内挿関数を
1
1(12)E13 = (1 + r2 ) E13A + (1 − r2 ) E13C
2
2
1
1
(13)
E23 = (1 + r1 ) E23D + (1 − r1 ) E23B
22

のように再定義する。ここで、 E13A , E23B , E13C , E23D はそれぞれ要
素内の自然座標 ( r1 , r2 ) における点 A＝（0，1），B＝（-1，0），C＝（0，-1），D＝（1，0）における面外せん断歪テンソルの成分である。式（12，13）のように面外せん断歪テンソルの内挿関数を再定義することにより、膜厚が非常に薄い場合でも Locking
を起こさずに安定に解析を行うことが可能である 15), 16)。また、有限回転を含む解析における収束性を向上させる目的で、野口らによって提案されている有限回転テンソル 17)によるディレクター
ベクトルの更新計算を行っている。
2.2 擬似分岐座屈解析
しわは分岐座屈後に生じる面外変形と捉えられるため、しわの解析には分岐座屈解析が必要となる。しかし、分岐座屈解析には分岐点の計算や座屈モードの計算が必要となるため、一般的に解析にかかる時間が長い。そのため、分岐座屈解析を行わずに分岐経路上の解を得ることができる擬似分岐座屈解析が、しわの解析にしばしば利用される 8), 9)。
本研究における擬似分岐座屈解析を行う手順は次の通りである。まず、初期不整として各節点の面外方向の変位に対して、解析前に初期変位を与える。この時、初期不整として与える面外方向の変位としては、大きさが膜厚の 10％以下でランダムな分布を持つ値を与えている。そして、その後は初期不整を含んだ
膜面について通常の非線形有限要素法解析を実行することで、
分岐座屈解析を行うことなしに、分岐経路上の平衡解を得ることができる。

なる。そこで、本研究ではt

F − tQt

Ft≤ R _ tol ,u − t −Δt utu≤ u _ tol(16)に示すような判定条件に従って、解析が十分に収束したかどうかの判定を行っている。ここで、R_tol や u_tol は収束を決定する許容値である。これらの値は収束解が得られるまでに要する解析時間の長さと得られる収束解の解析精度とのバランスを考慮して決定される。しかし、一般的に動的緩和法を用いて十分な精度の収束解を得るためには計算時間が長くかかるという問題がある。この収束性の問題は特に、対象とする問題が大規模な場合、非常に大きな問題となる。
3.2 動的緩和法の収束性を改善するための従来までの手法
これまでに、動的緩和法の収束性の問題を改善するための手法がいくつか提案されている。まず、Zhang ら 13)によって減衰マトリクスを適応的に変化させる適応減衰法が提案されている。
適応減衰法において減衰マトリクス C は
C = αC M(17)と与えられる。ここで、 α C は臨界減衰係数である。幾何学的非線形問題においては、臨界減衰係数の値は変位に応じて変化する。そのため、解析の各ステップにおいて臨界減衰を実現するためには臨界減衰係数の値も変位に応じて変化させる必要がある。Zhang らは各ステップにおける臨界減衰係数の値を、レーリー商を利用して
1

3． 動的緩和法
3.1 定式化
動的緩和法は強い非線形性を示す問題についても安定に解析することが可能なため、しばしば静的な反復法による解析では不安定となり易い問題について適用される。動的緩和法の基本的な考えは、外力によって生じた動的な変動を収束させることによって静的な平衡解を求めるというものである。以下、手法の詳細について述べる。まず、構造物の運動方程式は
 + C t u + t Q = t F
M tu

(14)

のように記述される。ここで、 u は変位ベクトル、 M は質量マトリ
クス、 C は減衰マトリクス、 Q は内力ベクトル、そして F は外力ベクトルである。もし今、減衰マトリクス C が存在し、その値が正 の値はであれば、変位速度ベクトル u と変位加速度ベクトル u
時間経過とともに 0 へと漸近していく。そして結果的に、動的な変動の収束解として静的な平衡解と同じ解t

Q = tF(15) の値は 0 へと漸近するものの、完全が得られる。この時、 u と u
 の値が十分小さくなったとこに0とはならない。そのため、 u と u
ろで解析を打ち切り、その時の値を収束解とする作業が必要と⎛ t uT t Q ⎞ 2
αC  2 ⎜ t T t ⎟⎝ u M u⎠(18)

によって概算することで適応的に変化させている。一方
Cundall ら 14)は運動減衰法と呼ばれる手法を提案している。運動減衰法は適応減衰法とは異なり、収束性を向上させるために減衰マトリクスによる減衰効果を利用しない。手法の具体的な実行手順は次の通りである。1) 減衰マトリクス C を 0 とした無減衰系について、式(14)を計算し動的解析を行う。2) 系の運動エネ
ルギが極大となった時に、全節点が持つ変位速度ベクトルの値を 0 とする。3) 2)の系の運動エネルギが極大となった時の変位ベクトルを初期値として、1)の無減衰系についての動的解析を再度実行する。4) 式(16)の収束判定条件によって解析が収束したと判定されるまで(1)-(3)のステップを繰り返す。以上の手順を踏むことによって過渡運動における全振動モードを効率的に減衰させることができ、短い解析時間で平衡解を得ることが可能となる。
これら従来までに提案されている適応減衰法や運動減衰法は、いずれの手法も動的緩和法における根本的な問題である動的な運動の解析に収束性が依存しているため、劇的な収束性の向上は難しいのが現状である。そのため、動的な運動の解析を必要としないような改善手法が求められている。

4.2 高収束な動的緩和法の提案
先に述べたように動的緩和法の収束性を改善するためには、解付近での局所的な変動を効率良く収束させるような改善が必要である。本研究では、この改善のために静的な反復解法、特に Newton-Raphson 法（以下 NR 法）を動的緩和法に組み合わせることを考える。NR 法は解析が解付近にある場合には高い収束性（2 次収束）を有しており、解付近での局所的な変動を効率良く収束させることが可能である。そのため、この NR 法が持つ解付近での高い収束性を利用することによって、動的緩和法の収束性を改善することを考える。また、NR 法と動的緩和法
は解近傍における収束性の違いだけでなく、いくつかの正反対の特徴を有している。

-4

6 10

収束解
-4

5 10

Y方向変位 [m]4． 本研究で提案する手法
4.1 動的緩和法の収束性が低い原因についての検討
動的緩和法の収束性を向上させるような新しい手法を考えるにあたり、まず本節では動的緩和法の収束性が低い原因について検討を行う。検討を行うにあたり、以下図 2 に示すようなトグル構造物の集中荷重問題を例にとって考える。解析に使用した
入力定数を Table 1 に示す。荷重としては P=100N の場合を考える。境界条件としては両端完全固定としている。また、解析には適応減衰法を組み込んだ動的緩和法を使用した。ここで、動的緩和法における解析の収束性について考えるため、荷重が加わっている中心節点の Y 方向変位の収束解が得られるまでの収束の様子を図 3 に示す。図から収束までに要した全解析ス
テップの約 1/4 のステップが終了した時点で、既に収束解の約98％の変形が得られていることがわかる。このことから、動的緩和法の収束性が低い原因は、大局的な変形終了後精確な解を得るために局所的な変動を収束させるのに時間がかかっているためであることがわかる。このことから、動的緩和法の収束性を向上させるためには、解付近での局所的な変動を効率良く収束させるような手法の改善が望ましいことがわかる。

-4

4 10

-4

3 10

98%
-4

2 10

-4

1 10

25%

0

0 10

0

1000 2000 3000 4000 5000 6000 7000 8000

ステップ数

図 3 負荷部の Y 方向変位の収束履歴
NR 法は解析が十分に解近傍にない場合には収束が保証されないが、動的緩和法では解析が解近傍にない場合でも解析の収束は保証されている。また、NR 法では座屈を含む解析では、計算の途中で接線剛性マトリクスが特異となり、解析が不安定となり易いという問題があるが一方、動的緩和法では座屈を含む解析についても安定に解析を行うことが可能である。以上のこと
から動的緩和法と NR 法は互いに補完し合うような関係にあり、両手法をうまく組み合わせることによって、それぞれの手法が持つ利点だけを利用するような手法が考えられる。
本研究では動的緩和法と NR 法を解析の状況に応じて切り替えるような手法を考える。解析が解遠方や座屈点近傍にある場
合には、高い安定性を持つ動的緩和法を用いて解析を行い、解析が平衡解の近傍にある場合には高い収束性を持つ NR 法を適用する。このような手法を考えることによって、従来の動的緩和法が持っていた解析の安定性はそのままに、解近傍における収束性を向上させることを考えられる。提案する手法における解析の流れを図4 に示す。まず、解析の最初の時点では解析は解近傍にないと考えられるため、動的緩和法を用いて安定に解近傍と思われる程度まで解析を収束させる。そして、解析が十分に解近傍まで収束されたら、動的緩和法から NR 法へと解析を切り替える。しかし、この時、解析の初期値となる動的緩和法の収束解が NR 法で安定に解析できる程十分に解近傍の結図 2 トグル構造物の集中荷重問題
Table 1 入力定数
Input Constant

Value

H
L
a
w
Young's Modulus
Poisson Ratio
Density

19.6 mm328.57 mm
6.17 mm
19.13 mm72.32 GPa
0.31130 kg/m3

図 4 提案する手法の解析の流れ

果となっている保証はない。そのため、NR 法の解析途中で解析が不安定となっていないかどうかの確認を毎ステップごとに行う必要がある。もし、動的緩和法の収束結果が、NR 法で安定に収束できるほど十分に解近傍になく、解析が不安定となっていると判定されたら、解析手法を NR 法から動的緩和法へと再び切り替える。そして、その時に動的緩和法の収束の判定条件を
前回の収束条件よりもさらに厳しく設定することによって、より解の近傍まで動的緩和法を用いて解析を収束させるようにする。
そして、再び動的緩和法によって解近傍まで解析を収束させた
後、NR 法へと切り替えて解析を行う。この動的緩和法と NR 法の切り替えを NR 法において平衡解が得られるまで繰り返す。
以上が提案する手法における解析の一連の流れである。このような流れで解析を行うことによって、動的緩和法が持つ解析の安定性は維持したまま、NR 法が持つ高い局所収束性を利用した解析を行うことが可能となると考えられる。また、提案する手法
において、動的緩和法には特別な修正を必要としないため、従来までの手法にそのまま適用することが可能である。そのため、
3.2 節で述べたような適応減衰法や運動減衰法などの動的緩和法の収束性を向上させる改善手法もそのまま利用することが可能である。
4.3 解の近傍判定
本研究で提案する手法において、重要となる点は解の近傍判定部分と NR 法における安定性の判別部分をいかに設計するかという点である。そこで、本節では解の近傍判定部分の詳細について述べる。
従来までの動的緩和法において、解析の収束の有無は式(16)によって判定される。そのため、十分な精度の解を得るためには、求める精度に応じてR_tolやu_tolといった許容値を十分に小さな値に設定する必要がある。しかし、本研究で提案する手法においては、動的緩和法は解の近傍まで解析を収束させる
ためだけに使用するため、従来までのように許容値として小さな値を設定する必要がない。むしろ、より大きな許容値を設定して、早く解析を NR 法へと切り替えるほうが、全体の解析時間を短くする目的では望ましい。そのため、収束の判定条件としては従来までと同じ式(16)を利用しながらも、あらかじめ設定する許容値としては、必要な精度よりも大きな値を設定する。しかし、設定した許容値で収束させた結果が、NR 法で安定に解析できるのに十分な程解近傍の結果であるかどうかを、解析の実行前に判断することは困難である。そのため、もし設定した許容値では
NR 法が安定に収束するほど解近傍まで収束できていなかった
場合には、NR 法から動的緩和法への解析の切り替え時に
R _ tol ( New) = α × R _ tol (Old )u _ tol ( New) = β × u _ tol (Old )(α < 1)( β < 1)の切り替えが生じないよう、ある程度小さな値を設定する必要がある。本研究では係数α、βの値を 0.1 に設定している。
4.4 NR 法における安定性判別
NR 法において収束解が得られない原因としては、解析の初期値が解遠方にある場合と、解析の途中で接線剛性マトリクスが特異となり、解析そのものが不安定となる場合が考えられる。まず、解析の初期値が解遠方にある場合については、途中で解析が不安定とならなくても、数十回以上ステップを繰り返しても解析が収束しないということがある。この場合、解析そのものは不安定になってはいないものの、本研究で期待している NR 法の高い収束性が十分に発揮されていないため、ある一定回数
以上ステップを繰り返しても収束解が得られていない場合には、解析が不安定になっていると判定して解析を動的緩和法に切り替えるようにしている。
次に、解析の途中で接線剛性マトリクスが特異となり、解析が不安定となる場合について考える。接線剛性マトリクスが特異となる原因の 1 つとして、解析が座屈点近傍にある場合が考えられる。特に本研究で考えているような膜構造物のしわ解析においては、非常に近い範囲に多数の分岐点が存在しているため、座屈点の近傍において解析が不安定となり易い。そこで、本研
究では解析の安定性の判定を、解析が座屈点を通過するかどうかで行う。具体的には、接線剛性マトリクスの行列式の符号変化を計算することによって、解析が座屈点を通過したかどうかの判定を行い、もし解析が座屈点を通過したと判定された場合には、解析が不安定となり易い状態にあるとして、動的緩和法に解析を切り替えるようにする。この時、接線剛性マトリクスの行列式の
値は

(

)

n

det ( K ) = det LDLT = det ( D ) = ∏ Dii

(20)

i =1

のように接線剛性マトリクスの三角分解時の対角マトリクスから計算することが可能である。ここで、 L は対角項が 1 の下三角マト
リクス、 D は対角マトリクスである。さらに、本研究のように接線剛性マトリクスの行列式の符号変化のみに注目する場合には、対角マトリクスにおける負の項の数が変化したかどうかを計算するだけで良い 18）。そのため、解析の安定性判別をNR法において毎ステップ行うことによる解析時間の増加は無視できるほど小さい。
5. 数値解析本研究で提案する手法の有効性を確認するため、数値解析を通して従来までの手法との比較を行う。

(19)の式に従って許容値をより小さな値に設定するようにする。このように設定することによって、動的緩和法と NR 法の切り替えを通して、許容値を NR 法で収束解が得られる程十分に小さな値に適応的に修正することが可能である。ただし、許容値の更新時における係数α、βの値としては、何度も動的緩和法と NR 法

5.1 トグル構造物の集中荷重問題
本問題は、4.1 節で従来までの動的緩和法の収束性が低い原因について考察した時に解析した問題であり、詳細は 4.1 節で述べたとおりである。また、以下の比較において解析結果を公平に比較するために、全ての手法で同じ解析精度の結果が得られるように収束判定の許容値を設定した。具体的には、提案する手法における NR 法の収束判定の許容値と従来までの動的緩和法における許容値に同じ値（R_tol＝10-11 、u_tol＝10-11）を用いて解析を行った。また、提案する手法において、解近傍判定の許容値の初期値は R_tol＝10-5，u_tol＝10-5 に設定した。加えて、提案する手法において実行する動的緩和法の収束性を向上させるために、運動減衰法を使用した。
まず、中心荷重 P の値を 100N から 800N まで変化させた時の中心節点の Y 方向変位についての解析結果を Table 2 に示す。従来の手法によって得られた結果と比べると、提案する手法で得られた結果と従来までの手法で得られた結果で、ほぼ同じ解析精度で解が得られていることがわかる。次に、それぞれの解析において解を得るまでにかかった計算時間を Table 3 に示す。値としては提案する手法においてかかった計算時間を 1とした時に、それぞれの手法で解を得るまでにかかった計算時間の割合を示している。結果から提案する手法を用いた場合、運動減衰法を用いた場合と比べて、計算時間が 1/3 から 1/30程
度に、適応減衰法を用いた場合と比べて 1/7 から 1/80 程度になっており、提案する手法を用いることで、劇的に計算時間を短縮できていることがわかる。また、Table 2 の結果からもわかるように本問題は P=700N の時に極限点型の座屈が生じるような問題となっているが、提案する手法で使用する NR 法に特別な安定化手法を施さなくても、座屈点を通過する問題を安定に解析することが可能であった。

5.2 矩形膜の単純せん断変形問題
問題として図 5 に示すような矩形膜の単純せん断変形について考える。本問題は上辺のせん断変位量δを与えると、膜面全体にしわが生じるような問題となっている。そのため、以下の解析
では、しわの表現理論として擬似分岐座屈理論を利用して解析を行った。境界条件としては、膜面の上辺は X 方向の自由度以
外固定、底辺は全自由度完全固定、それ以外の辺は自由として解析を行った。また、提案する手法において、解近傍判定の許容値の初期値は R_tol＝10-7，u_tol＝10-7 に設定した。
まず、提案する手法を用いて解析した時の膜面外方向の変形
結果を図 6 に示す。図において、色が薄くなっている部分が紙面手前に、色が濃い部分が紙面奥に変形していることを示しており、結果から膜面に 8 本のしわ状の変形が生じていることがわかる。次に、手法による解析結果の違いを確認するために、Y
方向中央部における膜面の断面形状の比較を図 7 に示す。

Table 2 中心節点の Y 方向変位の比較
P[N]

提案する手法
[m]

100
200
300
400
500
600
700
800

5.20E-04
1.09E-03
1.74E-03
2.49E-03
3.42E-03
4.77E-03
3.33E-02
3.39E-02

運動減衰法
との誤差
1.48E-12
3.47E-12
6.26E-12
1.06E-11
1.88E-11
4.38E-11
4.48E-10
1.34E-10

図 5 矩形膜の単純せん断変形

適応減衰法
との誤差
1.26E-11
4.80E-11
1.10E-10
8.22E-11
6.03E-11
1.17E-09
5.16E-10
1.25E-10

図 6 面外方向の変形結果
0.002

Table 3 計算時間の比較（提案する手法を 1 とした場合の割合）

提案する手法

運動減衰法

適応減衰法

1
1
1
1
1
1
1
1

2.8
4.0
4.5
7.0
11.6
5.4
2.7
30.7

7.3
12.1
12.4
18.8
39.3
29.3
7.8
76.1

面外方向変位 [m]

P
[N]
100
200
300
400
500
600
700
800

0.001

0

-0.001

提案する手法
運動減衰法
適応減衰法-0.002-0.00300.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

X [m]

図 7 Y 方向中央部における断面形状

Table 4 計算時間の比較
提案する手法 運動減衰法
解析ステップ数
計算時間
[sec]
時間比

適応減衰法

13453 (11)36993

3439341573.5 (44.6)

1499.1237442.1

1

2.6

414

（）内の値は全体の内のNR法の結果を示している

結果において、膜面の両端付近で手法による違いがわずかに見られるものの、それ以外の部分では、全ての手法で断面形状は良く一致していた。次に、それぞれの解析において解を得る
までにかかった計算時間を Table 4 に示す。結果から提案する手法を用いた場合、運動減衰法を用いた場合と比べて、計算時間が約半分以下に、適応減衰法を用いた場合と比べて約 1/400になっており、提案する手法を用いることで、劇的に計算時間を短縮することができた。
6． 結論
本研究では、動的緩和法の収束性を向上させる新しい手法の提案を行った。提案した手法は、従来までの動的緩和法に静的反復解法である Newton-Raphson 法を組み合わせるもので、解の近傍判定と安定性判別の 2 つの判定条件によって、解析の状況に応じて使用する手法を切り替えるような手法となっている。
そして、提案する手法を従来までの手法と数値解析を通して比較することによって、提案する手法が従来までの動的緩和法が持つ高い安定性を維持しながら、解析に必要となる時間を劇的に低減することを確認した。今後は、より複雑なしわが生じるような膜構造物の解析について本手法を適用していき、提案した手法の有効性の確認を行っていく予定である。

参考文献
1) Hossain, N., M., A., Jenkins, C., H., Woo, K. andIgawa, H., “Transverse Vibration Analysis for Partly
Wrinkled Membranes,” Journal of Spacecraft andRockets, Vol.43, No.3, pp.626-637, 2006.
2) Miller, R., K. and Hedgepeth, J., M., “An Algorithm for
Finite Element Analysis of Partly Wrinkled
Membranes,” AIAA Journal, Vol.20, No.12,
pp.1761-1763, 1982.
3) Stein, M. and Hedgepeth, J., M., “Analysis of Partly
Wrinkled Membranes,” NASA TN D-813, 1961.
4) Roddeman, D., G., Drukker, J., Oomens, C., W., J. andJanssen, J., D., “The Wrinkling of Thin Membranes:
Part I-Theory,” Journal of Applied Mechanics, Vol.54,
pp.884-887, 1987.
5) Roddeman, D., G., Drukker, J., Oomens, C., W., J. andJanssen, J., D., “The Wrinkling of Thin Membranes:
Part II-Numerical Analysis,” Journal of Applied
Mechanics, Vol.54, pp.888-892, 1987.
6) Lu, K., Accorsi, M. and Leonard, J., “Finite Element
Analysis of Membrane Wrinkling,” International
Journal for Numerical Methods in Engineering,
Vol.50, pp.1017-1038, 2001.
7) Nakashino, K. and Natori, M., C., “Efficient
Modification Scheme of Stress-Strain Tensor for
Wrinkled Membranes,” AIAA Journal, Vol.43, No.1,
pp.206-215, 2005.
8) Wong, Y., W. and Pellegrino, S., “Computation of
Wrinkle Amplitudes in Thin Membranes,” 43rd
AIAA/ASME/ASCE/AHS/ASC Structures, Structural
Dynamics,
andMaterials
Conference,
AIAA-2002-1369, 2002.
9) Iwasa, T., Natori, M., C. and Higuchi, K., “Evaluationof Tension Field Theory for Wrinkling Analysis with
Respect to the Post-Buckling Study,” Journal ofApplied Mechanics, Vol.71, pp.532-540, 2004.
10) Haseganu, E., M. and Steigmann, D., J., “Analysis ofPartly Wrinkled Membranes by the Method of
Dynamic Relaxation,” Computational Mechanics,
Vol.14, pp.596-614, 1994.
11) Kadkhodayan, M., Zhang, L., C. and Sowerby, R.,“Analysis of Wrinkling and Buckling of Elastic Plates
by DXDR Method,” Computers & Structures, Vol.65,
No.4, pp.561-574, 1997.
12) Zhang, W., Hisada, T. and Noguchi, H., “Post-Buckling
Analysis of Shell and Membrane Structures by
Dynamic Relaxation Method,” Computational
Mechanics, Vol.26, pp.267-272, 2000.
13) Zhang, L., C. and Yu, T., X., “Modified Adaptive
Dynamic Relaxation Method and Its Application to
Elastic-Plastic Bending and Wrinkling of Circular
Plates,” Computers & Structures, Vol.33, pp.609-614,
1989.
14) Cundall, P., A., “Explicit Finite Difference Methods inGeomechanics,” Proceedings of the EF Conference on
Numerical Methods in Geomechanics, Vol.1,
pp.132-150, 1976.
15) Dvorkin, E., N. and Bathe, K., J., “A Continuum
Mechanics Based Four-Node Shell Element for
General Nonlinear Analysis,” Engineering Computers,
Vol.1, pp.77-88, 1984.
16) Bathe, K., J. and Dvorkin, E., N., “A Formulation ofGeneral Shell Elements – The Use of Mixed
Interpolation of Tensorial Components,” International
Journal for Numerical Methods in Engineering,
Vol.22, pp.697-722, 1986.

17) 野口裕久，久田俊明，“有限回転増分を考慮した効率的シ
ェル要素の開発およびその評価，” 日本機械学会論文集（A 編），Vol.55，No.550，pp.943-950，1992．
18) 久田俊明、野口裕久、“非線形有限要素法の基礎と応用，”
丸善，1995.

WRINKLING ANALISIS FOR MEMBRANE STRUCTURE USING DYNAMIC RELAXATION METHOD

Munetaka Kashiwa *1)
Junjiro Onoda *2)
SYNOPSIS
Dynamic Relaxation (DR) method has been recently applied to investigate wrinkled membrane structures. DR method is very stable
and can solve stably problems including strong nonlinearity. However, generally speaking, the calculation cost to obtain the converged
solution is very high. In this paper, a new DR method that combines DR method and static iterative solution technique is proposed. The
proposed method can improve convergence speed of DR method without detracting the analysis stability of DR method. The
performance of the proposed method is verified through a comparison of numerical analyses for wrinkled membrane with the new
method and conventional methods.

*1）Graduate Student at University of Tokyo, Department of Aeronautical and Astronautical*2）Professor at Institute of Space and Astronautical Science, JAXA