一定の圧縮応力を許容する張力場モデルの定式化とそのインフレータブルビームの曲げ解析への適用

秋田 剛*1
石田良平*2

1.

梗

概

本論文では，一定の圧縮応力を許容する張力場モデルの定式化を示す．定式化では，許容する圧縮応力の大きさを表す圧縮パラメータを導入し，射影行列を使って修正弾性行列を導出する．また，増分型の修正弾性行列を導出し，その主軸座標系における成分を明示する．これにテンソルの座標変換式を適用すれば，一定の圧縮応力を許容する張力場モデルの接線剛性行列を，容易に求めることが可能である．本モデルをインフレータブルビームの曲げ解析に適用し，実験との比較から有効性を検証する．

1.

はじめに近年の宇宙開発の発展に伴い，ソーラーセイル(1)や宇宙アンテナ用のサンシールド(2)など軽量で収納効率の高い膜材を利用した軽量膜面宇宙構造物が数多く提案されている．それに伴い，リンクルを生じた膜面の非線形有限要素解析法も多数提案されている(3-6)．
多くの場合，リンクル現象は張力場モデルに基づいてモデル化される．張力場モデルでは，膜面の曲げ剛性を零と仮定し，膜面に圧縮応力が作用した際，リンクルの形成によって圧縮応力は完全に解放されるものと考える（図 1(a)）．このとき膜面の応力場は，一方向の主応力が正，もう一方向の主応力が零となる 1 軸引き張りの応力（張力場）でモデル化される．また，リンクルの形成に伴う面外変形は，等方弾性膜を考えた場合，張力に直交する方向の面内収縮変形としてモデル化される（図 1(b)）．このようなモデル化により，膜面のリンクル現象を座屈解析を行うことなく平面問題として簡便に扱うことが可能となる．
張力場モデルにおいては，圧縮応力の伝達を完全に無視するため，シェル理論に基づく厳密な膜面解析の力学応答との間に差違を生じることが知られている(7)．文献 8 では，許容する圧縮応力が，矩形膜のせん断座屈後のリンクルに与える影響について考察している．また，数値解析の収束性にも影響を与えることが知られている(9)．
著者等(10)は，片持ちのインフレータブルビームの曲げ特性を対象に，リンクル発生後の圧縮応力の存在が構造特性に与*1*2

図１張力場理論の基本的な仮定える影響に着目し，張力場モデルの理論解析と実験との比較を行った．比較結果から， (1) 張力場モデルで解析した場合，曲げ剛性が実験結果に比べ低下すること，(2) 一定の圧縮応力を許容することで，理論解析と実験が良く一致すること，を確認した．本論文では，「一定の圧縮応力を許容する張力場モデル」を汎用的な有限要素法に組み込むための定式化を示す．このような定式化は文献 8 で既になされているが，応力・ひずみ関係式が複雑なテンソル式で表されており，実際のプログラムの組み込みに労力を要する．本論文では，等方弾性宇宙航空研究開発機構 情報計算工学センター 招聘研究員
大阪府立大学工学研究科 航空宇宙工学分野

膜に限定した形ではあるが，主軸座標系で表した極めて簡便な応力・ひずみ関係式を示す．また，本研究で示した定式化を，文献 10 のインフレータブルビームの非線形有限要素解析に適用し，実験結果との比較し，有効性を検証する．

2.

一定の圧縮応力を許容する張力場モデル

2.1. 等方弾性膜の基礎式
本論文では等方弾性膜を解析の対象とする．平面応力状態にある等方弾性膜の x-y 座標系における応力，ひずみ成分は次式のように三次元ベクトルで表すことができる．

σ   x y  xy  , ε   x  y  xy (1)応力・ひずみ関係式は弾性行列 C ，ヤング率 E ，ポアソン比 を用いて1 0 
E 
 10  ε  Cεσ2 1  1  
0 0

2 

