張力場膜要素で求まる膜面形状について

岩 佐
名 取
樋 口古 谷梗貴 史∗通 弘∗∗健∗∗∗寛∗∗∗∗

概張力場膜要素で求まる膜面形状の妥当性について検討するため，張力場理論に基づく非線形有限要素法を用いて膜面の形状解析を行った．本研究では，矩形膜（ Polyimid Film ）の三隅を固定し残りの一隅に強制変位を導入した際の薄膜の挙動を解析した．解析結果と実験結果の比較より，張力場膜要素で求まる膜面形状は一軸引張領域（ wrinkled region ）と二軸引張領域（ taut region ）において実際の膜面形状を最小二乗法で一次近似した仮想的な中立面と最大リンクル高の30%以下のバラツキ誤差で一致する結果となった．本論文は，張力場膜要素を用いて計算される膜面形状と実際の膜面形状との対応関係について検証した結果を報告する．

1. はじめに薄膜の挙動は，局所的な座屈で生じるリンクリング現象を含め幾何学的非線形性が非常に強い現象である．このような幾何学的非線形性の強い薄膜の形状を詳細に求めるには，分岐座屈理論に基づく数値解析を行わなくてはならない．しかし，分岐座屈理論に基づく解析は薄膜の面積が大きくなるにつれ計算コストが高くなり，ソーラーセイルやサンシールドのような大規模な薄膜構造物の挙動解析をするにあたって適用が困難となる．そのため，次世代の大型膜面宇宙構造システムを実現するには，計算コストの小さい実用的な薄膜の挙動解析法の構築が必要となってくる．
張力場理論は，膜の持つ微小な曲げ剛性を零と仮定した通常の膜理論に加え圧縮に対する抵抗力も零とした簡易理論であり，計算コストが小さいことから薄膜の解析手法として従来から広く使用されてきた1−3) ．張力場理論を用いた薄膜の解析は，静的な解析だけでなくエアーバッグ ·パラシュート ·ソーラーセイル等の動的な挙動解析にも適用され，膜面の接触問題を扱ったものも報告されている4) ．このような張力場理論による薄膜の動的解析は，薄膜の応力状態だけでなく時々刻々変化する膜面形状を解析結果として出力するため，∗宇宙航空研究開発機構 招聘研究員·博士（工学）
宇宙航空研究開発機構 教授·工博∗∗∗
宇宙航空研究開発機構 助教授·工博∗∗∗∗
東京工業大学助教授·工博∗∗

大型膜面宇宙構造システムの動的な挙動を予測する実用的な解析法として期待されている．
しかし，張力場理論による解析は膜面に局所的に生じるリンクリング形状を詳細に解析することができず，曲げ剛性と圧縮に対する抵抗力が零の仮想的な薄膜の挙動を解析しているにすぎない．そのため，張力場理論に基づく実用的な薄膜の解析法を構築するためには，実際の薄膜の挙動との対応関係を詳細に検証することが重要となる．現在のところ，張力場理論で得られる膜面の応力場をシェル理論で導かれる結果と比較した論文は報告されているものの5) ，膜面形状といった薄膜の挙動そのものに着目して張力場理論の結果を検証した論文は少ない．
本研究では，はじめに張力場理論に基づく非線形有限要素法にViscous Relaxation法を適用し，解析の収束性を向上させた膜面の形状解析法を構築した．その後，リンクルやスラックを有する矩形のPolyimid Filmを対象に膜面の形状解析を行い，得られた結果を別途実施した実験結果と比較した．
本論文では，解析結果と実験結果の比較より，張力場理論で得られる膜面形状と実際の膜面形状との対応関係について検証した結果を報告する．

