座標仮定有限要素法を用いた張力構造の動的解析
本間 俊雄 *1
佐伯 裕介*2
坂中 玲子*2
梗

概

著者らは安定形態の座標値を直接未知量とした、ケーブルや膜材で構成する張力構造の幾何学的非線形解析手法を座標仮定有限要素法により提案し、裁断図解析や動的解析例を示してきた。解析手順は、平面上(自然状態)のケーブルや膜材を基準に、３次元空間上で構成する張力構造の応力･変形解析を可能としている。本論文は、これらの内容を発展させ、通常、座標仮定の定式化ができない梁要素を導入した張力構
造(複合構造)の動的解析例を示す。その際、膜材の裁断図解析を繰り返すことで、比較的荒い要素分割による裁断線から細かい要素分割による裁断線を決定する。動的解析では、細かい要素分割による裁断線を利用し、風荷重を想定した懸垂物を有する複合構造の挙動を調べる。

１．はじめに
著者らは、座標値を直接未知量とする幾何学的非線形項を考慮したケーブル･トラス要素の有限要素法による離散化式を示し、従来の変位仮定による構造解析と同水準であることを示している[1]。この定式化を座標仮定と呼び、膜シート(織布)に拡張して、張力構造解析一般に利用できることを明らかにした[2]。定式化と解析例では、張力構造のように想定形状
に対し、平面上の無応力の自然状態を基準に、初期張力導入時の設計曲面形状、静的あるいは動的へと連続した解析を行っている[3]。計算手順の特徴は、要素座標系から全体系に移行する際、座標変換が一切いらない点にあり、この解析法の利点を応用した曲面を滑らかに被う膜シートの立体裁断図解
析を示した[3,4]。裁断図解析では、応力指定･コンプライアン
ス最小化･形状指定の目的関数を設定した数値例を提示している。解析モデルは、曲率が小さなものと大きなものを採用
し、モデルの違いによる裁断図の差を考察した。また、空気膜構造やケーブルを導入した膜構造の裁断図解析にも応用している。時刻歴応答解析では、裁断図解析で得られたケーブ
ル補強膜構造の形態を利用し、非抗圧材に圧縮力が作用する
際の具体的な処理方法を整理した。さらにこれらの解析手順
による張力構造における設計支援システムのコンセプトをまとめている[4]。
本論文では、通常、座標仮定の定式化ができない梁要素を導入した張力構造(複合構造)の時刻歴応答解析例を示す。その際、膜材の裁断図解析を繰り返すことで、比較的荒い要素
分割を出発点とした裁断線の利用により、細かい要素分割時の裁断線決定手順を説明する。解析モデルは、ワイヤーで吊り下げられた懸垂物を有する複合構造であり、変動風速デー*1*2

鹿児島大学 工学部
建築学科
教授
工学博士
鹿児島大学 理工学研究科 建築学専攻 大学院学生

初期張力
境界条件
安定形態

不安定形態

立体裁断

+

膜材の自然状態

溶着･縫合＋接合

外力･境界条件導入

骨組構造

図１ 膜(張力)構造の解析概念図
タによる風荷重を作用させる。
２．解析内容について
ここでは、座標仮定有限要素法を用いた張力構造の解析内容･解析モデルを説明し、解析範囲を整理する。
図１に示すように、張力材であるケーブル要素と膜要素は、境界条件を設定していない平面上に無応力の自然状態で配置された状況を考える。次に、一般構造材であるトラス要素･梁要素･シェル要素･連続体要素などは、通常の有限要素解析
で行われるように境界条件を与えた初期構造形態を想定する。
この初期構造形態に、自然状態の張力材に張力を導入し、幾何剛性の付加後の全体構造を構成する。つまり、解析上のイ
ンプットデータは、材料諸情報を除けば、平面(x,y)上の自然状態の位置情報のみの張力材情報と、通常の安定形態にある
一般構造材を含む初期構造形態(座標値(x,y,z)＋境界条件)情報である。
座標仮定有限要素法による張力構造の形状解析では、境界条件を与えた状態で、張力材の幾何剛性項のみを用いた線形解析が可能である。得られる数値結果は直接構造形状の座標値となる。膜材の張力比を与えた異方張力の線形形状解析にも対応できる。得られた線形解析結果を初期値に、平衡方程
式に弾性係数を微小値に設定するとそのまま非線形形状解析に移行することができる。
応力･変形解析では、材料を線形弾性体とし、初期構造形態の位置情報を固定した状態で解析を行うと線形解析が可能である。幾何学的非線形性のみを考える場合には、total-Lagrange
法で対応できる｡しかし、材料非線形性(塑性･クリープ)に対応する場合には、update-Lagrange 法で扱わなければならない。
ただし、update-Lagrange 法においても、解析の特徴である座
標変換が一切要らない特性は保持される。ただし、トラス要素を除き、一般構造部材に対しては通常の座標変換処理が必要である。
以上、座標仮定有限要素法は、張力材や一般構造材で構成される複合構造の解析において、形状解析、裁断図解析、静的･動的解析が連続して行える。