れぞれのプロセスに対応したひずみベクトルに分解する．また，分解されたひずみから射影行列を用いて，リンクルを生じた膜面の応力・ひずみ関係を表す修正弾性行列を導出する．
いま，リンクルを生じて張力状態にある膜面の主応力値を 1 とする．まず，図 2(a)に示すように，最初のプロセスでは，膜面が張力 1 によって弾性変形したと考える．このとき，膜面は張力によって引張方向に，またポアソン効果によって収縮方向に変形する（弾性ひずみ）．次のプロセスでは，図2(b)に示すようにしわの発生によって膜面が収縮変形したと考える．張力場理論では，このプロセスの変形はエネルギを発生しない特別な変形モードとなり，この結果生じるひずみをリンクルひずみとよぶ．以上より，リンクルを生じた膜面の全体ひずみは，弾性ひずみとリンクルひずみに分解される．
次に，弾性ひずみとリンクルひずみをベクトル形式で表す．まず，式(1)の第 1 式より，一軸の張力状態にある膜面の応力ベク(2)となる．いま，主軸座標系から角度 回転した x-y 座標系を考えると，応力，ひずみ成分の主軸座標系から x-y 座標系への成分変換式を，以下のようにベクトル形式で表すことができる．

σ   1n1   2 n 2(3)Tn1  cos 2 sin 2 sin  cos  n 2  sin 2 cos 2  sin  cos  s1  cos 2 sin 2 2sin  cos  s 2  sin 2 cos 2 2sin  cos  Tε  1s1   2s 2(4)TT

ここに i ,  i (i  1, 2) はそれぞれ等方弾性膜の主応力，主ひずみであり 1   2 , 1   2(5)とする．最後に以下の式で定義されるベクトルn3 を導入する．

n3  n1  n 2(6)トル σ t は次式のように表される．

σ t   1  n1  0  n 2   1  n1(7)のように主軸角度 とひずみベクトルε の関係式が得られる．

2.2. 張力場モデルの修正弾性行列
ここでは，文献 6,11 に従い，通常の張力場モデルに基づいたしわの２つの形成プロセスを示し，膜面の全体ひずみをそ(8)

いま， ε w をリンクルひずみベクトルとすると，張力場理論の仮定から以下のように表現することができる．

εw   2w  s2

式(3)の第１式と式(6)の内積をとり式(2)を考慮すると

f ( )  n3T σ  n3T Cε  0

図２しわの形成プロセス(9)

ここで  2w は，非零のスカラー値である．式(8)と式(9)の内積を取り， s 2 と n1 が直交することを利用すると

ε w σ t    2 w 1  s 2 n1  0TT(10)となる．このことから， ε w に相当する変形は，エネルギを発生しない変形モード，つまり図 2(b)に対応する変形モードに対応することがわかる．いま，リンクル発生後の膜面の全体ひずみベクトルを ε と置く．このとき，弾性ひずみベクトル ε e は， ε と ε w の差分として次式のように表される．

2 E
 2w
1  2(27)(11)である．リンクル判定に混合基準(12)を適用した場合，リンクル状態で(28)
 2 w   2 1  0

ここで弾性ひずみは，純粋な弾性変形の結果生じたものなのでσ t と ε e の間に，次式のように通常の応力ひずみ関係が成り立つ．となるので，式(27)から式(26)の σ w は，二軸の圧縮応力状態に対応した応力ベクトルとなる．
式(21)の両辺に C をかけると，応力ベクトルの分解表現がεe  ε  ε w  ε   2w  s2σ t  Cε e  C  ε   2 w  s 2 (12)式(9)と式(12)を式(10)に代入して整理すると

 2 w  s 2T CεTs 2 Cs 2Q  s 2s 2 CTε e  Pε(14) s Cs T22P  I Q

(16)

(15)
(17)

ここで P, Q は以下の関係式を満たす射影行列である．

Q 2  Q (18)  CQ   CQ
T(19)PT CQ  0 (20)

式(14)，(15)，(17)から，全体ひずみの分解表現が次式のように与えられる．

ε  Pε  Qε  ε e  ε w(21)