2. 張力場理論による非線形有限要素解析
張力場理論は，膜の持つ曲げ剛性と圧縮に対する抵抗力をともに零と仮定した理論であり，その定式化にはいくつかの方法が提案されている1−3) ．薄膜の挙動解析は非線形性が強く，張力場理論に基づく膜要素（張力場膜要素）を用いた解析であっても定式化によって解析の安定性が異なる．本研究では，いくつか提案されている張力場理論の中でも解析の収束性が良好なRoddeman理論3)に基づく張力場膜要素を用いて解析を行った．
張力場膜要素は，膜面の初期状態を無応力状態とした場合，膜面に垂直に作用する荷重に対して抵抗力を持たない．これは，膜の曲げ剛性を零と仮定しているためである．そのため，解析初期の段階において重力のような微小荷重が膜面に垂直に作用すると，接線剛性マトリクスが特異となり解析が安定しない．この問題に対処する方法として，一般逆行列を用いる方法6) ，ペナルティーばねによる特異行列の回避法7) ，Viscous Relaxation法8)等が提案されているが，本研究ではプログラミングの簡便性ならびに解析の安定性からViscous
Relaxation法を適用した．
本解析は，汎用有限要素解析プログラムFeappv（ Finite
Element Analysis Program-personal version ）9) を用いて行った．解析で用いた張力場膜要素ならびにViscous Relaxation
法は，著者がFeappvへ別途追加した．以下，解析で用いた張力場膜要素ならびにViscous Relaxation法の定式化について示す．

2.1 張力場膜要素
Roddeman理論に基づく張力場膜要素は，リンクル領域内において変形勾配テンソルを修正することによって圧縮応力を解放する3) ．秋田は，Roddeman理論に基づく張力場膜要素の定式化を修正弾性マトリクスを用いて簡易的に行えることを示した10) ．修正弾性マトリクスは，射影行列を用いて通常の弾性マトリクスを修正したものであり，この修正によりリンクル領域内の圧縮応力を解放する．本解析では，定式化の簡便性ならびに解析の安定性を考慮し，秋田が提案した射影行列による修正弾性マトリクスを用いて張力場膜要素の定式化を行った．
Total Lagrange法における時刻tの修正弾性マトリクスt0 C̄
は，以下の式で表される．
 a2 − (1 − a1 )ν
E t
0 C̄ =a1 a2 ν1 − ν2 

0

ただし，

2a1 a2 νa1 a200 

0 

aG

(1)

a2はスラック領域内の修正パラメータで薄膜のひずみ状態に応じて以下のように設定される．
• ²2 + ν²1 ≥ 0 · · · taut region
a1 = 1, a2 = 1• ²1 ≤ 0 · · · slack region
a1 = 1, a2 ¿ 1• ²2 + ν²1 < 0, ²1 > 0 · · · wrinkled region
a1 ¿ 1, a2 = 1
a1,2をひずみ状態に関係なく 1 とした時には，通常の弾性マトリクスとなる．また，薄膜が完全に圧縮力を負担しないと仮定した場合には，リンクル領域でa1 = 0，スラック領域でa2 = 0である．
薄膜が圧縮力を完全に負担しないと仮定した解析は，広範囲にわたってスラックが生じるような問題に対して剛性マトリクスが特異となり安定しない10) ．そこで，本解析では，a1,2を以下のように設定することで，薄膜がリンクルおよびスラック領域において微小な圧縮応力を負担するように修正している．
a1 = 1 × 10−5 ：リンクル領域
a2 = 1 × 10−5 ：スラック領域
以上の修正弾性マトリクスを通常の膜要素の弾性マトリクスに置き換え，張力場膜要素の定式化を行った．

2.2 Viscous Relaxation法
張力場膜要素は膜の持つ曲げ剛性を零と仮定するため，解析初期の無応力状態において重力のような微小荷重が膜面に垂直に作用すると剛性マトリクスに特異性が生じ解析が安定しない．この問題に対処する方法として，本研究ではViscous Relaxation法を適用した．Viscous Relaxation法
は，静的なつり合い方程式に仮想的な減衰項を加えた一階の微分方程式を解くことによって，静的なつり合い解を求める準動的解析法の一つである．本方法は，解析初期の段階で特異性を生じる接線剛性マトリクスに仮想的な減衰項を加えることによって特異性を回避している．本研究では，Websterにより提案されたViscous Relaxation法8)を用いた．
一般に，静的なつり合い方程式はF (u) = P(3)で表され，uは変位ベクトル，F は内力ベクトル，P は外力ベクトルである．上式に Newton-Raphson 法を用いると，[K](t+∆t u −t u) =t+∆t P −t F(4)となる．ただし，

