摘要: 工程數學為大專校院之理工科系必修的重要學科, 它承接大一的微積分且為大二或大三許多理論學科的基礎。 根據筆者過去授課經驗發現, 許多學生因認為該學科較艱澀、 不知如何將數學理論與其他專業科目知識連結再應用於實際工程問題、 以及未列入國家考試科目而降低學習意願。以真實工程問題串連必修課程為教育部在推行「大專校院新工程教育方法實驗與建構計畫」 的重要主軸, 筆者以海陸域風力發電為主題獲得此計畫, 並在執行時連結土木系三門必修學科知識 --- 工程數學、 土壤力學、 與基礎工程。 文章由土木系學生在土壤力學與材料力學常見的學習紊亂開始, 闡述不同學科因需求與定義的差異而導致在應力分析之莫爾圓公式的差別。 對照力學教科書上基於自由體力平衡所得的莫爾圓公式, 由作者基於旋轉矩陣所產生之相似矩陣的特徵值性質之數學證明與前者結果一致。 藉由對比本文中在工程數學部分的一系列證明與力學教科書上的公式, 讀者可強化鏈結數學與力學的關聯性與建立相關之物理意義。 數學證明結果經 MATLAB 驗證並以圖形闡釋, 而矩陣特徵值與莫爾圓的觀念可應用在基礎工程的河海堤加載應力分析上、 並以軟體 FLAC 來解析。 本文基於力學莫爾圓問題, 透過工程數學矩陣特徵值來證明其相關性質且釐清常見的學習紊亂, 藉此說明力學為數學之延伸、 並可應用在真實的工程問題上。

關鍵字: 工程數學、 工程力學、 矩陣特徵值問題、莫爾圓、 新工程教育方法實驗與建構計畫。

1. 前言

工程數學為工學院學生的必修科目, 其內容通常包含： 微分方程式、 偏微分方程式、 拉普拉斯變換 (Laplace transform)、 傅立葉分析、 向量分析、 線性代數、 複變分析等章節, 而這些知識可做為學生未來追求進階科目所需的基礎。 以土木工程學系為例, 彈性力學、 結構動力學、 土壤動力學為研究所的核心課程, 而此些科目需構築於工程數學之上。 然而, 筆者於大學部工程數學課堂中所做調查指出 (圖 1), 修課學生普遍不知如何將工程數學的知識應用於其他學科之上, 且認為數學與力學學科知識是彼此獨立並無法交互串聯和應用。 再者, 土木工程學系畢業生若未來參加技師高考, 工程數學並非考試科目, 故學生對此科目的學習意願低落。 陳正宗 (2007) 指出, 如何讓學生不畏懼與有效學習工程數學一直是工程教育的重要課題。 因此, 本文透過工程領域的視角帶領讀者認識工程知識背後的數學架構與工程數學的應用。

教育部提出「以真實工程問題串連必修課程」的概念推動「大專校院新工程教育方法實驗與建構計畫」 (擷取自新工程教育方法實驗與建構計畫網站)。 2019年, 作者基於過往於美國離岸石油與天然氣業界工作經驗, 帶領逢甲大學土木工程學系以 「智慧營建在海陸域風力發電設施之理論設計與實務應用」為題獲得 B 類計畫補助; 更於 2021 年進一步地結合逢甲大學聚焦的研究領域, 以「智慧營建 深耕土木 永續防災」 為主體獲得 A 類計畫補助。 我國政府推動可再生能源產業的發展, 並預計於 2035 年前達到 10GW 的離岸風場建設, 故畢業生進入離岸工程領域是可預見的。 為了使學生更全面性地認識離岸風力發電產業, 作者於教材發展及課堂講解中, 提供全面性的規畫。