また，式(4)で  0 として式(13)に代入すると

 2 w   2 1(22)を得る．
最後に膜面の修正弾性行列を導出する．式(15)を式(12)に代入すると膜面の応力ひずみ関係式

σ t  Cε(23)  CPCσ  σ t  σ w  P T σ  QT σPT σ t  PT CPε  σ tPT σ w  PT CQε  0  CPε  PT Cεσ t  Cε

応力ベクトルの分解表現と修正弾性行列の物理的解
釈
リンクルひずみベクトルε w に弾性行列 C を作用させ(25)Tを得る．ここに

(26)

(32)

図３を参照に式(32)の物理的解釈を示す．図３左図のように，全体ひずみは弾性ひずみとリンクルひずみに分解される．まず，通常の弾性行列により全体ひずみを応力ベクトルに変換する（図３中央図）．ここで，応力ベクトルは一軸の張力状態にある膜面の応力ベクトル σ t と二軸の圧縮応力状態に対応した σ w に分解される．次に σ に PT を作用させ， σ w を除去し σ t を抽出する（図３右図）．これらの作用をあわせて，修T正弾性行列が C  P C で表現される．

2.3.

σ w   2  1 0(31)となるので，PT は σ から σ w を除去し，σ t を抽出する射影行列となることがわかる．
いま，式(19)，(17)，(24)を参照に式(23)を次式のように変形する．

(24)を定義する．式(4)，(16)で  0 としてものを，式(25)に代入
すると(30)と表すことができる．上式は，応力ベクトルの射影行列による分解表現(11)である．ここで，式(17)～(20)を参照に，式(23),(25)に PT を作用させると

 は修正弾性行列であり，通常の弾性行列とを得る．ここで C射影行列の積として表される．

σ w  Cε w  CQε(29)のように得られる．さらに式(19),(20)を参照すると，式(29)を(13)を得る．さらに式(13)を式(9)と式(11)に代入して整理すると，次のようにリンクルひずみと弾性ひずみのベクトル形式が得られる．

ε w  Qεσ  CPε  CQε  σ t  σ w

図３修正弾性行列の物理的解釈

2.4.
一定の圧縮応力を許容した修正弾性行列
いま， ( 0) を，許容する圧縮応力を表すパラメータ（以降で圧縮パラメータとよぶ）とする．式(25)の両辺に /  2 をかけて式(26)を考慮すると 
σc CQε    1 
2(33)

 0 σ2 w(34)は，ε によらない圧縮応力ベクトルとなり， の値に応じた定数値となる．式(23)と式(33)の両辺を足して，式(24)を考慮
すると，一定の圧縮応力を許容する張力場モデルの応力・ひずみ関係式が σ t  σ c  σ  C  P  Q  ε  Cε
2 (35)と表される．ここに C  C P  Q
2 

(36)

は修正弾性行列である．
いま，式(19),(17)を考慮して，式(35)の修正弾性行列をC   P T  QT  C
2



 2w

 2w0

0 
E 

2   1 0(38)のように表される．

が得られる．ここで

σc 
 v2E


 2w
 C  2w
 0

(37)と表し，式(37)の物理的解釈を，図４を参照に説明する．図３と同様に，全体ひずみは弾性ひずみとリンクルひずみに分解され，全体ひずみに弾性行列を作用させ，応力ベクトルに変換する（図４左図，中央図）．次のステップでは，まず σ にPT を作用させ σ t を抽出する．同時に QT を作用させ σ w を抽出し，さらに圧縮応力を一定にするために係数 /  2 をかけてスケーリングする（図４右図）．これらのベクトルを合わせて，一定の圧縮応力を許容したリンクル状態の応力ベクトルσ が得られ，修正弾性行列が式(37)のように表される．