ν²1 + ²2 i1−νha2 + (1 − a1 )Ga =2²1 − ²2(2)である．ここで，E：ヤング率，ν：ポアソン比，²1,2：最大最小主ひずみである．a1はリンクル領域内の修正パラメータ，[K] =³ ∂F ´¯¯¯∂u t u(5)である．張力場膜要素を用いた解析は曲げ剛性を零としているため，解析初期の段階で膜面に垂直な外力が作用すると式（ 5 ）の剛性マトリクス[K]が特異となり解を求めることができなくなる．そこで，以下に示すような仮想的な減衰項を加えた一階の微分方程式を考える．

ここで，[D]は仮想的な減衰項である．上式を解くと静的なつり合い解は，速度ベクトルが零になった時に得られる．具体的には，未知量t+∆t u，t+∆t u̇の間の関係式
ttu = u + (1 − α)∆t · u̇ + α∆t ·t+∆tu̇: Polyimid Film: TAB: Fixed: Displacement

TAB−1TAB−2(6)[D]u̇ + F (u) = Pt+∆t  

yTAB−3TAB−4x(7)図 1 解析モデル
表 1 薄膜のパラメータを用いて，式（ 6 ）を[K̄] · t+∆t u = t+∆t P − t F +1−α[D] · t u̇α(8)パラメータ諸量

ヤング率, E [MPa]2959.6

ポアソン比, ν0.34密度, ρ [g/cm3 ]h¡ 1¢i[K̄] = [K] +[D]α∆t(9)と展開する．そして，式（ 8 ）を満たすように変位を求める．
αは 0〜1 までの定数であり，本解析では収束性の観点か
ら α = 1とした8) ．減衰項[D]は，単位マトリックス[I]と変数C を用いて

1.42

サイズ [mm × mm]

200 × 200

膜厚, t [µm]25強制変位量, d [µm]

722

TAB-1

Amplitude [mm][D] = C[I](10)と表し，要素減衰として与えた．ここで，変数C は52.50-2.5-5TAB-3(11)TAB-40.050.1
x-location [m]で与えられ，βは0.150.050.1y-location [m]0

図 2 初期形状¡¢
β = avr diag[K](12)である．µは，解析時の収束性を加速させるための変数で，解析時の速度ベクトルに応じて決定される．初期値は0µ = 1 × 10−5とした8) ．

3. 解析モデル

0.20.150

C =µ·βTAB-2

解析モデルを図1に示す．本解析は，矩形のPolyimid Filmを対象とした．矩形膜の３隅（ TAB-1,2,3 ）を固定し，残りの１隅（ TAB-4 ）に強制変位を与えた際の膜面形状を解析した．強制変位は，図に示す通りx軸の方向に導入している．
Polyimid Filmの解析パラメータは表1に示すとおりである．
この解析モデルは，別途実施した模型実験11,12)に対応させたものである．解析モデルの節点数は10201点とし，節点の位置を模型実験における変位の計測点に一致させている．
解析で用いた矩形膜の初期形状は，実験で計測した値を用

表 2 TABの材料特性パラメータ諸量ヤング率, E [MPa]1000ポアソン比, ν0.34

サイズ [mm2 ]10 × 10板厚, t [µm]25いた．図2に初期形状を示す．初期応力状態は無応力状態とした．載荷方法は，重力（ 1G場）をx − y平面に対し鉛直下向に作用させ（初期自重解析），その後1G場を保った状態でTAB-4に強制変位（ d = 722µm ）を導入した．
実験では，強制変位導入時に矩形膜とTABとの接着面（図1：TABの領域）で微小変形が生じ，TAB周辺の引張応力が緩和される結果となった12) ．これは，矩形膜とTABとの接着を両面テープで行ったためである．この不完全固定条件で

TAB-1

Amplitude [mm]52.50-2.5-5TAB-2TAB-3TAB-4

0.050.1
x-location [m](a)0.050.150.1
y-location [m]0.000

図 4 膜面の応力分布

実験結果

Section (y=0.02 m)

TAB-2

TAB-30.2

Vertical Disp. [mm]5.0

Amplitude [mm]TAB-4(b)

0.15

0-2.5

解析結果

