第三章 數值計算方法與程序
3.6 數值程序
本文使用的增量迭代法之數值程序可以分為三個主要部分 1. 輸入與計算開始分析所需的資料
(a) 輸入結構與負荷資料。
(b) 選擇一個參考自由度，並設定期望此自由度達到的最大位移。
(c) 給定最大增量數、每個增量期望的迭代數與最大的迭代數、收 斂時的容許誤差。
(d) 形成剛度矩陣並求得(3.10)式中的R 。0
(e) 利用(3.1)-(3.3)式、(3.9)式與(3.10)式計算第一次增量的增量弧 長、增量負荷參數與初始增量位移。
2. 使用迭代法求增量的收斂解
(a) 利 用 已 知 的 增 量 位 移 求 得 目 前 元 素 的 變 形 ， 並 以(2.39)、
(2.45)式計算元素節點內力。
(b) 計算(2.54)式的不平衡力Ψ 。
(c) 檢查(2.55)式的收斂準則，若滿足則進行(e)；否則檢查迭代數，
如果小於給定之最大迭代數，則進行(d)；否則減少增量弧長並以
(3.1)-(3.3)式與(3.9)式計算新的增量位移與增量負荷，回到步驟(a)重 新此增量。
(d) 以(3.13)式與3.2節的方法計算增量位移修正量與增量負荷參數
修正量，然後回到步驟(a)。
(e)判斷系統切線剛度的行列式值是否變號，並使用二分法來偵測分 歧點，若找到第一個分歧點，則進行(f)，否則進入第三部分。
(f)使用逆冪法計算挫屈模態，進入第3部分，且此後第二部分改用 標準牛頓法進行迭代。
3. 計算下一次增量所需的資料
(a) 檢查參考自由度的位移或已進行的增量次數是否已達給定值。
若已達到給定值則停止分析工作；否則進行下一步工作。
(b) 計算(3.12)式中的切線剛度矩陣與FP。
(c) 若第二部分有進行(f)步驟，則調整所算出的挫屈模態向量的長 度，並以此挫屈模態作為下一增量的增量位移，擾動進入次要路徑。
若第二部分沒有進行(f)步驟，則以(3.1)-(3.3)式與(3.8)-(3.9)式計算下 一次增量的增量位移、增量弧長與增量負荷參數。
(d) 回到第二部份執行迭代工作。
第四章 數值分析與結果
在本章中，將先以文獻上的數值例題驗證本文第二章之板元素及第 三章所提的幾何非線性分析之數值程序與方法的正確性及可行性，再 以文獻上的例題探討膜結構受到剪位移負荷作用下的非線性行為。
例題一：圓柱殼受到位移負荷作用
圖4.1 為圓柱殼結構示意及其受之位移負荷圖，如圖4.1所示，AD， BC 兩邊為鉸接，AB，CD 兩邊為自由邊，本例題在圖 4.1 的中點 M 施加一個向下的位移負荷λM ，其中λM =λ(mm)，λ 是負荷參數 (Loading parameter)，由於結構為雙對稱結構，本例題分別取1/4個結 構，與 1/2個結構，還有完整結構進行分析。分析結構平衡路徑的主 要路徑時，本例題使用10×10元素網格，容許誤差值取1×10−4，取 1/4個結構分析時，使用 20 個增量，平均迭代次數約為4，取 1/2個 結構分析時，使用 20 個增量，平均迭代次數約為 4，取完整結構分 析時，使用 23 個增量，平均迭代次數約為 4。分析結構平衡路徑的 次要路徑時，本例題使用10×10元素網格，容許誤差值取1×10−4，取 1/2個結構分析時，使用 27個增量，平均迭代次數約為4，取完整結 構分析時，使用33個增量，平均迭代次數約為4。圖4.2為本例題所 得到 M 點的位移-負荷曲線圖的結果。由圖 4.2 可以看出本例題在取 1/2 個結構或整個結構作分析時，都可以偵測到平衡路徑上的分歧 點，並可以挫屈模態擾動進入結構的次要平衡路徑。但在取1/4個結 構分析時，無法找到平衡路徑上的分歧點，也無法找到結構的次要平 衡路徑，這應是因為若僅取1/4結構進行分析，在設定邊界條件的時 候，將圖4.1 上 M點的X 與 Y方向位移鎖住了，所以不會偵測出結 構平衡路徑上的分歧點。圖 4.3 為本例題與文獻的比較，其中RM 為
在位移負荷方向的反力，RM及 λ相當於文獻之力負荷及力負荷方向 的位移，由圖4.3 可以看出本例題的結果與文獻[18]的結果相當重合。
例題二：正方形膜結構邊界受剪位移負荷
圖4.4 所示為本例題分析之正方形膜結構與受負荷圖，本例題是參 考文獻[11]提出的例題，正方形膜結構邊長L=229 (mm)，厚度
0762 . 0
h= (mm)，楊氏係數E=3790 (MPa)，蒲松比ν=0.38。如圖 4.4，本例題AB，CD邊為自由邊，BC邊為固接(Fixed End)，AD邊 為固接，但可在X軸方向上移動，本例題在 AD邊上施加均勻的水平 位移λ (mm) (圖4.4)。λ稱為負荷參數(loading parameter)，在施加剪位 移λ的過程中，AD 與BC兩邊之間的距離是維持固定不變的。
本例題的挫屈負荷為切線剛度行列式值的正負號恰欲由正變負時 的位移負荷，40×40元素網格挫屈負荷值為6.56×10−4 (mm)，60×60 元素網格挫屈負荷值為6.55×10−4 (mm)，80×80元素網格挫屈負荷 值為6.549×10−4 (mm)，可以看出使用80×80元素網格挫屈負荷的值 收斂了，本例題使用80×80元素網格對膜結構進行分析，並與使用
40
40× 以及60×60元素網格的結果以及文獻[11,12]作比較。使用 40
40× 元素網格例題在結構挫屈後使用標準牛頓法進行迭代，平衡迭 代容許誤差為5×10−5，使用 110 個增量，每個增量平均迭代次數為 約為 16 次，使用60×60元素網格例題結構挫屈後使用標準牛頓法進 行迭代，平衡迭代容許誤差為5×10−5，使用 202 個增量，每個增量 平均迭代次數為6，使用80×80元素網格例題結構挫屈後使用標準牛 頓法進行迭代，平衡迭代容許誤差為1×10−4，使用 203 個增量，每 個增量平均迭代次數為5。圖 4.5-4-7為膜結構的第一個挫屈模態，可 以看出使用不同元素網格計算出來的挫屈模態型態非常相似。圖 4.8
為挫屈模態的剖面圖。圖4.9-4.11為膜結構正中央點I 點之位移wI與 AB 邊中點 E 點之位移w -E 負荷參數 λ 曲線圖，I 點與 E 點位置如圖 4.4所示，由圖 4.9-4.11可以看出wI隨著負荷參數λ增加一直在增加，
且使用40×40元素網格與60×60元素網格以及80×80元素網格wI的 方向相反。wE在初期會有不穩定的跳動現象，之後則趨於穩定，圖 4.12是使用60×60元素網格，wE在使用不同的增量負荷參數∆ (mm)λ 進行迭代的結果，可以看出wI曲線是趨於穩定的，但wE會有不同的 結果，這可能是膜結構自由邊邊界鬆弛了，所以使的結構受到差異極 小的負荷條件，但結構的變形差異卻很大，圖 4.13 是使用60×60元 素網格以及80×80元素網格的比較圖，可以觀察到wI的曲線是接近 的，但wE的曲線仍有差異。圖 4.14(a)為位移邊界 AD 之節點合力
∑
RXi –位移負荷參數λ 曲線圖，圖4.14(b)為位移邊界 AD之節點合 力∑
RYi–位移負荷參數λ 曲線圖，圖4.14中之反力∑
RXi，∑
RYi各為圖4.4 中的位移邊界AD之節點 X與Y 方向反力的和，由圖4.14 可看出不同網格下，位移邊界AD 之節點反力的合力
∑
RXi，∑
RYi– 位移負荷參數λ 的曲線相當接近。圖 4.15-4.17 為位移邊界 AD 在 X 與Y軸方向上的節點反力。由圖 4.15-4.17 可以看出越靠近A點的地 方節點反力越大，圖4.18-4.22是膜結構變形圖，由圖4.18-4.22可以 看出皺折的走向是呈AC的方向，膜結構最大拉應力的分佈的應是沿 著對角線AC的走向，隨著位移邊界 AD上的點越靠近 D點，拉應力 逐漸變小，故位移邊界上之節點力的分佈是合理的。圖4.18-4.22為等高線圖及透視圖，透視圖是使用視角53o繪製。膜
結構側方向最大正負位移，以及最大正負位移所在位置的初始座標如 表4.1 所示。由圖4.18-4.22可以看出本例題使用60×60與80×80元素 網格，膜結構的皺折數目較使用40×40元素網格多。圖 4.23-4.36 為
結構重要軸線在 Z 方向上的變形，由圖 4.23-4.36 可以看出隨著負荷
小位移大約相差4.33(mm)，以本例題與文獻[11,12]比較，可以看出本 例題以40×40元素網格分析的結果與文獻[11]數值例題的結果較接 近，都有三條突起的皺折，膜結構側方向最大正負位移的差也接近，
但是本例題以40×40元素網格分析的膜結構側方向最大正位移比起 文獻[11]的結果大了許多。本例題使用60×60與80×80元素網格分析 的結果與文獻[12]實驗結果較接近。主要有四條主要突起的皺折，膜 結構側方向最大正負位移的差也接近，但是本例題以80×80元素網格 分析的側方向最大正位移比起文獻[12]的結果大，最大負位移比起文 獻[12]的結果小，比較本例題使用80×80元素網格分析的結果與文獻 [12]實驗結果，可以發現兩者的皺折型態有相同的特徵，但細部仍有 差異。造成差異的原因有可能是本例題使用80×80元素網格還不夠 密，也有可能是文獻[12]在做此實驗時，有先將膜結構稍微彎曲一下，
才開始施加剪位移，這部分在本例題的數值計算中沒有被考慮到，也 有可能文獻[12]實驗真正的邊界條件與本例題所設定的邊界條件有差 異，也有可能是膜結構某些地方很容易鬆弛，差異極小的負荷條件也 會造成不同的變形結果。
第五章 結論與展望
本文以共旋轉有限元素推導法，殼元素來分析膜結構的皺折，由本 文分析之數值例題的結果，可得以下的結論：
(1)膜結構邊界受剪力側向挫屈後有很多不同平衡路徑，且膜結構受 力後的變形很敏感，差異很小的負荷條件，有可能會造成差異很大的 變形結果。
(2)由本文例題結果發現膜結構受剪位移時，若元素網格不夠密，則 膜結構皺折的數目太少。若欲得到精確的結果，應使用較密元素網格。
(3)本文對於平面膜結構產生側方向位移的機制為在分歧點加上挫 屈模態擾動，使膜結構有側方向的位移，本研究發現挫屈模態與膜結 構挫屈後皺折的的形狀差異很大。
(4)由本研究分析的結果與文獻上的實驗結果比較，可以發現兩者的 皺折型態有相同的特徵，但細部仍有差異，其原因可能是本研究分析 時，元素的數目不夠，文獻上的結構有初始不完美，且實驗操作上真 實的邊界條件可能與分析的邊界條件不同。
本文例題發現使用太密的元素網格，會出現不易收斂的情形，這可 能是因為本研究以平板三角殼元素來模擬薄膜結構的行為，可能因平 板殼元素沒有考慮撓曲變形與膜變形間的耦合及採用近似的幾何剛 度，造成收斂的困難。也可能是個數值方法上的問題，未來可考慮採 用其他元素或數值方法，若能解決這些困難，使用較密的元素網格，
有可能可以得到更精確的結果。
參考文獻
[1] J. W. Leonard, “Tension Structure.” McGraw-Hill Book Company;
1988.
[2] H. Ding, B. Yang, “The modeling and numerical analysis of wrinkled membranes.” Int J Numer Meth Engng. 2003; 58: 1785-1801.
[3] E. S. Lee, S. K. Youn, “Finite element analysis of wrinkling membrane structures with large deformations.” Finite Elements Anal Design. 2006; 42: 780-791.
[4] R. Rossi, M. Lazzari, R. Vitaliani, And E. Onate, “Simulation of light-weight membrane structures by wrinkling model.” Int J Numer Meth Engng. 2005; 62: 2127-2153.
[5] J. Månsson, J. Söderqvist, “Finite element analysis of thin
[5] J. Månsson, J. Söderqvist, “Finite element analysis of thin