2.5. 増分型の応力・ひずみ関係式
式(35)に示される応力・ひずみ関係式は，C が主軸角度 の関数となるので，非線形式になる．一般に，非線形の応力・ひずみ関係式を持つ物質の有限要素解析を行う際には，線形化された増分型の応力・ひずみ関係式を求める必要がある．
本節では，増分型の修正弾性行列を導出する．
式(35)の両辺の時間微分をとると，

    2σ  Cε  C  P  Q
 ε  2 CQε2 2(39)となる．上式の右辺第 2 項，第 3 項をそれぞれ ε と行列の積の形に変形することで，増分型の修正弾性行列を導出することができる．

(A) 右辺第 2 項の変形式(16)より，射影行列 Q の時間微分は Q2 1  2 E

 n s  s n  C
T3 22T3(40)と表される．ただし上式を等出する際，以下の関係式を利用した．

sT2 Cs 2 E1  2(41)s 2  2n3(42)

式(40)の右から ε をかけ，式(7)を考慮すると

 
Qε2 1  2  sT2 Cε En3(43)が得られる．いま， f ( )  0 の時間微分から

 f ε   f  ε (44)が求められる．式(7)から，上式の分子と分母の項はそれぞれ

図４ 一定の圧縮応力を許容した場合の修正弾性行列
最後に，式(37)の射影行列に  0 を代入して整理すると，主軸座標系における修正弾性行列の成分がf nT3 Cε(45)

f
f n3T   s 2  s1  Cε n3 (46)となる．ここで，式(46)の導出に n3   s 2  s1 なる関係式を用いた．式(45)，(46)を式(44)に代入すると

 nT3 C

 s2  s1  CεT

最後に式(56)の右辺の各項に  0 を代入して整理すると，主軸座標系における増分型の修正弾性行列の成分が
E 0Cinc   0 0

0 0

(47)

εが得られる．式(47)を式(43)に代入すると

    n nT CεQε3 3(48)となる．ここに

2 1  2 EsT2 Cε s2  s1  CεT(49)である．式(17)，(48)から式(39)の右辺第 2 項は   C  P  Q
 ε   1   CQε2  2    1   Cn3nT3 Cε 2 (50)と表される．

(B) 右辺第 3 項の変形式(27)，(13)，(41)より

 2  sT2 Cε  sT2 Cε(51)となる．ここで，式(7)，(42)より上式の右辺第 2 項は 0 となることを考慮して，式(51)を式(39)の右辺第 3 項に代入すると

 2CQε   2  CQε   sT2 C  ε222(52)を得る．ここで，式(9)，(14)より

CQε  Cε w   2 wCs 2(53)となるが，これを式(52)に代入して式(27)を考慮すると

 2 1  2CQεCs 2sT2 Cε 222 E(54)を得る．
(C) 増分型の修正弾性行列式(50)，(54)を式(39)に代入すると，増分型の応力・ひずみ関係式がσ  Cinc ε(55)として得られる． ここに，  1  2 T Cinc  C  C 1   n3nT3 s s C
2 E 2 2   2 
は増分型の修正弾性行列である．

(56)


0
0

E1     1 2  1   2  (57)のように表される．
本モデルで非線形有限解析を行う際には，まず主軸座標系で式(57)を評価し，次に解析に用いる座標系にテンソルの座標変換式を用いて変換し，構成則テンソルを求めて接線剛性行列を計算すればよい．詳しい計算手順は文献 13 を参照してほしい．

3.

インフレータブルビームの曲げ解析への適用