図 3 載荷後の膜面形状の比較
生じる引張応力の緩和は，TABの材料定数を修正することによって簡易的に再現できることがシェル要素を用いた解析により判ったため12) ，本解析ではTABの材料定数としてこの修正したパラメータを用いた．表2に解析で用いたTABの修正後のパラメータを示す．

4. 解析結果と実験結果の比較
図3に解析結果と実験結果の比較を示す．実験結果において矩形膜辺上に計測エラーが生じているが，これはレーザー変位計の焦点が膜面から外れたために生じた結果であり11,12) ，本研究では検討の対象外としている．解析結果と実験結果の比較より，張力場膜要素で得られる膜面形状は矩形膜対角線に生じるリンクリング現象を再現できないものの，膜面全体の変形形状を定性的に再現しているのが判断できる．
図4に膜面の応力分布を示す．図に示す応力分布は，一軸引張領域（ wrinkled region，図中グレー領域），二軸引張領域（ taut region，図中ホワイト領域）および微小応力領域（ slack region，図中ブラック領域）の３領域を示している．
一軸引張領域（ wrinkled region ）および二軸引張領域（ tautregion ）は，最大主応力が0.4M P a以上の領域とし，微小応力領域（ slack region ）は最大主応力が 0.4M P a以下の領域とした．なお，引張領域と微小応力領域との境界値0.4M P aは初期自重解析で生じる最大主応力値を参考に設定したものであり，本研究では重力によるたるみ状態を二軸引張状態（ taut ）と考えず微小応力状態（ slack ）として扱っている．

00.05

0.1
0.15
x-location [m]Section (y=0.10 m)

5.0
Vertical Disp. [mm]0.1x-location [m]

0.05

0.0

0.1
y-location [m]0.2

Analysis
Experiment2.5
0.0-2.5-5.0

00.05

0.1

0.15

0.2

x-location [m]Section (y=0.16 m)
5.0
Vertical Disp. [mm]0.05

Analysis
Experiment2.5-5.0

0.15
00.05

0.10
0.15
x-location [m]

TAB-1

5
2.50-2.5
-5

0.20
0.15
0.10
0.05 y-location [m]

0.2
0.15

01
0.8
0.6
0.4
0.20Analysis
Experiment2.50.0-2.5-5.0

00.050.10.15x-location [m]0.2

図 5 断面の変位分布の比較
図より，応力分布をy軸方向に見てみると，y = 0.00m〜0.05mの領域では二軸引張領域（ taut region ）が支配的となっており，y = 0.05m〜0.10mの領域では一軸引張領域（ wrinkledregion ）が支配的となっている．そして，y = 0.10m〜0.20mの領域では，微小応力領域（ slack region ）が現れ，一軸引張領域（ wrinkled region ）と混在する状態となっている．
この応力分布は図3（ a ）に示す実験結果と概ね対応しており，張力場膜要素で得られる応力状態の妥当性が判断できる．
次に，図4に示す三つの応力状態が支配的となる断面の変形形状に着目し ，張力場膜要素で得られる膜面形状と実際の膜面形状との対応関係について詳細に検討した．図5に
y = 0.02m, 0.10m, 0.16mにおける断面の変位分布の比較をVertical Disp. [mm]taut, wrinkle

00.05

0.1
0.15
x-location [m]1512.5107.552.50-2.50.2（ a ）最大主応力分布による引張領域の推定
Section (y=0.10 m)

5.0

Experiment
Neutral Surface

2.50.0-2.5-5.000.05

0.1

0.15

0.2

x-location [m]（ b ）最小二乗法による中立面の設定

図 6 引張領域の推定と中立面の設定
Section (y=0.10 m)

5. 張力場膜要素で得られる膜面形状について

5.0
Vertical Disp. [mm]Analysis
Neutral Surface

2.50.0-2.5-5.000.05

0.1
0.15
x-location [m]0.2図 7 実際の膜面の中立面との比較
0.40.2MPa
0.4MPa1.0MPa0.3∆E[mm]

