跳轉至

外力虛功與全域方程式的組裝

在內力虛功的離散化中，已將弱形式左側整理為元素內力向量 \(\boldsymbol{q}^e\)（UL 法）或 \(\boldsymbol{Q}^e\)（TL 法）。本章由外力虛功導入元素節點外力向量 \(\boldsymbol{F}^e\)，再經由依全域節點編號重新排列並彙總元素節點物理量的組裝操作，得到 FrontISTR 非線性結構分析所需求解的節點位移非線性方程組。

外力虛功的元素分解

虛功原理的右側可依元素分解為由體積力（物體力）以及力學邊界上的指定表面力所構成的外力虛功。為了以矩陣形式表示形狀函數與有限元素近似中導入的位移插值，令節點 \(\alpha\) 的形狀函數 \(N_\alpha^e\) 排列在對角線上形成 \(d \times d\) 區塊 \(\boldsymbol{N}_\alpha\)，並將其橫向排列成 \(\boldsymbol{N} = [\boldsymbol{N}_1, \ldots, \boldsymbol{N}_{n_e}]\)，使 \(\delta\boldsymbol{u} = \boldsymbol{N}\, \delta\boldsymbol{u}^e\)。將其代入以參考配置表示的外力虛功，可得

\[ \delta W^{\mathrm{ext}} = \sum_e \delta\boldsymbol{u}^{eT} \boldsymbol{F}^e, \qquad \boldsymbol{F}^e_\alpha = \int_{\Omega^e_0} \boldsymbol{N}_\alpha^T \rho_0 \boldsymbol{g}\, dV + \int_{\Gamma^e_{0t}} \boldsymbol{N}_\alpha^T \bar{\boldsymbol{t}}_0\, d\Gamma_0 \]

其中，元素節點外力向量排列為 \(\boldsymbol{F}^e = (\boldsymbol{F}^{eT}_1, \ldots, \boldsymbol{F}^{eT}_{n_e})^T\)。如此可將外力虛功整理成與內力側相同的「元素節點向量 × 測試函數」形式（即使以現行配置表示，只需作 \(dV \to dv\)、\(\rho_0 \to \rho\)、\(\bar{\boldsymbol{t}}_0 \to \bar{\boldsymbol{t}}\) 的替換，即得到相同形式）。

元素節點物理量的組裝

將各元素得到的節點物理量 \(\boldsymbol{Q}^e_\alpha, \boldsymbol{F}^e_\alpha\) 彙總到依全域節點編號排列的全域向量中。令元素 \(\Omega^e\) 的元素內節點編號 \(\alpha\) 所對應的全域節點編號為

\[ \mathrm{gdx}(e, \alpha) = i_g \]

如此，元素節點物理量便與全域節點物理量的對應成分一致（例如 \(\boldsymbol{u}^e_\alpha = \boldsymbol{u}_{i_g}\)）。一般而言，節點 \(i_g\) 會由多個元素共享，因此定義由全域節點編號為 \(i_g\) 的 \((e, \alpha)\) 組合所構成的集合

\[ \mathcal{E}(i_g) = \{ (e, \alpha) \mid \mathrm{gdx}(e, \alpha) = i_g \} \]

利用此集合將總和改寫為 \(\sum_e \sum_\alpha = \sum_{i_g} \sum_{(e,\alpha) \in \mathcal{E}(i_g)}\)，即可得到跨越全部 \(n_g\) 個節點的節點內力與全域內力向量

\[ \boldsymbol{Q}_{i_g} = \sum_{(e,\alpha) \in \mathcal{E}(i_g)} \boldsymbol{Q}^e_\alpha, \qquad \boldsymbol{Q} = (\boldsymbol{Q}^T_1, \ldots, \boldsymbol{Q}^T_{n_g})^T \]

其中，\(\boldsymbol{Q}_{i_g}\) 相當於作用於節點 \(i_g\) 的元素節點內力之合力；若無外力作用且處於平衡狀態，則為 \(\boldsymbol{0}\)。UL 法也可用相同程序得到 \(\boldsymbol{q}_{i_g}, \boldsymbol{q}\)，且其數值滿足 \(\boldsymbol{q} = \boldsymbol{Q}\)，因此以下除非需要區分，統一使用 \(\boldsymbol{Q}\) 表示。全域外力向量 \(\boldsymbol{F}\) 也由相同的彙總方式得到。

實作上不會明確建立集合 \(\mathcal{E}(i_g)\)，而是在元素迴圈中加總至對應成分。

將全域內力向量 Q 初始化為 0：Q_{i_g} = 0  (i_g = 1, ..., n_g)
for e = 1 to (元素數)
    for α = 1 to n_e
        i_g = gdx(e, α)
        Q_{i_g} += Q^e_α
    end for
end for

全域外力向量 \(\boldsymbol{F}\) 也以相同程序建構。這種將元素節點物理量加總並儲存至以全域節點編號編號的向量與矩陣中的操作，稱為組裝（assemble）。對於關聯兩個節點編號的二階張量（例如剛性矩陣），可利用集合 \(\mathcal{E}^2(i_g, i_h) = \{ (e, \alpha, \beta) \mid \mathrm{gdx}(e, \alpha) = i_g\ \mathrm{and}\ \mathrm{gdx}(e, \beta) = i_h \}\) 進行相同類型的組裝操作（具體建構請參閱切線剛性矩陣）。

待求解的非線性方程式

將內力與外力的組裝結果代入虛功原理，並利用其對滿足幾何邊界條件的任意測試函數 \(\delta\boldsymbol{u}^n\) 皆成立，可得

\[ \boldsymbol{Q}(\boldsymbol{u}^n) - \boldsymbol{F}(\boldsymbol{u}^n) = \boldsymbol{0} \]

在增量分析（增量分析架構）的脈絡中，恢復時間下標 \(_{n+1}\)，並省略表示全域節點向量的上標 \(^n\)，則待求解的方程式為

\[ \boldsymbol{Q}(\boldsymbol{u}_{n+1}) - \boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{0} \]

因此，求取時刻 \(t_{n+1}\) 之節點位移 \(\boldsymbol{u}_{n+1}\) 的離散化邊界值問題，可歸結為將此關於位移的非線性方程式與幾何邊界條件一併求解。方程式的線性化與切線剛性矩陣的建構請參閱切線剛性矩陣，迭代求解法請參閱 Newton-Raphson 法。

相關項目