y

Le

O

z, Z

s
2 (s2)1 (s1)xy, Y

1 (X1, Y1, Z1)
2 (X2, Y2, Z2)
1
2

O

安定形態

x, X 自然状態ケーブル要素

図２ ケーブル要素の自然状態と安定形態の座標関係

y3 (x3, y3)
x’

Se

φ x
O

Z

2 (x2, y2)

1 (x1, y1)

1 (X1, Y1, Z1)
2 (X2, Y2, Z2)
3 (X3, Y3, Z3)

x
1

y, Y

3

2

安定形態

O
自然状態膜要素x, X

３．基礎式と離散化式
3.1 基礎式
仮想仕事の原理より、応力τとひずみγを変形後の位置ベク

図３ 膜要素の自然状態と安定形態の座標関係

トル X で表現した場合の平衡方程式 F ( X, f , λ ) = 0 は次のよ

KG ( X) =∫ Ω δγ ( X ) τ ( X ) d Ω − δ X λ f = 0
T

T

(1)

ただし、f：荷重モードベクトル, λ：荷重パラメータ, δγ：ひ
ずみ増分である。
3.2 張力要素
ひずみ-座標関係は次の Green ひずみを採用する。
γ = ⎡⎣ γ ij ⎤⎦ , γ ij =*T ∂τ(6a,b)

動的解析では、次の非線形運動方程式を考える。

うに与えられる[2]。
δ XT F ( X , f , λ ) =∂B *T∫ Ω ∂X τ d Ω , K S ( X ) = ∫ Ω B ∂X d Ω

⎞
X Xβ
1⎛
− δij ⎟
⎜ δαβ α⎟xi x j2 ⎜⎝⎠( i, j, α, β = 1, 2, 3)(2)

ここで、xi ( ≡ x )：変形前の直交直線座標, Xi ( ≡ X )：変形後の直交直線座標, δij：Kronecker’s Delta である。全体座標系における変位 u と座標 x, X との関係は次式の通りとなる。
ui = X i − xi ( u = X − x )(3)

ひずみとひずみ増分関係は次式のように表現できる。
γ = B ( X ) X + C (C:定数ベクトル),

δγ = B * ( X ) δX

(4a,b)

式(1)は、従来の変位ではなく、変形後の座標値を未知量とする有限要素法による離散化のための基礎式である。なお、離散化式を解くための接線剛性行列 K t は、X を節点位置(座標値)ベクトルとすると次式で与えられる。
∂F= KT ( X ) = K G ( X ) + K S ( X )∂X(5)

ここで、KG , KS は各々幾何剛性行列と線形＋大変位剛性行列であり、変位仮定と同形式の次式で与えられる。

()

 + CX
 + F X, f , λ * = S ( t )MX(7)

ここで、M：質量行列, C：減衰行列, S：時間 t に依存した動的荷重ベクトル, F：式(1)で定義された静的な平衡方程式, λ*: λの設定値 である。
以上の基本式に対して、ケーブル要素と三角形要素の離散化式は以下の通りである。
○ケーブル(トラス一次)要素
ケーブル材はトラス一次要素を採用する(図２参照)。最終的な代表ケーブル要素の離散化平衡方程式は次式を得る。
⎡ 1⎤
Ee Ae T
G Xe ⎢
XeT GXe + Ce ⎥ − λfe = 0
2
Le⎢⎣ 2 Le
⎥⎦

(

)

(

)

(8)

ただし、ケーブル要素 e の Ee：弾性係数, Ae：断面積(要素内一定)とする。接線剛性行列 k Te は、次式で与えられる。
k Te = k Ge + k Se =Ae TAEG τ e + e e G T Xe XeT G
Le
Le3

(9)

なお、各行列･ベクトル成分は以下の通りである。
γe =

(
2L
1

e

1
Ce = −
2

2)

XeT GXe + Ce
T,