張力場膜要素で得られる膜面形状は，微小応力領域（ slackregion ）において実験結果と大きな違いが生じ るものの，二軸引張領域（ taut region ）では実際の膜面形状とよく一致し，一軸引張領域（ wrinkled region ）ではリンクル波を概ね２等分する結果となった．そこで，一軸引張領域（ wrinkled region ）と二軸引張領域（ taut region ）に着目し，張力場膜要素で得られる膜面形状と実際の膜面形状との対応関係について検討した．
本検討では，はじめに実際の膜面における応力状態を張力場膜要素で得られる最大主応力分布を用いて推定し（図6（ a ）），実際の膜面に対し一軸引張領域（ wrinkle region ）或いは二軸引張領域（ taut region ）と推定された領域の変形形状から最小二乗法を用いて中立面を計算した（図6（ b ））．その後，計算で求めた実際の膜面形状の中立面と張力場膜要素で求めた膜面形状との対応関係について検討した．一般に，初期形状を平坦とした解析では，張力場膜要素で求まる膜面形状はリンクル波の中立面（リンクル波を２等分する面）に対応する．しかし，初期形状が平坦でない場合に対し，同様の結果が成立するかどうか検討した論文は著者が知るところ皆無であり，本研究では実験結果を用いることによって実現象と張力場膜要素で求まる膜面形状との対応関係について検討した．
図7に，断面y = 0.10mにおける実際の膜面の中立面と張力場膜要素で求まる膜面形状との比較を示す．図に示す中立面は，実際の膜面形状に対し一次近似（直線近似）で求めた結果である．図より，張力場膜要素で得られる膜面形状は，一次近似で求めた中立面に対し概ね一致しているのが判る．
図8に，一次近似で求まる中立面と張力場膜要素で求まる膜面形状との差∆E を示す．同図は，各断面（ y = 0.02m〜0.18m ）で生じる誤差を示しているが，y = 0.00m, 0.20mの２断面に関しては実験結果に計測エラーが生じているため

Exp-Disp
Prin. Stress

Major Prin. Stress [MPa]Section (y=0.10 m)
5.0
2.5
0.0-2.5-5.0-7.5-10.0-12.5

Vertical Disp. [mm]示す．同図より，二軸引張領域（ taut region ）が支配的となる y = 0.02mの断面の変形形状は，解析結果と実験結果でよく一致しており，張力場膜要素で得られる膜面形状の妥当性が判断できる．また，一軸引張領域（ wrinkled region ）が支配的となるy = 0.10mの断面の変形形状は，解析結果が実際の膜面形状（リンクル波）を概ね２等分する結果となっている．そして，一軸引張領域（ wrinkled region ）と微小応力領域（ slack region ）が混在するy = 0.16mの断面の変形形状は，一軸引張領域（ wrinkled region ）と推定される領域（ x = 0.00m〜0.07m ）では 解析結果が実際の膜面形状（ リンクル波）を概ね２等分する結果となっているのに対し ，微小応力領域（ slack region ）と推定される領域（ x = 0.07m〜0.18m ）では解析結果が実験結果よりも鉛直変位が小さく再現性がよくない．この微小応力領域（ slackregion ）における膜面形状の違いは，別途実施したシェル要素による解析結果でも確認されており，TAB-4に導入した強制変位が実験において正確にx軸方向に沿っていなかったことが原因と考えている12) ．

0.2
0.1
0.00.000.050.100.15
y-location [m]0.20

図 8 各断面の誤差の比較検討対象外としている（図3参照）．なお，同図は一軸· 二軸引張領域と微小応力領域を規定する最大主応力の境界値を0.2M P a，1.0M P aとした場合の結果も同時に示してある．
図に示す誤差∆Eは，張力場膜要素で求まる膜面形状と中立面とのバラツキ誤差を示しており，断面毎に以下の式で計算している．
vu
Nx0u 1 Xt(wtf t,n − wapr,n )2(13)∆Ey =Nx0n=1

ここで，Nx0 は断面yにおける一軸引張領域（ wrinkled region ）或いは二軸引張領域（ taut region ）内の節点数，wtf t,nは張力場膜要素で計算した同領域内における節点nの面外変位，

1000.4MPaRate [%]

8060402000.00

0.050.100.15
y-location [m]0.20

図 9 微小応力状態の占める割合
Section (y=0.10 m)
Vertical Disp. [mm]5.0