根據作者授課經驗, 由於不同力學教科書中 (如土壤力學或材料力學) 的應力分析章節彼此採用不同的符號系統與公式表示式, 這使得學生易發生學習混亂。 不同的教科書中所顯示的應力分析公式存在著些許差異, 此導致學生需背誦不同的公式。 再者, 應力分析在工程力學科目採用基於力平衡條件 (Equilibrium Condition) 所推導的莫爾圓 (Mohr's circle) 公式解析、 而在工程數學的範疇裡則屬矩陣的特徵值 (Eigenvalue) 問題。 然而對於上述章節而言, 力學與數學間的知識並未被良好地統整, 此使得學習者認為上述二法是相互獨立的觀念並難以串聯學科間的知識 (陳正宗等人, 2019; 陳正宗等人, 2020)。 此外, 實際工程專案的內容是高度複雜的, 故培養學生具有操作相關數值模擬軟體的能力是必須的。 然而, 修課生將理論課程知識應用至數值軟體的能力較為薄弱。 因此, 本研究將基於力學應力分析的觀點應用工程數學中矩陣分析理論進行驗證, 並將所得計算結果經由 MATLAB 數值運算再次進行比較。 最後, 本文所得之各項理論結果將被連結至真實的工程問題並結合 FLAC 計算結果進行全面性地成果展現。 綜合上述, 以下為工程科系學生在應力分析上常見的學習困難：

1.1. 定義不同導致學習混亂

於力學分析領域中, 「應力轉換公式」與「應力莫爾圓」的觀念是極為關鍵的; 而在「土壤力學」與「材料力學」中, 此些概念係基於「力平衡」的觀念進行推演。 然而, 不同學科間所採用的張應力、壓應力、剪應力方向之定義不完全相等, 進而導致所對應的應力轉換公式不一致, 可參考表 1。 對於初學者而言, 此常導致發生混淆的現象。

1.2. 無法串聯學科間知識

為了分析微小單元體的受力情況, 「應力莫爾圓」與「應力轉換公式」可分別被應用並進一步地獲得此受力單元的主應力 (principal stress) 與主軸方向。 另一方面, 此議題亦可藉由分析此受力單元的應力矩陣獲得相同的結果, 而在工程數學矩陣分析章節中, 特徵值與特徵向量 (Eigenvector) 問題恰可用來解析此類問題。 因此, 力學與數學間的關係未被直接地連結與建立。 為了強化學習者鏈結數學與物理力學的關聯性, 本文藉由不同的學科知識闡述微小單元的受力情況並將所得結果相互連結。

1.3. 缺乏與數值軟體的連結

課堂間教師所傳授的專業知識係基於數學模型與力學理論, 而學生操作數值軟體係為學習成果的展現。 本研究分別採用應力莫爾圓與矩陣分析的觀點描述單元體的受力行為, 並將理論分析結果與 MATLAB 運算結果進行交互驗證。 此外, 此些理論結果亦被應用至實際的工程問題並進一步地使用高階數值模擬軟體FLAC進行闡述。 文中藉由河、 海堤或水壩的加載問題指出, 受力單元將發生主軸旋轉的現象, 而此亦呼應相關力學理論的分析結果。

2. 土壤力學與材料力學應力分析

2.1. 多套符號系統定義與併行

由於不同學科的需求與考量不同, 在座標系統與應力正負方向上的定義也有不同。 舉例來說, 土壤力學探討地表下土壤的受力行為, 故縱座標以向下為正。 再者, 又因土石材料多受壓力而無法承受張力, 故定義壓應力為正 (圖 4); 而剪應力的定義更有順時針為正或逆時針為正之區別, 由表 1可見, 現行力學教科書所採用的應力轉換公式之定義相當繁雜。 然而, 此尚未能完全地表達符號系統定義之混亂程度。 由於力學計算需搭配合適的座標系統方可進行, 而座標系統與應力正負方向之定義將進一步地使學習者困惑, 可參考圖 2 至圖 4。

2.2. 應力轉換公式

基於圖 2 至圖 4 所繪製的自由體圖 (free body diagram) 並採用力平衡條件, 則可獲得相關應力轉換公式, 式 \eqref{1} 或式 \eqref{2} :

\begin{align} &\left\{\begin{array}{l} \sigma_\theta =\dfrac{\sigma_x+\sigma_z}{2}+\dfrac{\sigma_x-\sigma_z}{2}\cos 2\theta +\tau_{xz}\sin 2\theta, \\[7pt] \tau_\theta =-\dfrac{\sigma_x-\sigma_z}{2}\sin 2\theta +\tau_{x z}\cos 2\theta, \end{array}\right.\label{1}\\[7pt] &\left\{\begin{array}{l} \sigma_\theta =\dfrac{\sigma_x+\sigma_z}{2}+\dfrac{\sigma_x-\sigma_z}{2}\cos 2\theta -\tau_{xz}\sin 2\theta,\\[7pt] \tau_\theta =\dfrac{\sigma_x-\sigma_z}{2}\sin 2\theta +\tau_{x z}\cos 2\theta. \end{array}\right.\label{2} \end{align}

值得注意的是, 圖 2 至圖 4 顯示三種不同的定義型式, 然而其所獲得的應力轉換公式之表示式僅存在兩種。 利用圖 2 與圖 4 所獲得的應力轉換公式為式 \eqref{1}, 而基於圖 3 之結果則為式 \eqref{2}。

2.3. 應力莫爾圓

利用應力轉換公式可獲得應力莫爾圓表示式, 而無論採式 \eqref{1} 或式 \eqref{2} 所得之結果皆為：

\begin{align} \Big(\sigma_\theta -\frac{\sigma_x+\sigma_z}{2}\Big)^2+\tau_\theta^2=\Big(\frac{\sigma_x-\sigma_z}{2}\Big)^2+\tau_{xz}^2.\label{3} \end{align}

在這之中, 雖然應力莫爾圓的數學表示式皆為式 \eqref{3}, 但於繪製圖形時其所採用的座標系統亦存在著不同, 可見圖 5 與圖 6。 藉由此些圖形可見, 兩者之縱軸方向恰好相反, 圖 5 係為縱軸向下為正、而後者為縱軸向上為正。 此外, 利用應力莫爾圓可知, 若微小單元體處於平面應變 (plane strain) 狀態下, 則其所對應的最大主應力 $(\sigma_1)$ 與最小主應力 $(\sigma_3)$ 大小可表示為式 \eqref{4} :

\begin{align} \left\{\begin{array}{l} \sigma_1=\dfrac{\sigma_x+\sigma_z}{2}+\sqrt{\Big(\dfrac{\sigma_x-\sigma_z}{2}\Big)^2+\tau_{xz}^2}, \\[7pt] \sigma_3=\dfrac{\sigma_x+\sigma_z}{2}-\sqrt{\Big(\dfrac{\sigma_x-\sigma_z}{2}\Big)^2+\tau_{xz}^2}.\end{array}\right.\label{4} \end{align}

3. 工程數學觀點

3.1. 應力矩陣

考慮在平面應變條件下之土體單元受力狀態為例, 如圖 7 所示, 其對應的二維應力矩陣 (stress matrix) 可表示為：

\begin{align} \mathbf{\tau} =\left[\begin{array}{ccc}\sigma_x&~&\tau_{xz}\\ \tau_{zx}&&\sigma_z \end{array}\right]. \label{5} \end{align}

再者, 若土體單元受力作用後仍保持靜力平衡, 則式 \eqref{6} 須被滿足 :

\begin{align} \tau_{xz}=\tau_{zx}. \label{6} \end{align}

由式 \eqref{5} 與式 \eqref{6} 可知, 此微小單元所對應的應力矩陣將屬於實對稱矩陣 (real symmetric matrix)。

3.2. 旋轉矩陣

旋轉矩陣 (rotation matrix) 可用來描述一位置向量發生旋轉的現象, 而此矩陣的類型又可分為兩種形式, 分別為座標系統的旋轉與位置向量的旋轉 (Arfken, 1985; Wolfram MathWorld, n.d.; 陳正宗等人, 2019)。 考量 $x$-$y$ 平面上存在一位置向量 $(x,y)$, 其繞座標原點逆時針旋轉 $(\theta_p)$ 後所對應的位置向量為 $(x',y')$, 如圖 8所示, 則兩位置向量可利用式 \eqref{7} 表示 :

\begin{align} \left[\begin{array}{c} x'\\ y'\end{array}\right]={\mathbf R}_{xy} \left[\begin{array}{c} x\\ y\end{array}\right]=\left[\begin{array}{ccc}\cos\theta_p&~&-\sin\theta_p\\ \sin\theta_p&&\cos\theta_p\end{array}\right]\left[\begin{array}{c} x\\ y\end{array}\right].\label{7} \end{align}

另一方面, 若考量原座標系統以原點為中心逆時針旋轉 $(\theta_O)$ 後, 其所對應的位置向量為 $(x',y')$, 圖 9, 則兩位置向量可利用式 \eqref{8} 表示 :

\begin{align} \left[\begin{array}{c}\overline x'\\ \overline y'\end{array}\right] =\overline {\mathbf R}_{xy} \left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{ccc}\cos\theta_O&~&\sin\theta_O\\ -\sin\theta_O&&\cos\theta_O\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right].\label{8} \end{align}

再者, 大地工程領域常採用縱座標軸 ($z$) 向下的座標系統, 即 $z$-$x$ 平面, 故若考量有一位置向量 $(x,z)$ 繞座標原點逆時針旋轉 $(\theta_p)$ 後所對應的位置向量為 $(x',z')$, 如圖 10 所示, 則兩位置向量可利用式 \eqref{9} 表示 :

\begin{align} \left[\begin{array}{c} x'\\ z'\end{array}\right]={\mathbf R}_{xz} \left[\begin{array}{c}x\\ z\end{array}\right]=\left[\begin{array}{ccc}\cos\theta_p&~&\sin\theta_p\\ -\sin\theta_p&&\cos\theta_p\end{array}\right] \left[\begin{array}{c}x\\ z\end{array}\right].\label{9} \end{align}

另一方面, 若考量原座標系統以原點為中心逆時針旋轉 $(\theta_O)$ 後, 其所對應的位置向量為 $(x',z')$, 圖 11, 則兩位置向量可利用式 \eqref{10} 表示 :

\begin{align} \left[\begin{array}{c}\overline x'\\ \overline z'\end{array}\right]={\overline{\mathbf R}}_{xz} \left[\begin{array}{c}x\\ z\end{array}\right]=\left[\begin{array}{ccc} \cos\theta_O&~&-\sin\theta_O\\ \sin\theta_O&~&\cos\theta_O\end{array}\right]\left[\begin{array}{c}x\\ z\end{array}\right].\label{10} \end{align}

值得注意的是, 對於位置向量繞原點旋轉之旋轉矩陣受到座標系統的影響, 兩者 (式 \eqref{7} 與 \eqref{9}) 的表示式將有所不同; 而同樣的現象亦發生於座標系統旋轉所對應的旋轉矩陣 (式 \eqref{8} 與 \eqref{10})。 另一方面, 本文主要探討並應用的議題為工程數學、 土壤力學與基礎工程學科知識的連結, 而土壤力學和基礎工程領域中多採用 $x$-$z$ 平面 (圖 7) 進行探討。 綜上所述, 本文後續論述將以式 \eqref{9} 與式 \eqref{10} 進行。

3.3. 利用應力矩陣與旋轉矩陣獲得應力轉換公式

微小單元受力作用後所對應的應力矩陣之分量 (component) 會受到座標系統的影響, 如座標系統旋轉後所對應的應力分量將有所改變。 另一方面, 土體單元發生旋轉的現象亦將使應力矩陣的元素 (element) 發生變化。 值得注意的是, 利用系列座標系統描述單一土體單元的受力狀態或考量土體單元發生連續地旋轉所對應不同的應力分量之行為, 此些概念可被稱為應力轉換。 於土壤力學或材料力學課程中, 應力轉換公式係利用力平衡的概念進行推演, 然而此議題亦可利用矩陣運算求得。 考量初始土體單元之應力矩陣為式 \eqref{11} :

\begin{align} \mathbf{\tau} =\left[\begin{array}{ccc}\sigma_x&~&\tau_{xz}\\ {\tau}_{xz}&~&\sigma_z\end{array}\right],\label{11} \end{align}

考慮以點位旋轉矩陣 (${\mathbf R}_{xz}$) 對應力矩陣 ($\mathbf\tau $) 所形成的相似矩陣 (${\mathbf\tau}_\theta $), 則可表示為式 \eqref{12} (Das, 2008; 陳正宗等人, 2019):

\begin{align} \mathbf{\tau}_\theta =\,&\left[\begin{array}{ccc}\tau_{\theta_{11}}&~&\tau_{\theta_{12}}\\ {\tau}_{\theta_{21}}&~&\tau_{\theta_{22}}\end{array}\right] =\mathbf{R}_{xz} \mathbf{\tau} \mathbf{R}_{xz}^T,\label{12}\\ {\hbox{其中：}} \tau_{\theta_{11}}=\,&\frac{(\sigma_x+\sigma_z)}{2}+\frac{(\sigma_x-\sigma_z )}{2}\cos 2\theta_p+\tau_{xz}\sin 2\theta_p,\label{13}\\ \tau_{\theta_{12}}=\,&\tau_{\theta_{21}}=-\frac{(\sigma_x-\sigma_z )}{2}\sin 2\theta_p+\tau_{xz}\cos 2\theta_p,\label{14}\\ \tau_{\theta_{22}}=\,&\frac{(\sigma_x+\sigma_z )}{2}-\frac{(\sigma_x-\sigma_z )}{2}\cos 2\theta_p-\tau_{xz}\sin 2\theta_p.\hskip 2cm~\label{15} \end{align}

值得注意的是, 式 \eqref{13} 與式 \eqref{14} 的表示式即為式 \eqref{1} 所列的應力轉換公式。 另一方面, 若考慮以座標旋轉矩陣 $(\overline {\mathbf R}_{xz})$ 對應力矩陣 ($\mathbf\tau $) 所形成的相似矩陣 (${\mathbf\tau}_{\theta'}$), 則可表示為式 \eqref{16} :

\begin{align} \mathbf{\tau}_\theta' =\,&\left[\begin{array}{ccc}\tau_{\theta_{11}}'&~&\tau_{\theta_{12}}'\\ {\tau}_{\theta_{21}}'&~&\tau_{\theta_{22}}'\end{array}\right] ={\overline{\mathbf{R}}_{xz}} \mathbf{\tau} \overline{\mathbf{R}}_{xz}^T,\label{16}\\ {\hbox{其中：}} \tau_{\theta_{11}}'=\,&\frac{(\sigma_x+\sigma_z)}{2}+\frac{(\sigma_x-\sigma_z )}{2}\cos 2\theta_O-\tau_{xz}\sin 2\theta_O,\label{17}\\ \tau_{\theta_{12}}'=\,&\tau_{\theta_{21}}'=\frac{(\sigma_x-\sigma_z )}{2}\sin 2\theta_O+\tau_{xz}\cos 2\theta_O,\label{18}\\ \tau_{\theta_{22}}'=\,&\frac{(\sigma_x+\sigma_z )}{2}-\frac{(\sigma_x-\sigma_z )}{2}\cos 2\theta_O+\tau_{xz}\sin 2\theta_O.\hskip 2cm~\label{19} \end{align}

在這之中, 式 \eqref{17} 與式 \eqref{18} 之表示式則為式 \eqref{2} 所顯示的應力轉換公式。

3.4. 應力分量的變化軌跡即為應力莫爾圓

若將式 \eqref{12} 與式 \eqref{16}中元素表示式以極座標系統進行表示, 則可進一步地將矩陣運算結果與應力莫爾圓進行連結。 考量以下條件：

\begin{align} r=\,&\sqrt{\Big(\frac{\sigma_x-\sigma_z}{2}\Big)^2+\tau_{xz}^2},\label{20}\\ {\hbox{與}} &\left\{\begin{array}{l}\dfrac{\sigma_x-\sigma_z}{2}=r\cos\alpha\\ \tau_{xz}=r\sin\alpha \end{array}\right.,\label{21} \end{align}

其中, $\alpha$ 為一個角度。 因此, 式 \eqref{12} 及式 \eqref{16} 可利用式 \eqref{20} 與式 \eqref{21} 改寫為式 \eqref{22} 與式 \eqref{23}

\begin{align} {\mathbf\tau}_\theta =\,&\left[\begin{array}{ccc}\frac{(\sigma_x+\sigma_z )}{2}+r\cos(\alpha -2\theta_p )&r\sin(\alpha -2\theta_p) \\ r\sin(\alpha -2\theta_p )&\frac{(\sigma_x+\sigma_z )}{2}-r\cos(\alpha -2\theta_p )\end{array}\right]\nonumber\\ =\,&\left[\begin{array}{ccc}\sigma_c+r\cos t&~&r\sin t\\ r\sin t&~&\sigma_c-r\cos t\end{array}\right],\label{22}\\ {\hbox{及}} {\mathbf\tau}_\theta'=\,&\left[\begin{array}{ccc}\frac{(\sigma_x+\sigma_z )}{2}+r\cos(\alpha +2\theta_O )&~&r\sin(\alpha +2\theta_O)\\ r\sin(\alpha +2\theta_O)&~&\frac{(\sigma_x+\sigma_z )}{2}-r\cos(\alpha +2\theta_O )\end{array}\right]\nonumber\\ =\,&\left[\begin{array}{ccc}\sigma_c+r\cos t'&~&r\sin t'\\ r\sin t'&~&\sigma_c-r\cos t'\end{array}\right],\label{23} \end{align}

其中, $t=\alpha -2\theta_p$ 及 $t'=\alpha +2\theta_O$。 因此, 微小單元的受力情況可顯示於 $\sigma$ -$\tau$ 平面, 即：

\begin{align} &\left\{\begin{array}{l}\sigma =\sigma_c+r\cos\overline t,\\ \tau =r\sin \overline t.\end{array}\right.\label{24}\\ {\hbox{並可得：}} &(\sigma -\sigma_c)^2+\tau^2=r^2.\label{25} \end{align}

因此, 若將土體單元受力條件顯示於 $\sigma $-$\tau$ 平面, 則依據式 \eqref{25} 可知其應力分量的變化軌跡為圓心在 $\sigma$ 軸上的圓形方程式 (圖 12), 且此結果即為應力莫爾圓的數學表示式 (式 \eqref{3})。 再者, 應力矩陣 ($\mathbf\tau $) 可利用點位旋轉矩陣 $({\mathbf R}_{xz})$ 與座標旋轉矩陣 (${\overline {\mathbf R}}_{xz}$) 獲得兩相異的應力轉換公式, 但其於極座標系統下之表示式是一致的 (式 \eqref{25})。 此外, 應力矩陣經由點位旋轉矩陣 (${\mathbf R}_{xz}$) 或座標旋轉矩陣 (${\overline {\mathbf R}}_{xz}$) 運算後所得的應力莫爾圓皆相同, 其圓心與半徑皆不變, 而此結果即為陳正宗等人 (2019) 透過張量分析所得之不變量。

3.5. 應力矩陣之特徵值即為主應力大小

土體單元受特定外力作用後, 其所對應的正向應力與剪應力分量變化行為可藉由應力莫爾圓進行全面性地闡述。 另一方面, 土體單元的受力情況亦可採用應力矩陣進行表示, 然而在不同的座標系統下其所對應的應力分量將會有所不同。 因此, 為了更進一步地連結一特定的應力莫爾圓, 其所對應的應力矩陣分量變化行為, 矩陣運算中的特徵值問題將被採用。 為了計算式 \eqref{11} 所對應的特徵值 ($ \sigma $), 式 \eqref{26} 所列條件須被滿足 :

\begin{align} \left|\begin{array}{ccc}\sigma_x- \sigma &~&\tau_{xz}\\ \tau_{xz}&~&\sigma_z- \sigma\end{array}\right|=0, \label{26} \end{align}

並可得：

\begin{align} \sigma =\frac{(\sigma_x+\sigma_z )\pm \sqrt{(\sigma_x+\sigma_z )^2-4(\sigma_x \sigma_z-\tau_{xz}^2 )}}{2},\label{27} \end{align}

其中, 應力矩陣屬於實對稱矩陣, 因此其所對應的特徵值將屬於實數。 再者, 式 \eqref{27} 可進一步地表示為：

\begin{align} \sigma =\sigma_c\pm r.\label{28} \end{align}

因此, 應力矩陣所對應的特徵值即為此受力單元的主應力大小 (式 \eqref{29})

\begin{align} \left\{\begin{array}{l}\sigma_1=\sigma_c+r,\\ \sigma_3=\sigma_c-r.\end{array}\right.\label{29} \end{align}

3.6. 相同應力莫爾圓上有著相同特徵值

若進一步地利用極座標的分析結果表示微小單元體的應力矩陣：

\begin{align} {\mathbf\tau}_\theta =\left[\begin{array}{ccc} \sigma_c+r\cos(\overline t)&~&r\sin(\overline t)\\ r\sin(\overline t)&~&\sigma_c-r\cos(\overline t)\end{array}\right].\label{30} \end{align}

為了計算式 \eqref{30} 所對應的特徵值, 可利用式 \eqref{31} 進行求解

\begin{align} \left|\begin{array}{ccc}\sigma_c+r\cos(\overline t)- \sigma &~&r\sin(\overline t)\\ r\sin(\overline t)&~&\sigma_c-r\cos(\overline t)- \sigma\end{array}\right|=0,\label{31} \end{align}

並可得

\begin{equation} \begin{aligned} \sigma =\,&\frac{2\sigma_c\pm \sqrt{4\sigma_c^2-4(\sigma_c^2-r^2 )}}{2}\\ =\,&\sigma_c\pm r. \label{32} \end{aligned} \end{equation}

因此, 對比式 \eqref{28} 與式 (32), 顯示於相同的應力莫爾圓上其所對應的應力矩陣皆有相同的特徵值。

3.7. 特徵方向所旋轉的角度恰為應力莫爾圓上旋轉角之半

式 \eqref{30} 所對應的特徵向量可利用式 \eqref{33} 及式 \eqref{34} 進行求解 :

\begin{align} \left[\begin{array}{ccc}\sigma_c+r\cos(\overline t)-(\sigma_c+r)&~&r\sin(\overline t)\\ r\sin(\overline t)&~&\sigma_c-r\cos(\overline t)-(\sigma_c+r)\end{array}\right]\left[\begin{array}{c} x_{11}\\ x_{21}\end{array}\right]={\mathbf 0},\label{33}\\[6pt] \left[\begin{array}{ccc}\sigma_c+r\cos(\overline t)-(\sigma_c-r)&~&r\sin(\overline t)\\ r\sin(\overline t)&~&\sigma_c-r\cos(\overline t)-(\sigma_c-r)\end{array}\right] \left[\begin{array}{c}x_{13}\\ x_{23}\end{array}\right]={\mathbf 0},\label{34} \end{align}

並可得單位化後之特徵向量式 \eqref{35}、式 \eqref{36}, 分別為最大主應力 ($\sigma_1$) 與最小主應力 ($\sigma_3$) 所對應之特徵向量。 此外, 材料單元受力後將存在最大主應力 ($\sigma_1$)、 中主應力 ($\sigma_2$)、 及最小主應力 ($\sigma_3$), 並可分別獲得其主軸方向 ${\mathbf u}_1$、 ${\mathbf u}_2$、 ${\mathbf u}_3$。 基於平面應變分析觀點, 土體單元受力後最大主應力 $(\sigma_1)$ 與最小主應力 $(\sigma_3)$ 的分析是較受到關注的, 而在此條件下中主應力方向 $({\mathbf u}_2)$ 將不受到外力作用的影響其方向固定。

\begin{align} {\mathbf u}_1=\,&\left[\begin{array}{c}u_{11}\\ u_{21}\end{array}\right]=\left[\begin{array}{c}\cos\Big(\frac{\overline t}2\Big)\\[5pt] \sin\Big(\frac{\overline t}2\Big)\end{array}\right],\label{35}\\ {\hbox{及}} {\mathbf u}_3=\,&\left[\begin{array}{c}u_{13}\\ u_{23}\end{array}\right]=\left[\begin{array}{c}-\sin\Big(\frac{\overline t}2\Big)\\[5pt] \cos\Big(\frac{\overline t}2\Big)\end{array}\right].\label{36} \end{align}

式 \eqref{35}與式 \eqref{36} 分析結果指出, 最大主應力與最小主應力所對應的特徵向量可組成旋轉矩陣。 此外, 特徵方向旋轉的角度 ($\overline t/2$) 恰為應力莫爾圓上旋轉角 ($\overline t$) 之半, 圖 13。

3.8. 具備相同特徵向量之莫爾圓性質

若有兩受力單元分別為 $A$ 與 $B$, 則其所對應的應力矩陣可依序被表示為式 \eqref{37} 與式 \eqref{38} :

\begin{align} {\mathbf\tau}_A=\,&\left[\begin{array}{ccc} \sigma_{cA}+r_A\cos ({\overline t}_A)&~&r_A\sin ({\overline t}_A)\\ r_A\sin ({\overline t}_A)&~&\sigma_{cA}-r_A\cos ({\overline t}_A)\end{array}\right],\label{37}\\ {\mathbf\tau}_B=\,&\left[\begin{array}{ccc}\sigma_{cB}+r_B\cos ({\overline t}_B)&~&r_B\sin ({\overline t}_B)\\ r_B\sin ({\overline t}_B)&~&\sigma_{cB}-r_B\cos ({\overline t}_B)\end{array}\right],\label{38} \end{align}

並進一步地考量兩受力單元之特徵向量相同, 即

\begin{align} \overline t={\overline t}_A={\overline t}_B.\label{39} \end{align}

因此, 此二受力單元所對應的單位化特徵向量可以式 \eqref{40} 進行表示 :

\begin{align} \left[\begin{array}{ccc} {\mathbf u}_{1\overline t}&~&{\mathbf u}_{3\overline t} \end{array}\right]= \left[\begin{array}{ccc} \cos\left(\frac{\overline t}2\right)&~&-\sin\left(\frac{\overline t}2\right)\\[5pt] \sin\left(\frac{\overline t}2\right)&~&\cos\left(\frac{\overline t}2\right)\end{array}\right].\label{40} \end{align}

由於 $A$ 與 $B$ 受力單元之應力矩陣具有相同的特徵方向, 故兩者可進行同步對角化：

\begin{align} \left[\begin{array}{ccc} {\mathbf u}_{1\overline t}&~&{\mathbf u}_{3\overline t}\end{array}\right]^{-1} {\mathbf \tau}_A \left[\begin{array}{ccc} {\mathbf u}_{1\overline t}&~&{\mathbf u}_{3\overline t} \end{array}\right]=\,&\left[\begin{array}{ccc} \sigma_{1A}&~&0\\ 0&~&\sigma_{3A}\end{array}\right],\label{41}\\ {\hbox{與}} \left[\begin{array}{ccc}{\mathbf u}_{1\overline t}&~&{\mathbf u}_{3\overline t}\end{array}\right]^{-1} {\mathbf\tau}_B \left[\begin{array}{ccc} {\mathbf u}_{1\overline t}&~&{\mathbf u}_{3\overline t}\end{array}\right]=\,&\left[\begin{array}{ccc}\sigma_{1B}&0\\ 0&\sigma_{3B}\end{array}\right].\label{42} \end{align}

因此, 若進一步地將受力單元 $A$ 與 $B$ 分析結果連結至 $\sigma $-$\tau$ 平面, 則可表示為圖 14。 綜合上述與圖 14, 若兩應力矩陣具備相同之特徵向量, 其對應之莫爾圓角度必然相同。

4. 論例

此章節以土壤力學的應力分析為例, 座標系統與莫爾圓正負方向定義分別如圖 4 與圖 6 所示。 考量一微小矩形單元受力情況如圖 15所示, 則此單元所對應的應力矩陣可表示為式 \eqref{43} :

\begin{align} &{\mathbf\tau} =\left[\begin{array}{ccc} -4&~&-2\\ -2&~&6\end{array}\right],\label{43}\\ {\hbox{並可得}} &\left\{\begin{array}{l}\sigma_c=1.0\\ r=\sqrt{29}\end{array}\right.,\label{44}\\ {\hbox{與}} &\quad \overline t=201.8^\circ.\label{45} \end{align}

依據式 \eqref{29}, 此微小受力單元所對應的最大主應力、最小主應力, 此亦為應力矩陣之特徵值大小, 可顯示為式 \eqref{46} :

\begin{align} &\left\{\begin{array}{l}\sigma_1=1+\sqrt{29}\\ \sigma_3=1-\sqrt{29}\end{array}\right. ,\label{46}\\ {\hbox{而其所對應的特徵向量分別為式 \eqref{47} 與式 \eqref{48}。}} {\mathbf u}_1=\,&\left[\begin{array}{c} \cos\Big(\frac{\overline t}2\Big)\\[7pt] \sin\Big(\frac{\overline t}2\Big)\end{array}\right]=\left[\begin{array}{c}-0.189\\ 0.982\end{array}\right],\label{47}\\[5pt] {\mathbf u}_3=\,&\left[\begin{array}{c} -\sin\Big(\frac{\overline t}2\Big)\\[7pt] \cos\Big(\frac{\overline t}2\Big)\end{array}\right]=\left[\begin{array}{c} -0.982\\ -0.189\end{array}\right].\label{48} \end{align}

另一方面, 微小單元的受力狀況亦可採用極座標表示式進行表示, 而圖 16 與圖 17 即為此受力單元所對應的應力莫爾圓。 圖 16 採用平面原點 (origin of planes) 描述作用面變化的情形, 即應力莫爾圓上點位與平面原點之連線代表著微小單元作用面的方向。 相較於一般教科書, 本文提出應力原點 $O_s$ 的概念, 即該應力點與應力原點的連線即為該面之正向應力方向, 如圖 17 所示。 若進一步地利用 MATLAB 驗證理論分析結果, 其結果被顯示於表 2。 表中所顯示的 8 個應力矩陣皆位於同一個應力莫爾圓上, 而數值計算結果顯示其所對應的特徵值大小皆相同並可良好地符合理論計算結果。 對於表 2 中顯示的第一個應力矩陣, 其特徵向量亦滿足式 \eqref{47} 與式 \eqref{48}。 再者, 特徵向量即顯示作用面的法向量或正向應力的作用方向, 故表 2 中 No. 1 最大主應力所對應的特徵向量將與圖 16 中 $\vec{O_p D}$ 向量垂直, 且與圖 17 中 $\vec{O_s D}$ 向量平行。

5. 基礎工程

5.1. 平面應變

平面應變條件常應用於大地工程領域, 並可用以簡化力學分析模型。 考量一土壤單元於 $x$-$z$ 平面處於平面應變狀態, 則其應力矩陣可表示為式 \eqref{49} :

\begin{align} {\mathbf \tau} =\,&\left[\begin{array}{ccccc}\sigma_x&~&0&~&\tau_{xz}\\ 0&~&\sigma_y&~&0\\ \tau_{xz}&~&0&~&\sigma_z\end{array}\right]=\left[\begin{array}{ccccc}\sigma_x&~&0&~&\tau_{xz}\\ 0&~&\nu(\sigma_x+\sigma_z )&~&0\\ \tau_{xz}&~&0&~&\sigma_z\end{array}\right],\label{49} \end{align}

其中, $\nu$ 為材料泊松比 (Poisson's ratio)。 顯然地, 式 \eqref{49} 所對應的中主應力與所對應的特徵值即為

\begin{align} \left\{\begin{array}{l} \sigma_2=\nu (\sigma_x+\sigma_z) \\[4pt] {\mathbf u}_2=\left[\begin{array}{c}0\\[-3pt] 1\\[-3pt] 0\end{array}\right]\end{array}\right..\label{50} \end{align}

因此, 於平面應變條件下, 可將三維問題化簡為二維問題進行簡化分析。

5.2. 平面應變下的加載問題

建築物通常構築於土壤之上, 對於土壤而言此為外力加載議題。 土壤單元受力分析應屬於三維問題, 然而若基礎的長度遠大於基礎寬度則可利用平面應變將問題簡化分析, 圖 18。 基於彈性力學理論, 水壩 (圖 18) 或條形基礎 (strip footing) 加載問題可簡化為圖 19 進行分析。 圖 19中P點因均佈鉛直應力作用所造成的應力增量可表示為：

\begin{align} \Delta \sigma_x=\,&\frac q\pi [\alpha -\sin\alpha\cos(\alpha +2\delta )],\label{51}\\ \Delta \sigma_z=\,&\frac q\pi [\alpha +\sin\alpha\cos(\alpha +2\delta )],\label{52}\\ {\hbox{與}} \Delta \tau_{xz}=\,&\frac q\pi [\sin\alpha \sin(\alpha +2\delta )].\hskip 4cm ~\label{53} \end{align}

根據式 \eqref{53} 計算結果進行無因次分析後可得圖 20 與利用數值模擬軟體 FLAC 分析結果圖 21。 由圖中可見當土壤受力後造成部分 $x$-$z$ 平面上出現剪應力 ($\tau_{xz}$), 即表示在此區域內土體主應力方向沿y軸 (垂直紙面向外) 出現旋轉的現象。

6. 結論與建議

本研究藉由執行教育部「新工程教育方法實驗與建構計畫」, 串聯計畫中的三門土木工程學系必修學科：「土壤力學」、「工程數學」、與「基礎工程」。 由土壤力學與材料力學的莫爾圓應力分析開始, 進入工程數學中的矩陣單元, 推導證明基於旋轉矩陣所產生之相似矩陣的特徵值與特徵向量之相關性質, 其結果與土壤力學或材料力學教科書上的莫爾圓應力分析中, 基於自由體力平衡所得之公式一致。學習者可藉由本文中在工程數學的一系列證明, 對比土壤力學與材料力學教科書上的公式及其物理意義, 強化鏈結數學與物理力學的關聯性 (圖 22)。

本研究所得各項分析成果, 其結論與重點如下：

I. 利用應力矩陣與點位旋轉型式的旋轉矩陣可獲得式 \eqref{1} 之應力轉換公式, 而此結果與基於力平衡條件所得公式一致。 另一方面, 若採用座標旋轉型式的旋轉矩陣則可獲得式 \eqref{2}, 且此亦與力平衡條件所得公式相同。 此外, 若兩者顯示於 $\sigma $-$\tau$ 平面, 則其軌跡皆為圓形, 而此即為應力莫爾圓。

II. 微小受力單元所對應的應力矩陣, 其特徵值即為主應力大小。 再者, 對於一固定的應力莫爾圓, 其所對應的應力矩陣皆有著相同的特徵值大小。

III. 若將應力矩陣所對應的單位特徵向量依特徵值由大而小排列, 則此即為旋轉矩陣。

IV. 若不同莫爾圓具有相同的旋轉角則具備相同的特徵向量, 且特徵向量的旋轉角恰為應力莫爾圓旋轉角之半。 此外, 若矩陣擁有相同的特徵方向, 則可進行同步對角化。

誌謝

本論文係教育部新工程教育方法實驗與建構計畫「智慧營建在海陸域風力發電設施之理論設計與實務應用」(MOE-NEEMEC-B-108-02)、 「智慧營建 深耕土木 永續防災」 (MOE-NEEMEC-A-110-11)、 科技部專題研究計畫「台灣離岸浮動式風機平台與自升式平台船之錨錠基礎極限承載力分析」 (MOST 110-2221-E-035-029) 與逢甲大學精進教師教學方案 「由新開發教材導入實際工程問題與電腦計算思維來學習工程數學」 (20C52324) 及「結合國家雲端地質資料庫與串聯跨學科知識的土壤力學專題式學習」 (21C52336) 之研究成果, 承蒙教育部、 科技部與逢甲大學補助相關經費與研究人力使本研究得以順利完成, 僅致謝忱。

參考文獻

本文作者紀昭銘任教逢甲大學土木工程學系,林正山投稿時就讀於逢甲大學土木水利工程與建設規劃博士學位學程。