⎡H 0 0 ⎤⎢⎥G = ⎢0 H 0⎥⎢⎥⎢⎣ 0 0 H ⎥⎦

⎡1

−1⎤

, H = ⎢⎢−1 1 ⎥⎥ ,⎣⎦, Xe = ⎣⎢ X1 X 2 Y1 Y2 Z1 Z 2 ⎦⎥○三角形膜要素
最終的な代表膜要素の離散化平衡方程式は次式を得る(図３参照)。
⎡ 1⎤
1
QT Xe De ⎢XeT QXe + Ce ⎥ − λfe = 0
2
4 Se⎣⎢ 8Se
⎦⎥

(

)

(

)

(

) (

1 T
1
Q τe +QT Xe De XeT Q
4 Se
16 Se3

i’

)

(11)

i

(12a,b)

y

wi
ui

vi

θ ezi

j’

wj

Le
uj

ここで、第１項が幾何剛性行列、第２項が線形+大変形剛性行列である。各行列成分は以下の通りである。
T1
Q = ⎡⎣bˆ aˆ cˆ ⎤⎦ , Ce = ⎢⎣ −1 −1 0⎥⎦T2

θ zj

Le + ue(10)

ただし、Se は要素面積である。膜要素の接線剛性行列 k Te は次の通りに求められる。
k Te =

y

z

j

i’ Rz
i
z

θ yj

vj

i’ R y

x

i

x

j
Ry

j’

θ eyj

θ yi

x-z 投影面

j’
θ ezj

θ zi

x-y 投影面

θ eyi

Rz

j

x

図４ 剛体変形除去法における梁要素の座標系aˆ , bˆ , cˆ は次の内容である。
⎡b 0 0 ⎤⎡a 0 0 ⎤⎡ c 0 0 ⎤⎢⎥aˆ = ⎢⎢ 0 a 0 ⎥⎥ , bˆ = ⎢ 0 b 0 ⎥ , cˆ = ⎢⎢ 0 c 0 ⎥⎥⎢ 0 0 b ⎥⎢⎣ 0 0 a ⎥⎦⎢⎣ 0 0 c ⎥⎦⎣⎦(13a-c)

ただし、 a = [ a1 a2 a3 ] , b = [b1 b2 b3 ] とおくと、 a = aT a ,b = bT b , c = aT b + bT a で与えられる。

De は構成関係行列を表し、膜部材を直交異方性材と仮定す
ると次式で与えられる。
D e = R eT D e + R e⎡ E1t⎢1 − υ υ
1 2⎢⎢ υ Et+De = ⎢ 2 1
1 − υ1υ2
⎢
⎢ 0
⎢
⎣

υ1E2t
1 − υυ
1 2，E t
2

1 − υυ
1 2
0

⎤
0 ⎥
⎥
⎥
0 ⎥
⎥
G12t ⎥⎥⎦(14a-c)

各々繊維方向 1, 2 の縦弾性係数、ν1 , ν2：各々繊維方向 1, 2 に対するポアソン比, G12：せん断剛性 である。
3.3 一般構造材
一般構造材として、梁要素を取り上げる。曲げ力が生じる要素では、座標仮定による定式化ができないため、変位仮定
有限要素技術による定式化を示す。ここでは、部材の変形に追従する移動座標系を設定する剛体変位除去法を用いる[5]。
しかし、三次元問題において変形後の要素節点を結ぶ軸線方
向とこれに直交する二方向に移動座標系を設定し、変位から剛体変形成分を完全に除去した場合、剛性行列の表現が煩雑になる。本論文では、変位から剛体変形成分の一部を取り除く近似的な定式化を採用する[6]。
○移動座標系
初期状態からの外力及び変形の移動状態をパラメータ t で表現する。いま、部材 i - j が外力を受け、t = t* において i’- j’に変形したものとする(図４参照)。部材に設定した局所座標
o-xyz における要素節点変位ベクトルを de とする。de から軸方向に対しては剛体変位成分を、回転角に対しては部材角を除去して定義される移動座標系における要素節点変位ベクトルの成分d e を、次のように与える。
wi θ xi θ yi θ zi

ujd eT = ⎢⎣ue θex θeyi θ ezi θ eyj

θ ezj ⎥⎦viue =( u j − ui + Le ) + ( v j − vi ) + ( w j − wi ) − Le ,222

θexi = θ xj − θ xi , θeyi = θ yi − R y , θeyj = θ yj − R y ,θezi = θ zi − Rz , θezj = θ zj − Rz

(16a-f)