Analysis
Neutral Surface

2.50.0-2.5-5.0

0

0.05

0.1

0.15

0.2

x-location [m]（ a ）二次近似

Section (y=0.10 m)

Vertical Disp. [mm]5.0Analysis
Neutral Surface

2.50.0-2.5
-5.0

0

0.05

0.1
0.15
x-location [m]0.2（ b ）三次近似

図 10 高次関数で近似した中立面との比較
100.00Third order
Second orderFirst order∆E/wmax

10.00
1.00
0.10
0.01
0.00

0.05

0.10
0.15
x-location [m]0.20

図 11 最大リンクル高に対するバラツキ誤差
Section (y=0.16 m)
5.0
Vertical Disp. [mm]waprは一次近似で求めた中立面における節点nの鉛直方向の位置をそれぞれ表している．
図に示すように，一次近似による中立面と張力場膜要素による膜面形状とのバラツキ誤差∆Eは，y = 0.00mの断面から離れるに従い大きくなる傾向にある．これは，y = 0.00mの断面から離れるに従い微小応力領域の占める割合が大きくなり，一軸·二軸引張領域の変形形状に微小応力領域（ slackregion ）の与える影響が大きくなるためと考えられる．また，一軸·二軸引張領域と微小応力領域を規定する最大主応力の境界値が大きいほどバラツキ誤差が小さくなり，張力場膜要素で得られる膜面形状が一次近似で求めた中立面に近付く．これは，最大主応力の境界値が大きいほど一軸·二軸引張領域内の平均張力が大きくなり，引張領域内の変形形状に微小応力領域（ slack region ）の与える影響が小さくなるためと考えている．
図9に，各断面における微小応力領域の占める割合を示す．
同図は，引張領域と微小応力領域の境界値を0.4M P aとした場合の結果を示している．図に示すように，微小応力領域はy = 0.00mから離れるに従い増加しており，これが図8のバラツキ誤差に影響を与えているものと考えられる．しかし，図8の0.4M P a の結果と図9を詳細に比較してみると，図8のバラツキ誤差∆Eは断面y = 0.14m〜0.20m の領域において減少傾向を示すのに対し，図9 の微小応力領域の占める割合は単調増加しており完全な対応関係が見られない．これは，張力場膜要素による膜面形状と一次近似による中立面とのバラツキ誤差∆Eが微小応力領域以外にも影響を受けることを示しており，今後シェル要素による解析結果と比較することによりバラツキ誤差に与える影響を詳細に検討していく必要がある．
図10に中立面を表す近似関数の次数を変化させた場合の結果を示す．図に示す結果は，引張領域と微小応力領域の境界値を0.4M P aとした場合である．図に示すように，近似曲線の次数が大きくなると張力場膜要素で求まる膜面形状と中立面との差は大きくなるのが判る．
図11に張力場膜要素で求まる膜面形状のバラツキ誤差∆Eを最大リンクル高wmax（中立面に対する最大リンクル高）で正規化した結果を示す．図に示す通り，近似曲線の次数が小さいほど，張力場膜要素で求まる膜面形状は中立面に近付き，そのバラツキ誤差は一次近似の場合で最大リンクル高の30%以下となっているのが判断できる．
最後に，図12に初期形状から最終形状へ移行する際の膜面形状を示す．同図は，張力場膜要素による解析結果を示している．図に示すように，初期の段階でたるんだ状態にある薄膜に荷重を載荷すると，薄膜は鉛直上向き方向へ移動し，最終形状へ収束していく．この時，一軸引張領域或いは二軸引張領域（図中x = 0.00m〜0.07m ）において，張力場膜要素で計算される膜面の最終形状は実際の膜面形状を最小二乗法で一次近似した仮想的な中立面に対し，最大リンクル高の30%以下のバラツキ誤差内で一致する結果となっている．

100%
66%
33%0%

2.50.0-2.5-5.000.05

0.1
0.15
x-location [m]図 12 薄膜の変形過程

0.26. まとめ

2) Ding,H., Yang,B., Lou,M. and Fang,H. :張力場膜要素を用いて得られる膜面形状と実際の膜面形状との対応関係について検討するため，矩形膜を対象とした模型実験の数値シミュレーションを行った．本研究で得られた結果を以下にまとめる．