3.1. 解析モデル
解析対象のインフレータブルビームは，図 5 に示すような全長 l ，半径 r の片持ちの円筒チューブとする．円筒チュー
ブは周辺部分の薄膜と先端部分の円盤により構成され，空気圧 p によって剛性を与えられる．
先端の円盤にせん断荷重 Fが作用すると，荷重が小さい領域では円筒チューブは片持ちはりとして機能する．荷重が大きくなると固定点近辺の下腹部の膜が圧縮応力に耐えられなくなり，しわを生じて崩壊する．
著者等(10)は本モデルに対し，通常の張力場モデルによる理論解析を行い、実験結果と比較した．その結果，解析ではインフレータブルビームの剛性が過小評価され崩壊荷重が実験値よりも低くなることを確認した．また，微小な一定の圧縮応力を許容する張力場モデルの理論解を導出し，許容する圧縮応力の大きさを調整することで，理論解と実験とがよく一致することを確認した．
ここでは，第 2 章で導出した圧縮パラメータを調整して文献 10 と同様の検証を行い，より汎用的な有限要素解析で一定の圧縮応力を許容した場合の有効性を示す．
3.2. 解析条件
インフレータブルビームを構成する膜材の物性値は，ヤン
グ率 1 GPa，ポアソン比 0.3 とする．またビームの全長を 0.6ｍ，半径を 0.03ｍ，膜厚を 60μｍとする．検証に用いる解析ケースは， p  1600, 2100,3000, 4000 Pa の 4 ケースとする．解析では，荷重増分法によりビーム先端にせん断荷重を与えて，変位・荷重曲線を求める．ケースごとに圧縮パラメータαを 6.0  105 Paから 5.0  104 Paの範囲でスイープさせ，変位・荷重曲線の実験との RMS 誤差をもっとも小さくするαの値を探索する．なお，荷重増分法は文献 13 に基づき，埋め込み座標を利用した幾何学的非線形有限要素法を用いて行った．

3.3. 実験結果
ケースごとに，文献 10 の実験値との RMS 誤差を最小とするαを求めた結果が図 6 である．図の横軸がインフレータ
ブルビームの内圧，縦軸が圧縮パラメータを示す．縦軸の値を正にするために，－αをプロットしている．図より内圧の増加に対して－αが単調に増加することがわかる．図 7～図10 に，各ケースで，RMS 誤差を最小とするαを用いて求めたインフレ－タブルビームの荷重・変位曲線を示す．図中に実験結果（Experiment data）および圧縮を許容しない通常の張力場理論の結果（Tension field theory）もあわせて示す．
52x 10 [Pa]

1.5

10.5
1500

2000

図 5 インフレ－タブルビームモデル

2500
3000
Pressure [Pa]

3500

4000

図 6 内圧とαの関係図
35303025

Applied force [gf]Applied force [gf]25201510201510500

Modified tension field theroy
Tension field theory
Experimental data
24681012
Tip displacement [mm]14Modified tension field theroy
Tension field theory
Experimental data50016

図 7 実験結果と解析結果（p = 1600 Pa）

51015
Tip displacement [mm]20

図 8 実験結果と解析結果（p = 2100 Pa）

507045604050
Applied force [gf]Applied force [gf]3530252015302010

Modified tension field theroy
Tension field theory
Experimental data

5004051015
Tip displacement [mm]20

図 9 実験結果と解析結果（p = 3000 Pa）

Modified tension field theroy
Tension field theory
Experimental data

100051015
Tip displacement [mm]2025

図 10 実験結果と解析結果（p = 4000 Pa）

図の横軸が，インフレータブルビームの先端変位，縦軸がせん断荷重である．図より，どのケースでも微小な一定の圧縮応力を許容することで，リンクルを生じる前の線形領域からリンクル発生後の非線形領域まで，解析結果と実験値がよく一致することがわかる．一方で，通常の張力場理論による結果は，線形領域から非線形領域に移行する際の荷重が，実験結果よりも小さくなることがわかる．

結論と今後の課題
本研究では，一定の圧縮応力を許容する張力場モデルを汎用的な有限要素法に組み込むために必要な，修正弾性行列の定式化を示した．定式化では圧縮パラメータを導入し，修正弾性行列の力学的な解釈を，射影行列を用いて示した．
また，主軸座標系における増分型の修正弾性行列を明示した．これにより，一般的な非線形有限要素法への組み込みが容易となる．本研究で示した解析方法を，インフレータブルビームの曲げ解析に適用し，実験との比較からその有効性を示した．
今後の課題として，圧縮パラメータの合理的な設定方法の構築があげられる．このために，実験面からのアプローチと理論面からのアプローチが考えられる．前者に関しては，実験ケースを増やし，応答曲面を利用した圧縮パラメータの設定法を検討する予定である．後者に関しては，より現象に忠実なシェル理論に基づく解析結果と本論文の解析結果を比較し，圧縮パラメータの設定法に関して考察する予定である．
4.