ここに部材角 Ry, Rz は次式の通りである。
⎡ cos2 φ− sin φ cosφ ⎤sin2 φ⎢⎥2
Re = ⎢ sin φcos2 φsin φ cosφ ⎥22⎢2sin φ cosφ −2sinφ cosφ cos φ − sin φ ⎥⎣⎦φ：繊維 1-方向(x’-方向)と平面座標 x とのなす角度、E1, E2：

deT = ⎢⎣ui

幾何学的な関係から、これら二つの変位ベクトルには次の関
係がある。

vj

wj

θ xj θ yj θ zj ⎥⎦(15a,b)⎛ w j − wi ⎞ ,⎛ v j − vi ⎞R y = − tan −1 ⎜⎟ Rz = tan −1 ⎜⎟⎜ u j − ui + Le ⎟⎜ u j − ui + Le ⎟⎝⎠⎝⎠(17a,b)式(16)は、t によって変化する移動座標系における要素節点変位の関係である。式(16)を t で微分すると節点変位増分ベクトルの座標変換式が得られる。
d e = Te d e

(18)

ここで、 a = da / dt を意味する。
次に節点力の移動座標系への変換を定義する。局所座標o-xyz で測った節点力ベクトル feを、移動座標系における節点
力ベクトル fe との関係は次式が成立する。
fe = TeT fe(19)

式(19)を変形し、両辺を t で微分すると次式を得る。
  T fe = TeT fe + T
e fe

(20)

fe , fe の成分は式(15)の変位成分に対応させたものとする。

また、
式(20)の右辺第２項は次のように表現することができる。
 Tf = kTe e
G 2e d e(21)

ここで、kG2e：剛体変位したことにより生じる幾何剛性行列
である。
○離散化式
材軸上の任意位置における変位は、移動座標系による要素
節点変位を用いて、軸線方向一次、法線方向三次の多項式により近似する。ひずみは式(2)に対応する Green ひずみを用いる。部材断面力とひずみの関係は次式で表される。

N = EAε , M x = GJ κ x , M y = EI yκ y , M z = EI z κ z(22a-d)

ここで、N：部材軸力、Ms：s 軸周りの部材モーメント、ε：
部材軸ひずみ、κs：s 軸周りの部材曲げひずみ、E：弾性係数、
A：断面積、Is：s 軸周りの断面二次モーメントである。これ
らの式を式(1)に対応する通常の変位仮定における仮想仕事式に導入すると、次式のような移動座標系における代表要素 eの離散化式を得る。
k e (d e ) d e − λ fe = 0(23)

ここで、 k ：部材剛性行列である。
式(23)の増分型の平衡方程e式は次式で与えられる。
k Te (d e ) d e − λ fe = 0 ,k Te = k Ge + k Se(24a,b)式(24b)は式(5)に対応する。
移動座標系から全体座標系へ変換すると、最終的な系全体の代表要素 e における平衡方程式及び接線剛性行列が次式のように得られる。
Fe ( d 0 e , f0 e , λ ) = L eT TeT k e H e (L e d 0 e ) − λ f 0e = 0()

k 0Te = L eT TeT k Te Te + k G 2 e L e

(25)
(26)

ここで、d 0 e：全体座標系の要素節点変位ベクトル、f 0 e：全体
座標系の要素節点力ベクトル、L e：固定局所座標と全体座標間の変換行列、Te：移動座標と固定局所座標間の変換行列、
H e：移動座標系への変換式、k0Te：系全体の接線剛性行列 である。
４．張力材要素と一般構造要素との関係
通常、張力構造は張力構造材のみで構成されることはなく、問題によっては複合構造として捉える必要がある。
ここでは、張力部材の座標値による離散化式(8)-(11)と一般構造部材である梁要素の変位仮定による離散化式(25),(26)を同時に扱う方法を示す。
対象構造物における変形前･後の全体座標x eとX eの関係は、局所座標系要素節点変位ベクトル de を用いることにより次式で与えられる。
(27)

L e xe + d e = L e Xe

式(27)を式(25)に代入すると次式を得る。
Fe ( X e , f0 e , λ ) = L eT TeT k e H e ( L e ( X e − x e ) ) − λ f0 e = 0(28)

同様に式(26)も同形式の式を得る。以上により、梁の代表要素
における安定形態の座標値で表現された平衡方程式(28)及び
接線剛性行列は、座標仮定によるケーブル･膜の離散化式と同等に扱える。したがって、系全体の平衡方程式 F ( X, f , λ ) = 0と接線剛性行列 KT は次式で与えられる。
F ( X, f , λ ) =m1+ m 2Ae =1
m1+ m 2Fe ( X e , f0 e , λ ) = 0

