平面應力問題
平面應力問題的平面應力應變關係
\[\begin{Bmatrix}\sigma_{xx}\\\\\sigma_{yy}\\\\\tau_{xy}\end{Bmatrix}=\frac{E}{1-\gamma^2}\begin{bmatrix}1&\gamma&0\\\\\gamma&1&0\\\\0&0&\frac{(1-\gamma)}{2}\end{bmatrix}\begin{Bmatrix}\varepsilon_{xx}\\\\\varepsilon_{yy}\\\\\gamma_{xz}\end{Bmatrix}
\]
平面應變問題的應力-應變關係
\[\begin{Bmatrix}\sigma_{xx}\\\\\sigma_{yy}\\\\\tau_{xy}\end{Bmatrix}=\frac{E}{(1+v)(1-2v)}\begin{bmatrix}1-v&-v&0\\\\-v&1-v&0\\\\0&0&\frac{(1-2v)}{2}\end{bmatrix}\begin{Bmatrix}\epsilon_{xx}\\\\\epsilon_{yy}\\\\\gamma_{xy}\end{Bmatrix}
\]
小變形下的應力-位移關係
\[\begin{aligned}&\varepsilon_{xx}=\frac{\partial u}{\partial x}\\&\varepsilon_{yy}=\frac{\partial\nu}{\partial y}\\&\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial\nu}{\partial x}\end{aligned}
\]
寫成矩陣形式
\[\begin{Bmatrix}\varepsilon_{xx}\\\\\varepsilon_{yy}\\\\\gamma_{xy}\end{Bmatrix}=\begin{bmatrix}\frac{\partial}{\partial x}&0\\\\0&\frac{\partial}{\partial y}\\\\\frac{\partial}{\partial y}&\frac{\partial}{\partial x}\end{bmatrix}\begin{Bmatrix}u\\\\\nu\end{Bmatrix}
\]
簡寫
\[\{\epsilon\}=[L]U
\]
對於n節點的單元，未知位移的節點近似插值函式為
\[u=N_{1}u_{1}+N_{2}u_{2}+\cdots+N_{n}u_{n}\\\nu=N_{1}\nu_{1}+N_{2}\nu_{2}+\cdots+N_{n}\nu_{n}
\]
矩陣形式
\[\begin{Bmatrix}u\\v\end{Bmatrix}=\begin{bmatrix}N_1&0&|&N_2&0&|&\ldots&|&N_n&0\\0&N_1&|&0&N_2&|&\ldots&|&0&N_n\end{bmatrix}\begin{Bmatrix}u_1\\\nu_1\\\nu_2\\u_2\\\nu_2\\\vdots\\u_n\\\nu_n\end{Bmatrix}
\]
簡寫
\[\{U\}=[N]\{a\}
\]
用上式中U代替 \(\{\epsilon\}=[L]U\)中的U
\[\{\epsilon\}=[L]U\\
\qquad=[L][N]\{a\}\\
\qquad=[B]\{a\}\\
\]
其中
\[[B]=\begin{bmatrix}\frac{\partial N_1}{\partial x}&0&|&\frac{\partial N_2}{\partial x}&0&|&\ldots&|&\frac{\partial N_n}{\partial x}&0\\\\0&\frac{\partial N_1}{\partial y}&|&0&\frac{\partial N_2}{\partial y}&|&\ldots&|&0&\frac{\partial N_n}{\partial y}\\\\\frac{\partial N_1}{\partial y}&\frac{\partial N_1}{\partial x}&|&\frac{\partial N_2}{\partial y}&\frac{\partial N_2}{\partial x}&|&\ldots&|&\frac{\partial N_n}{\partial y}&\frac{\partial N_n}{\partial x}\end{bmatrix}
\]
矩陣[B] 被稱為應變矩陣，可以透過插值函式\(N_i(x,y)\)求導得到
為了得到單元單元上的載荷與節點位移之間的關係，我們用虛功原理
\[\int\limits_{V_{e}}\delta\{\epsilon\}^{\mathrm{T}}\{\sigma\} \mathrm{d}V=\int\limits_{V_{e}}\delta\{U\}^{\mathrm{T}}\{b\} \mathrm{d}V+\int\limits_{\Gamma_{e}}\delta\{U\}^{\mathrm{T}}\{t\} \mathrm{d}\Gamma+\sum\limits_{i}\delta\{U\}_{(\{x\}=\{\overline{x}\})}^{\mathrm{T}}\{P\}_{i}
\]
式中：
\(\{\epsilon\}\)為應變向量；\(\{\sigma\}\)為應力向量；\(\{U\}\)為位移向量
\(\{b\}\)為體力向量；\(\{t\}\)為表面力向量
\(\{P\}_i\)為作用在 \(\{x\}=\{\overline{x}\}\) 處的集中力向量
d\(V\)為單元體積；d\(\Gamma\) 為表面力\(\{t\}\)所作用單元的單元邊界
應變的變分 δ{ε}和位移的變分 δ{U} 可以分別表示為
\[\delta\{\epsilon\}=\delta([B]\{a\})=[B]\delta\{a\}\\
\delta\{U\}=\delta([N]\{a\})=[N]\delta\{a\}
\]
應力應變關係
\[\{\sigma\}=[D]\{\epsilon\}=[D][B]\{a\}
\]
代入虛功原理方程
\[\int\limits_{V_{e}}\delta\{a\}^{\mathrm{T}}[B]^{\mathrm{T}}[D][B]\{a\} \mathrm{d}V=\int\limits_{V_{e}}\delta\{a\}^{\mathrm{T}}[N]^{\mathrm{T}}\{b\} \mathrm{d}V+\int\limits_{\Gamma_{e}}\delta\{a\}^{\mathrm{T}}[N]^{\mathrm{T}}\{t\} \mathrm{d}\Gamma+\sum\limits_{i}\delta\{a\}^{\mathrm{T}}[N_{((x)=(\bar{x}))}]^{\mathrm{T}}\{P\}_{i}
\\
原方程\\
\int\limits_{V_{e}}\delta\{\epsilon\}^{\mathrm{T}}\{\sigma\} \mathrm{d}V=\int\limits_{V_{e}}\delta\{U\}^{\mathrm{T}}\{b\} \mathrm{d}V+\int\limits_{\Gamma_{e}}\delta\{U\}^{\mathrm{T}}\{t\} \mathrm{d}\Gamma+\sum\limits_{i}\delta\{U\}_{(\{x\}=\{\overline{x}\})}^{\mathrm{T}}\{P\}_{i}
\]
注意，對於平面單元，單元體積 d\(V\)和單元邊界 dГ 可以分別寫成 d\(\nu=t\)d\(A\)和 d\(\Gamma=t\)d\(l\) ,式中\(t\)
為單元的厚度，d\(A\)為單元的無限小面積，d\(l\)為單元的無限小邊界長度。
由於\(\delta\{a\}\) 是節點位移的變分，因此關於座標獨立，可以把它從積分符號中拿出來並消除，方程就變為
\[\left[\int\limits_{\mathcal{A}_{c}}[B]^{\mathrm{T}}[D][B]t\mathrm{d}A\right]\{a\}=\int\limits_{\mathcal{A}_{c}}[N]^{\mathrm{T}}\{b\}t\mathrm{d}A+\int\limits_{\mathcal{L}_{c}}[N]^{\mathrm{T}}\{t\}t\mathrm{d}l+\sum\limits_{i}[N_{(\{x\}=\{\overline{x}\})}]^{\mathrm{T}}\{P\}_{i}
\]
寫成矩陣形式
\[[K_{e}]\{a\}=f_{e} \\
\]
\[[K_e]=\left[\int\limits_{A_e}[B]^\mathrm{T}[D][B]t\mathrm{d}A\right]
\]
上式為剛度矩陣，其中
\[\{f_{e}\}=\int_{A_{e}}[N]^{\mathrm{T}}\{b\}t\mathrm{d}A+\int_{L_{e}}[N]^{\mathrm{T}}\{t\}t\mathrm{d}l+\sum_{i}[N_{((x)=\{\overline{x}\})}]^{\mathrm{T}}\{P\}_{i}
\]
上式為單位力向量
空間離散
常應變三角形(CST)單元
單元的插值函式
\[\begin{aligned}
N_{1}(x,y)& =m_{11}+m_{12}x+m_{13}y \\
N_{2}(x,y)& =m_{21}+m_{22}x+m_{23}y \\
N_{3}(x,y)& =m_{31}+m_{32}x+m_{33}y
\end{aligned}
\]
其中
\[\begin{gathered}
m_{11} =\frac{x_{2}y_{3}-x_{3}y_{2}}{2A} m_{12} =\frac{y_{2}-y_{3}}{2A} m_{13} =\frac{x_{3}-x_{2}}{2A} \\
m_{21} =\frac{x_3y_1-x_1y_3}{2A} \text{m22} =\frac{y_3-y_1}{2A} \text{m23} =\frac{x_1-x_3}{2A} \\
m_{31} =\frac{x_{1}y_{2}-x_{2}y_{1}}{2A} m_{32} =\frac{y_{1}-y_{2}}{2A} m_{33} =\frac{x_2-x_1}{2A}
\end{gathered}
\]
\[A=\frac{1}{2}\text{det}\begin{bmatrix}1&x_1&y_1\\1&x_2&y_2\\1&x_3&y_3\end{bmatrix}
\]
位移場
三角形單元的位移場可以表示為
\[u=N_{1}u_{1}+N_{2}u_{2}+N_{3}u_{3}\\\nu=N_{1}\nu_{1}+N_{2}\nu_{2}+N_{3}\nu_{3}
\]
矩陣表示
\[\begin{Bmatrix}u\\\nu\end{Bmatrix}=\begin{bmatrix}N_1&0&|&N_2&0&|&N_3&0\\0&N_1&|&0&N_2&|&0&N_3\end{bmatrix}\begin{Bmatrix}u_1\\\nu_1\\u_2\\\nu_2\\u_3\\\nu_3\end{Bmatrix}
\]
簡寫為
\[\{U\}=[N]\{a\}
\]
應變矩陣
\[\{\epsilon\}=[B]\{a\}
\]
其中
\[\begin{gathered}[B]=\begin{bmatrix}\frac{\partial N_1}{\partial x}&0&|&\frac{\partial N_2}{\partial x}&0&|&\frac{\partial N_3}{\partial x}&0\\0&\frac{\partial N_1}{\partial y}&|&0&\frac{\partial N_2}{\partial y}&|&0&\frac{\partial N_3}{\partial y}\\\frac{\partial N_1}{\partial y}&\frac{\partial N_1}{\partial x}&|&\frac{\partial N_2}{\partial y}&\frac{\partial N_2}{\partial x}&|&\frac{\partial N_3}{\partial y}&\frac{\partial N_3}{\partial x}\end{bmatrix}\end{gathered}
\]
代入得
\[[B]=\begin{bmatrix}m_{12}&0&|&m_{22}&0&|&m_{32}&0\\0&m_{13}&|&0&m_{23}&|&0&m_{33}\\m_{13}&m_{12}&|&m_{23}&m_{22}&|&m_{33}&m_{32}\end{bmatrix}
\]
備註：矩陣[B]與笛卡兒座標系的\(x\)和 y 座標無關，它只是節點座標的函式，並且在整個單元內為常量，因此在整個單元內應變向量也是常量。這就是為什麼這種單元被稱為“常應變三角形單元”的原因，本書中將常應變三角形單元簡稱為 CST(Constant Strain Triangle)單元。
由於矩陣[B]和[D]都是常數矩陣,因此剛度矩陣可表示為
\[[K_e]=[B]^\mathrm{T}[D][B]tA_e
\]
其中 \(A_e\)為單元的面積
單位力向量
體力
考慮體力{b}是重力引起的，因此
\[\int\limits_{A_{\epsilon}}[N]^{\mathrm{T}}\{b\}t \mathrm{d}A=t\int\limits_{A_{\epsilon}}\begin{bmatrix}N_{1}&0\\0&N_{1}\\N_{2}&0\\0&N_{2}\\N_{3}&0\\0&N_{3}\end{bmatrix}\begin{Bmatrix}0\\0\\-\rho g\end{Bmatrix}\mathrm{d}A=t\begin{bmatrix}0\\-\int\limits_{A_{\epsilon}}N_{1}\rho g \mathrm{d}A\\0\\-\int\limits_{A_{\epsilon}}N_{2}\rho g \mathrm{d}A\\0\\-\int\limits_{A_{\epsilon}}N_{3}\rho g \mathrm{d}A\end{bmatrix}
\]
這裡採用式(8.18)和式(8.19)給出的三角形積分公式，計算在整個單元上的插值函式的積分。運用這些公式，可將積分結果整理如下：
\[\int_{A_{e}}N_{1}\rho g \mathrm{d}A=\rho g\int_{A_{e}}N_{1}^{1}N_{2}^{0}N_{3}^{0} \mathrm{d}A=\rho g\frac{1!0!0!}{(1+0+0+2)!}2A_{e}=\rho g\frac{A_{e}}{3}\\\int_{A_{e}}N_{2}\rho g \mathrm{d}A=\rho g\int_{A_{e}}N_{1}^{0}N_{2}^{1}N_{3}^{0} \mathrm{d}A=\rho g\frac{0!1!0!}{(0+1+0+2)!}2A_{e}=\rho g\frac{A_{e}}{3}\\\int_{A_{e}}N_{3}\rho g \mathrm{d}A=\rho g\int_{A_{e}}N_{1}^{0}N_{2}^{0}N_{3}^{1} \mathrm{d}A=\rho g\frac{0!0!1!}{(0+0+1+2)!}2A_{e}=\rho g\frac{A_{e}}{3}
\]
代入得
\[\int_{A_e}[N]^\mathrm{T}\{b\}t \mathrm{d}A=-\frac t3\begin{Bmatrix}0\\\rho gA_e\\0\\\rho gA_e\\0\\\rho gA_e\end{Bmatrix}
\]
可以發現，單元自重在各節點是平均分配的。
表面力
如圖 所示的單元作用有大小為\(q\)的均布載荷，該均布載荷作用在單元的邊 2-3 上，並
且與 \(x\) 軸的夾角為\(\theta\) 。因此，表面力向量可以表示為\(\{t\}=\{-q\cos\theta,-q\sin\theta\}^\mathrm{T}\circ\)
則表面力可以表示為
\[\int\limits_{L_{n}}[N]^{\mathrm{T}}\{t\}t \mathrm{d}l=\int\limits_{L_{23}}\begin{bmatrix}0&0\\0&0\\N_2&0\\0&N_2\\N_3&0\\0&N_3\end{bmatrix}\begin{Bmatrix}-q\cos\theta\\-q\sin\theta\\-q\sin\theta\end{Bmatrix}t \mathrm{d}l=t\int\limits_{L_{23}}\begin{Bmatrix}0\\0\\-N_2q\cos\theta\\\\-N_2q\sin\theta\\\\-N_3q\cos\theta\\\\-N_3q\sin\theta\end{Bmatrix}\mathrm{d}l
\]
注意,在邊 2-3 上 \(N_1=0\)
採用式(8.18)給出的沿三角形邊長的積分公式來計算沿長度方向的積分。運用這個公式上式可以表示為
\[\int\limits_{L_{e}}[N]^{\mathrm{T}}\{t\}t \mathrm{d}l=t\begin{Bmatrix}0\\0\\-q\cos\theta L_{2-3}/2\\-q\sin\theta L_{2-3}/2\\-q\cos\theta L_{2-3}/2\\-q\sin\theta L_{2-3}/2\end{Bmatrix}
\]
從上式可以看出,節點2和節點3平均分配了作用在它們之間的均布載荷\(qL_{2-3}\)。
集中力
% 檔案：CST_COARSE_MESH_DATA.m
%
% 下列變數被宣告為全域性變數，以便被構成程式的所有函式（M檔案）使用
%
global nnd nel nne nodof eldof n % 全域性變數宣告
global geom connec dee nf Nodal_loads % 包括節點數、元素數、節點每元素數、每節點自由度數等
format short e % 設定MATLAB的輸出格式為科學計數法
% 定義模型引數
nnd = 21; % 節點數
nel = 24; % 元素數
nne = 3; % 每個元素的節點數
nodof = 2; % 每個節點的自由度數
eldof = nne * nodof; % 每個元素的自由度數
% 節點座標x和y
geom = zeros(nnd, 3); % 初始化幾何矩陣
% ... 此處省略了具體的節點座標賦值 ...
% 元素連線資訊
connec = zeros(nel, 3); % 初始化連線矩陣
# ... 此處省略了具體的元素連線資訊 ...
% 材料屬性
E = 200000.; % 彈性模量，單位MPa
vu = 0.3; % 泊松比
thick = 5.; % 梁厚度，單位mm
% 形成平面應力的彈性矩陣
dee = formdsig(E, vu); % 呼叫formdsig函式，該函式未在程式碼中給出
% 邊界條件
nf = ones(nnd, nodof); % 初始化自由度矩陣
% ... 此處省略了具體的邊界條件賦值 ...
% 計算自由度的總數
n = 0;
for i = 1:nnd
for j = 1:nodof
if nf(i, j) ~= 0
n = n + 1;
nf(i, j) = n;
end
end
end
% 載荷
Nodal_loads = zeros(nnd, 2); % 初始化節點載荷矩陣
% ... 此處省略了具體的載荷賦值 ...
%%%%%%%%%%%% 輸入結束 %%%%%%%%%%%%%