A Two-Viable Pa-rameter Membrane Model for Wrinkling Analysis of MembraneStructures, AIAA 2002-1460, SDM Conference, April, 2002
3) Roddeman,D.G, Drukker,J., Oomens C.W. and Janssen, J.D,The

WrinklingofThin

Membranes:Part1-Theory,

Part2-

Numerical Analysis, Journal of Applied Mechanics, Vol. 54,
pp.884-892, 1987.

• 張力場膜要素を用いた膜面形状解析は，リンクリング現象を再現できないものの，膜面形状の全体的な傾向を再現できている．

4) Kang,S., and Im,S., Finite Element Analysis of Dynamic Re-• 張力場膜要素で計算される膜面の応力分布は，実験結果と良い対応関係にある．

5) Iwasa,T., Natori,M.C., and Higuchi,K., Evaluation of Tension

sponse of Wrinkling Membranes, Comp. Math. Appl. Mech.
Engrg., 173, 227-240, 1999

Field Theory for Wrinkling Analysis with Respect to the PostBuckling Study, Journal of Applied Mechanics, Vol. 71, Transaction of American Society of Mechanical Engineers, 2004.

• 張力場膜要素で求まる膜面形状は，一軸·二軸引張領域（ wrinkled, taut region ）において実際の膜面形状を最小二乗法で一次近似した仮想的な中立面と最大リンクル高の30%以下のバラツキ誤差で一致する結果となった．

6)

半谷裕彦, 川口健一, 形態解析, 培風館, 1991.

7) Felippa,C.A, Traversing Critical Points with Penalty Springs,Transient/Dynamic Analysis and Constitutive Laws for Engineering Materials, C2, pp.1-8, 1987.
8) Webster,R.L.,• 張力場膜要素で求まる膜面形状と仮想的な中立面とのバラツキ誤差は微小応力領域等による影響と考えられるが，実験結果との比較のみでは定量的に解明することができなかった．今後は，シェル要素による解析結果と比較することによって，詳細に検討していく必要がある．

11, pp.137-145, 1980.
9) Taylor,R.L.: FEAPA Finite Element Analysis Program Version7.3 User Manual
10)

秋田 剛, 射影行列による膜面構造のリンクリング解析に関する研究,
東京大学博士論文, 2004.

11)

岩佐貴史, 名取通弘, 樋口健, リンクルに関する実験とそのスペクトル
評価, 構造工学論文集, Vol. 50B, P235-P242, 2004.3.

12)

岩佐貴史, 名取通弘, 樋口健, 古谷寛, 矩形膜のリンクリング現象に関する実験とその数値解析, 膜構造研究論文集, No.19, 2005.

参考文献
1) Stein,M. and Hedgepeth,J.M. :Analysis of Partly Wrinkled

Membranes, NASA TN D-813, 1961.7

On the Static Analysis of Structures with

Strong Geometric Nonlinearity, Computer & Structures, Vol.

Investigation on membrane surface features obtained from nite element analysis with tension eld theory
Takashi Iwasa∗)
M. C. Natori∗∗)Ken Higuchi∗∗∗)
Hiroshi Furuya∗∗∗∗)SYNOPSIS
Validity of the nonlinear nite element analysis with tension eld theory was investigated for predicting
the membrane surface feature including wrinkles and slacks. A squared wrinkled membrane, held xedat three corners and subjected to in-plane shear force at the remaining corner, was studied. Comparing
analysis results with experimental ones, it is presented that the membrane surface features in taut andwrinkled regions given by the analysis are in roughly agreement with the virtual neutral surface given by
the least mean square approximation of the actual wrinkled membranes. It is also shown that membranesurface variation in the analysis against the virtual neutral surface is within one third of the maximum
amplitude of the wrinkles.
∗)Dr.-Eng. Invited Scientist, Japan Aerospace Exploration Agency)Dr.-Eng. Prof. Japan Aerospace Exploration Agency∗∗∗)Dr.-Eng. Assoc. Prof. Japan Aerospace Exploration Agency∗∗∗∗)Dr.-Eng. Assoc. Prof. Tokyo Institute of Technology∗∗