∂F= K T = A k Te ( X e )
∂X
e =1

(29)
(30)

ここで、m1：張力材要素の数、m2：安定形態の位置を設定した構造材要素の数 である。 A, A はベクトルと行列に対応させたアセンブル記号である。
式(29),(30)は、座標値を直接未知量とした基礎式になっている。従って、運動方程式も結果的に式(7)と同形式で与えられる。非線形運動方程式の解法は文献[7]に従う。
５．裁断図解析5.1 最適化問題の定式化
最適化では、次の応力指定の釣合形状と設計者の要求する想定形状との偏差量を最小化する膜シートの裁断図決定の形
状指定問題を扱う。
OP

minimizen( ) ∑ 12 ( X i − Xˆ i )T ( X i − Xˆ i )(31)

σ L ≤ σ ≤ σU

(32)

f x* =

i =1

subject to(e = 1, 2,3," , m)

ここで、n：最適化対象となる総自由度数, x*：自然状態の膜シート上の設計変数(座標値)ベクトル、X i ：計算で得られる系全体の節点座標値, Xˆ i：想定形状における系全体の節点
座標値, σ L：膜応力の下限値, σU ：膜応力の上限値 である。
なお、最適化手法には逐次二次計画(Sequential Quadratic.
Programming：SQP)法を採用する。
5.2 裁断図解析におけるリメッシュ
サスペンション膜構造や空気膜構造の裁断図の解析例は数多くある(文献[4]の文献参照のこと)。しかし、動的解析に移行する際、膜材の要素分割に関する検討は十分なされていない。従来、初期裁断線上の節点を設計変数とすると、最適結果の裁断線が滑らかな曲線にならないことがある。あるいは要素数を増やすと設計変数が増大し解が収束しないことが起こる。そこで、パラメトリック曲線を導入し、制御点を設計変数とすることで設計変数の数を減らす操作が行われている。しかし、パラメトリック曲線を用いた裁断線は、節点数を増やすと必ずしも解の収束状況がよくなるとは限らない。
これは裁断線が逆にパラメトリック曲線に支配されてしまう
結果と考えられる。ここでは、要素分割数と裁断図解析の関係を調べ、リメッシュ手順を示す。解析モデルは HP(Hyperbolic paraboloid)曲面を構成する膜構造を想定し、ライズの違いによる要素分割数と解の収束状況及び固有周期の違いを把握する。
5.3 リメッシュ手順
裁断図解析において要素分割を増加させる場合、次のリメ
ッシュ手順に従う。
step-1 裁断線にパラメトリック曲線を採用し、荒い要素分割
で裁断図解析を実施する。
step-2 得られた裁断線に対し、裁断線を保持した状態で要素を細かくリメッシュする。
step-3 リメッシュ後、再び裁断図解析を実施する。
step-4 固有振動数を調べ、一定値になるよう十分な要素分割

になるまで step-2 に進む。
5.4 解析モデルと数値結果
解析モデルは、周囲境界を全て固定とした HP 曲面を構成するサスペンション膜構造である。図５に解析モデルの鳥瞰
図･膜帯接続情報･要素分割図(case-1：23 節点, 32 要素)を示す。
リメッシュ手順ではモデル図中 R = 5, 6, 7 m (ライズ比:
0.25, 0.3, 0.35)の３タイプとする。膜材料の諸量は表１に示す通りである。膜応力の上下限値は、主応力に対して次のように設定する。

: 裁断線上の未知座標値(境界節点): 裁断線上の未知座標値(縫合線上節点)(33)

20 000

要素分割モデルは、初期要素分割モデル(case-1)の他に、
case-2(43 節点, 64 要素), case-3(77 節点, 128 要素), case-4(149
節点, 256 要素)の計４つ準備した。case-2 以降の裁断図解析における初期値はリメッシュ手順で説明したように 1 段階前の解を用いる。なお、連続して静的･動的解析を行うため、裁断図解析では全体構造で計算した。
各タイプのモデルを用い、種々の条件で解析を行った。パ
ラメトリック曲線は、スプライン曲線を採用する[8]。表２に解析結果の状況をまとめた(○：解析可, ×：解析不可)。表２にはリメッシュ手順の他、case-4 の要素分割で裁断線を直線とした初期値にスプライン曲線利用による解析と利用しない裁断線上の全節点座標値を設計変数とする解析の結果も示している。case-4 の要素分割状態で十分であるとの判断は、
図６に示す固有周期と要素分割関係を調べた結果、一定値に収束した理由による。リメッシュ手順による解析結果は、R =6 のとき、各段階の裁断図を図７に、最終結果における釣合
形状と主応力図は図８,９に示す。主応力図に対し、主応力の平均値･最大値･最小値を表３にまとめた。最小主応力の最小値は境界部で発生しているが、膜面全体では滑らかな力の流れを示している。なお、参考に R = 6 のとき、裁断線上の節点の座標値を設計変数とした裁断図を図 10 に示す。
裁断線がジグザグになっていることが判る。ただし、R = 5 のときは、比較的滑らかな裁断線が得られた。なお、case-1 の裁断線をリメッシュ手順によらず、単純にそのまま case-4 のごとく再分割した場合、釣合形状の応力状態は大きくばらつき、制約条件を満足しない結果が得られる。

12 000b) 膜帯接続情報図