参考文献
1) 奥泉信克: 多粒子系モデルによるソーラーセイル膜面の展開解析, JSASS-2009-4320, 2009.
2) Ewing, A. P., Back, J. M., Schuettpelz, B. M., andLaue, G. P.,: James Webb Space Telescope Sunshield
Membrane Assembly, AIAA Paper 2009-2156, 2009.
3) Miller, R. K., Hedgepeth, J. M., Weingarten, VI., Das,
P. and Kahyai, S.: Finite Element Analysis of Partly
Wrinkled Membranes, Computers and Structures,Vol. 20, pp.631-639, 1985.
4) Nakashino, K. and Natori, M. C.: Efficient
Modification Scheme of Stress-Strain Tensor for
Wrinkled Membranes, AIAA Journal, Vol. 43,
pp.206-215, 2005.
5) Miyazaki, Y.: Wrinkle/Slack Model and Finite
Element Dynamics of Membrane, International
Journal for Numerical Methods in Engineering, Vol.
66, pp. 1179-1209, 2006.
6) Akita, T., Nakashino, K., Natori, M. C., and Park, K.
C.: A Simple Computer Implementation of Membrane
Wrinkle Behaviour via a Projection Technique,International Journal for Numerical Methods inEngineering, Vol. 71, pp.1231–1259, 2007.
7) Iwasa,T., Natori,M.C., and Higuchi,K.,.: Evaluation ofTension Field Theory for Wrinkling Analysis withRespect to the Post-Buckling Study, Journal ofApplied Mechanics, Vol. 71, Transaction of American
Society of Mechanical Engineers, 2004.
8) Inoue, S.,: Prediction Methods of Wrinkling inThin-Membrane,
Proceedings
ofthe27thInternational Symposium on Space Technology andScience, 2009-c-20s, 2009.
9) Jarasjarungkiat, A., Wuchner, R., and Bletzinger, B.
U., : Efficient Sub-Grid Scale Modeling of Membrane
Wrinkling by a Projection Method, Computer
Methods in Applied Mechanics and Engineering, Vol.
198, pp. 1097-1116, 2009.
10) 菅沼和敬, 秋田剛, 石田良平: インフレータブルビームの曲げ特性について, Vol. 75, pp. 388-393, 2009.
11) Akita, T., and Natori, M. C.,: Sensitivity Analysis
Method for Membrane Wrinkling Based on theTension-Field Theory, AIAA Journal, Vol. 46, pp.
1516-1527, 2008.
12) Kang, S., and Im, S.,: Finite Element Analysis ofWrinkling Membranes, Journal of Applied Mechanics,Vol. 64, No. 2, Transaction of American Society of
Mechanical Engineers , pp. 263-269, 1997.
13) 秋田剛,: 射影行列による膜面構造のリンクリング解析に関する研究, 東京大学博士論文, 2004.

A FORMULATION OF TENSION FIELD MODEL WITH CONSTANT COMPRESSIVE STRESSESANDITS APPLICATION TO BENDING ANALYSIS OF INFLATABLE BEAMS
Takeshi AKITA *1)
Ryohei ISHIDA *2)

SYNOPSISA formulation of a modified tension field model, which can allow small constant compressive stresses in wrinkled membranes, ispresented. The lower bound of compressive stresses allowed is controlled by a constant parameter, the choice of which has a great
influence on mechanical behaviors of wrinkled membranes. Both a modified stress-strain relation and its incremental form, which isneeded for nonlinear finite element analysis, are derived by using projection approach in explicit form. A buckling behavior of an
inflatable cylindrical beam is investigated by using the presented tension field model to show the effectiveness of the presentedformulation.
*1） Researcher, JAXA’s Engineering Digital Innovation Center (JEDI), Japan Aerospace Exploration Agency (JAXA).
*2） Osaka Prefecture University