20 000T
T
σ L = ⎣⎢ 2 2 0 ⎦⎥ , σU = ⎣⎢5 5 0 ⎦⎥a) 鳥瞰図3000 3000 3000 3000c) 初期裁断図及び要素分割モデル(case-1)

図５ 解析モデル 鳥瞰図･要素分割図 unit: mm
表１ 膜材定数表
膜厚
縦弾性係数
ポアソン比
せん断剛性
単位質量

t = 0.8 mm
Ex･t = 645.3 Ey･t = 213.8 ( N/mm)
νxy = 0.9055
νyx = 0.3
Gxy･t = 55.9 (N/mm)
1.215×10-6 (kg/mm2)

表２ リメッシュ手順による解析結果
従来の方法

ライズ

リメッシュ
手順

スプライン曲線

境界設計変数

H=5
H=6
H=7

○
○
○

×
×
×

○
○
×

○：解析可, ×：解析不可

5.6
5.5

R=7
R=6
R=5

5.4
5.3

Period (sec)６．動的解析
6.1 解析モデル
文献[3,4]では、既に座標仮定有限要素法による張力材のみで構成された張力膜構造の時刻歴応答解析の数値例が示されている。ここでは、図５のサスペンション膜構造のリメッ
シュ手順による解析結果(case-4; R = 6 m)(図 7)を利用し、懸垂物がある骨組材と張力材で構成する複合構造の動的解析を示す。解析モデルでは、図 11 に示すように膜境界部及び脚部に骨組を配置する。懸垂物は、骨組材から 6 本のワイヤー(ケーブル)で吊るすモデルを想定する( r = 4.5 m)。骨組及びケーブルの材料諸量は表４, ５に示す。なお、懸垂物は設備荷重

5.2
5.1
5.0
4.9
4.8
4.7
4.6

1

2

3

4

case No.
図６ リメッシュ手順による各段階の固有周期

case-2：43 節点, 64 要素

case-1：23 節点, 32 要素

case-4：149 節点, 256 要素

case-3：77 節点, 128 要素

図７ リメッシュ手順による各段階の最適裁断図(R = 6)

図９ case-4 の最適裁断図による
主応力図(R = 6)

図８ case-4 の最適裁断図による
釣合形状(R = 6)

図 10 case-4 の裁断線上の節点を設計変数とした最適裁断図(R = 6)

表３ 主応力値
主応力
目標応力
平均値
最大値
最小値

t σmaxt σmin2 ≤ t σmax , t σmin ≤ 5

4.59
5.0
3.3

2.77
4.57
2.0
unit : N/mm

図 11 時刻歴応答解析の構造モデル

表４ 骨組材定数表
断面積
断面二次モーメント
弾性係数
単位質量

unit:mm

2

80 883(mm )545 180 000 (mm4)210 000(N/mm2)-67.850 × 10 (kg/mm3)Y

表５ ケーブル(ワイヤー)材定数表
断面積
弾性係数
単位質量

210
(mm2)
210 000(N/mm2)-4
5.688 × 10 (kg/mm)を考慮して 1.5 ton の質量を設定した。
6.2 解析手順と設定荷重
膜シート部は式(10),(11)を用い、平面上にある裁断図を利用

O

X

図 12 風圧係数分布図

風向

12

22.5

膜応力 (N/mm )変動風速(m/sec )

25

20
17.5

t σmaxσmax
t σmin
σmin

9

6

3

15
0

5

10

15

20

25

30

0

5

10

15
time (sec )time(sec )2530

図 14 風荷重作用時の中央部膜応力(主応力)時刻歴応答
21.0

100
50
0
-50
-100
-150
-200
-250
-300

20.8

軸力(kN )

変位(mm )

図 13 風速データ(付図 1,点 10 位置)

20

Ux
Uy
Uz

0

5

10

15

time(sec )

20

25

20.6
20.4
20.2
20.0

30

図 15 風荷重作用時の中央部節点変位時刻歴応答する。骨組部は初期形態の位置情報から、式(28)を用いて座標値を未知量に変換する。懸垂物を支えるケーブル部は、初期形態の位置情報に基づき、式(8),(9)より構成する。懸垂物の扱
いは、自重に相当する静的荷重を鉛直方向に作用させ、ケーブル張力導入による幾何剛性保持の状態で剛性を評価する[9]。
１本のケーブルにおいては16 分割するモデルを採用し、非抗圧材に対する処理は行わない。以上の手続きにより、膜シート部、骨組部、ケーブル部に対し、式(29),(30)を用いて、構造全体の平衡方程式及び接線剛性行列を構成する。荷重は風を想定する。風荷重は、パワースペクトル･コヒーレンス･
フェイズを設定したとき、それらを多地点の風速変動の時系
列外力で与えられる。概要は付録に記す[10,11]。解析に用いた構造モデルの風圧係数分布は図12 の通りである[12]。
なお、減衰項は、Rayleigh 減衰を用い、減衰定数は 2%とした。時間
刻みは 1/500 秒としている。
6.3 解析結果
解析結果は、一例として風荷重作用時における膜部、ケーブル部の代表的な力学的挙動を示す。
図 13 は使用した変動風
速データ(付図１:点 10 のデータ)であり、図 14-16 は、風荷重作用時の膜中央部の応力、同じく節点変位及びケーブル軸
力の時刻歴応答結果である。
７．まとめ
張力材は初期張力を導入することにより構造剛性を保持する。このような張力材要素と一般構造要素を複合させた解析には工夫が必要である。本論では、座標仮定有限要素法を用

0

5

10

15

20

25

30

time(sec )図 16 風荷重作用時のケーブル軸力時刻歴応答
いた複合構造の離散化定式化を示し、膜材の裁断図解析と連続して動的解析の数値例を示した。張力材要素は、平面上に配置した自然状態を、一般構造要素は初期構造形態を基準にそれぞれを統合させた計算手順である。なお、要素分割が荒い場合と細かい場合では、時刻歴応答挙動が大きく異なったため、膜材の裁断図解析において要素分割の検討を行い、固有値解析の導入により再分割の手順を示した。時刻歴応答解
析は、変動風速データによる計算を行った。解析対象をケーブルで吊された懸垂物を有する複合構造とし、懸垂物を含む
力学的挙動が捉えられることを確認した。ただし、数値解析
上の現状の問題点は以下の通りである。
① 設計変数の数を削減するため、裁断線をスプライン曲線に置き換えた。また、要素形状を三角形(要素内ひずみ一定)とし、裁断線以外の節点は初期値に固定している。これらの影響で境界部付近の膜応力に乱れが見られる。
② 張力材に圧縮力が作用するとき、非抗圧材の処理が繰り返し行われ、計算時間の負荷が増大する。
これらの問題は、要素形状を矩形にする、非抗圧材の処理の回数を制限し、次ステップに処理を任せるなどにより解決することができると考えている。
以上、張力構造一般(複合構造)の座標仮定を基準とした有限要素法による基本的な解析ツールを開発することができた。
一般構造要素としては、梁要素だけでなく、シェル要素や連続体要素にも適用可能である。なお、ここで示した解析例では膜面中央部の変動が大きく、設定膜張力の変更やケーブル

補強の検討が必要である。このように、今後、張力構造の解析対象モデル毎あるいは張力構造における懸垂物の安全性や構造物への影響など、詳細な力学的挙動の検討を行っていきたい。

付録
変動風速シミュレーション
1. 風荷重 Wc (N)Wc = q H Cc Ac(a)

ここで、Cc :風力係数, Ac:受圧面積(m2)であり、速度圧 qH は次式

謝辞
風速変動シミュレーションでは、日本大学教授 岩谷祥美博士より助言を得ました。ここに深謝いたします。

で与えられる。
qH =1
ρU H 2
2( N/m2 )参考文献

ここで、ρ:空気密度 1.22 ( kg/m3 ), U H :設計風速 である。

[1] T. Honma, N. Ataka：Geometrically Nonlinear Structural Analysis by
FEM Using the Coordinate Value on a Deformed Body, Information,
7(5), 569-584, 2004
[2] 本間俊雄, 安宅信行：座標値を未知量とした有限要素法による張力構造の解析と評価, 膜構造論文集, 18, 15-21, 2005
[3] 本間俊雄, 合田雄策, 安宅信行：座標値を未知量とした有限要素技術による張力構造解析の一方法, 日本建築学会構造系論文集,
602, 161-169, 2006
[4] 本間俊雄, 森哲也, 坂中玲子:膜構造の裁断図解析と静的･動的な応力･変形解析及び発想･設計支援システムについて, 膜構造論文
集, 21, 1-13, 2007
[5] A. Neuenhofer, F.C. Filippou, Evaluation of nonlinear frame finite
element models, Journal of structural engineering, 123(7), 958-966,
1997[6] 山本憲司, 大森博司, 皆川洋一：座屈荷重を目的関数とする空間構造の形状最適化に関する研究,日本建築学会構造系論文集, 564,
95-102, 2003
[7] 藤井文夫,瀧諭,萩原伸幸,本間俊雄,三井和男：非線形構造モデルの動的応答と安定性,計算工学シリーズ 2, コロナ社, 2003
[8] 市田浩三, 吉本富士市：スプライン関数とその応用, 教育出版,
1979[9] 大木洋治, 元結正次郎：ケーブルによって吊り下げられた質点の繰り返し外力下での挙動, シェル･空間構造の耐震･耐風･耐雪･耐火を考えた性能設計と解析論文集(平成 11 年度京都大学防災研究
所共同研究集会), 1-7, 1999
[10] 岩谷祥美：任意のパワースペクトルとクロススペクトルをもつ多次元の風速変動のシミュレーション, 日本風工学会誌, 11, 5-18,
1982[11] 日本建築学会編：建築物荷重指針･同解説, 12-65, 2004[12] N. Esquillan, Y. Saillard：Hanging Roofs, proceedings of the IASS colloquium on
hanging roofs, continuous metallic shell roofs and superficial lattice roofs, Paris

2. 設計風速U H = U 0 K D E H k rw = 18.9 ( m/sec )

(b)

(c)

U 0 :基本風速 20 m/sec (地表面粗度区分Ⅱ, 地上 10 m における再現期間 100 年に対する値), K D : 風向係数 1.0, E H :風速の鉛直分
布係数 0.924 (鉛直分布係数の基準高さ H における値), 再現期間換
算係数 k rw :1.0
3. 変動風速シミュレーション
設計風速 18.9 m/sec, 時間刻み 0.1 sec, 点の数 10, 基本高さ
(点１の高さ) 0.7 m, 点間隔 0.7 m の条件により、統計処理に基づき、各点の変動風速の計算[11]後、各節点に相当する風速を線形補間により決定する(付図１)。得られた風速は式(b)に代入し、速度圧を求め､式(a)より風荷重を算出する。

9-11, 1962

Z
X

O
付図１ 解析モデルと変動風速モデルとの対応

DYNAMIC ANALYSIS ON TENSION STRUCTURE
USING FINITE ELEMENT TECHNIQUE WITH COORDINATES ASSUMPTION
Toshio Honma *1) , Yusuke Saeki *2) , Reiko Sakanaka *2)
SYNOPSIS
We have shown examples of a cutting pattern analysis and a dynamic analysis using the finite element technique of which an unknown quantity is thecoordinate value on a stability form. This procedure can examine the mechanical behavior of the tension structure in the 3 dimensional space based on the stateof the tension material (membrane and cable) that exists on the plane. In this paper, the example of the dynamic analysis on the tension structure (hybrid
structure) to take the beam element to this finite element technique is shown. The beam element cannot do the discretization formulation by the finite element
technique with coordinates assumption usually. However, it is shown that it is possible by a simple contrivance. The decision of the cutting line obtains themodel of a detailed element breakdown by starting from a comparatively rough element breakdown, and repeating the cutting pattern analysis. The dynamic
analysis model is a hybrid structure with a mass suspended by wires, using the obtained element breakdown. The dynamic external forces use the fluctuating
wind velocity data.
*1)*2)

Professor, Dept. of Architecture & Architectural Engineering, Faculty of Engineering, Kagoshima Univ., Dr. Eng.
Graduate student, Graduate School of Science and Engineering, Kagoshima Univ.