2023年8月9日 星期三

有限元素法入門簡介(21) 有限元素法的數學模型


在開始使用CAE前必須確實學好的有限元素法(11)


簡直就像是循環賽?!使用剛性矩陣來思考


做出加入節點間彈簧常數的矩陣,看起來就像循環賽一樣!這也是有限元素法分析模擬的基礎。


黑箱的內部就是這個!


透過這一系列的文字,並且使用了免費軟體來說明了很多實際模擬的知識。不知道各位讀者是否有實際體驗過了嗎?如果還沒有親自操作過的話,請務必參照前文跑看看。


雖說好漢不提當年勇,但可以使用免費軟體來跑有限元素法的分析模擬,在20年前是簡直像是做夢一樣。而且有限元素法模型可以直接用電腦顯現出來...。以前的話,如果要檢查模型是否有問題,還要特別用繪圖機將網格圖印出來看才行。而且如果節點的座標值的十位數錯亂,造成繪圖機亂印。結果用掉了繪圖機一整捲的紙,因此還挨了罵。


言歸正傳,本系列的文章也將走入尾聲。也該具體說明一下有限元素法模擬的原理了。


以下就正式開始說明。


請試著想成是彈簧的大集合


到目前為止,為了簡明易懂,所以說明上關於數學方面或是材料力學方面不太適切之處。這都是為了讓大家對於有限元素法的架構有個概念,所以就請大家手下留情、不要太過追究。


以下使用之前也介紹過的簡單模型為例,來說明有限元素法的程式到底在做甚麼。


其實在有限元素法中,有專門用來表示各個元素特性的數學式子。這些數學式子雖然複雜,但簡單來說,請想成這些式子都在表現元素的「硬度」就好了。元素是由節點所構成的,每個節點又有六個自由度。如果這樣講有點看不懂的話,請務必回頭看看前文的說明。


另外,以下的說明也必須考慮3D立體的狀態才行,但這次為了容易理解還是用2D平面的方式來解釋。


藉由表現元素特性的數學式,可以定義出節點之間的「硬度」。也可以把這個「硬度」稱呼為「剛性」。也可以試著用「彈簧」來做為模式化的表示。



圖1 用彈簧來表現有限元素。


原本這裡面也需要控制六個自由度的彈簧,但這次也加以省略,只考慮模型的平面上、作用在XY平面的彈簧而已。這樣有限元素模型,就變成了節點與彈簧的聚合體。雖是我自己的亂想像,總覺得上面這個彈簧圖看起來很像是金屬做的刷子...。


然後先來思考節點1的部分。節點1和節點2、5、6藉由彈簧而彼此相連著。所以當節點1想要動時,就被連在節點1上的彈簧所限制住。而這個限制的強度,就是由彈簧的強度、換句話說就是彈簧常數所決定。


說到彈簧常數,對,就會讓人想到「虎克定律」。虎克定律只有一條公式而已。那麼該如何使用虎克定律來表示有限元素法中節點間彼此相連的複雜狀況呢?


這時候,就是矩陣登場的時候了。


歡迎來到矩陣的世界


以下的部分就稍稍有點小困難了。因為虎克定律只有一條公式,該怎麼用來表現連接節點與節點之間的各個彈簧呢?這裡就需要矩陣的思考方式了。雖然有點小複雜,還是請耐著性子看下去喔。我想這裏多少可以讓大家看到工程與數學相結合瞬間的美麗世界。


為了整理出節點與節點之間的關係,用來表示彼此要全部都對戰過的循環圈賽表格就很適合。既然一個節點有六個自由度,則讓節點1在X平移方向的自由度稱為1X、節點1在Y平移方向的自由度稱為1Y,然後就能使用這種「節點編號+自由度方向」的形式,做出一個類似巡迴賽對戰表的東西。至於旋轉方向的自由度則是以類似θx的形式來表現。


在這次的範例中,由於共有八個節點,所以就會有6*8=48個的項目出現。接下來就要在這48個項目的「循環賽對戰表」中填入彈簧常數。本來這個填寫過程會更加困難的,但這裡就盡量用最簡略的方式來思考彈簧常數即可。經過這樣大略整理之後,就完成了填入節點之間彈簧常數的表格了。我們就把這個表稱為「剛性矩陣」。循環圈賽表中,橫列和縱列的隊伍數是相同的。在剛性矩陣中,也是縱橫的數量、換句話說就是行與列的數量是一樣的。以圖1為例,就需要48*48行列的矩陣。


以上這樣說明,可能會惹怒原本就熟悉有限元素法、或是具有寫作有限元素法程式技巧的朋友,但這樣的說明是基於建立概念為最優先的方針下做出的選擇,尚請見諒。想要知道更詳細正確內容的人,請務必找有限元素法的專門書籍來研究。



圖2 很像循環圈賽勝敗關係表的剛性矩陣


再來,請將這個剛性矩陣整體想成虎克定律中的彈簧常數。不過,在這之前,要先請大家複習一下矩陣的運算方式。矩陣運算照理說應該在國中或高中數學裡有教過才對。請務必一面回憶其內容、一面看下去。



圖3 矩陣的運算


請看一下圖3最上方的矩陣乘法計算。2X2的矩陣和2X2的矩陣相乘,結果還是2X2的矩陣。


至於第二排的矩陣計算,則是2X2與2X1矩陣相乘,結果則是2X1的矩陣。


那麼第三排的矩陣運算會有甚麼結果?結論是2X2和3X1的矩陣無法相乘,因為行和列的數目對不起來。矩陣運算時,被乘矩陣的列數要與相乘矩陣的行數要相同才行。


那麼就請大家把矩陣運算的原則記在心中,然後來看看如何用矩陣來表現虎克定律吧。


圖4 用矩陣表示虎克定律


接下來我們就要看看相當於施力的「f」矩陣。在有限元素法中該怎麼表現負載這個東西呢?就像之前在使用免費軟體實作一樣,施加在有限元素模型上的負載,就相當於施加在節點某個方向上、某個大小的力。因為同時牽涉到大小與方向,就是一種向量(Vector)。


這時候請回想一下剛性矩陣的構成項目。其各個項目是以「節點編號+方向」的方式來表示。在對應的項目空格裡,如果填入相當於大小的數字,就能定義出成為向量的力。像圖1這樣的例子,其力量矩陣就是48x1的矩陣。


然後我們用同樣的方式來看看伸長量「u」這個矩陣。「伸長量」在有限元素法中也可以說是變形量,換句話說就是「要求解的答案」。伸長量和力同樣,就是節點在某個方向能移動多少的量。所以這個變形量其實也是一種向量,而在圖1這個例子中,也是48x1矩陣。如果「某個節點」在「某個方向」上不能動的話,也就是「其變形量為零」,就是一種拘束條件了。


因此就像圖4所表示的矩陣相乘一樣,因為行列數相符,就知道這個f=Kxu的矩陣乘法是成立的。


剛性矩陣就是從有限元素模型中組合起來的。而負載也是透過負載條件設定、拘束也是利用拘束條件來設定。所以就把負載輸入48X1矩陣中應該所在方向的空格中。而拘束條件也依照同樣的方法在必須要固定的空格中設為0。


在這裡我們來順便看一下某一個給設計者使用的分析模擬軟體(譯註:顯然是Ansys Workbench)的負載條件設定指令。是不是有很多看起來很便利的指令排在一起?這裡想要請大家理解的是,不管你使用這個指令表中的哪一個指令,其實結果上都還是等於在節點上定義向量而已。


圖5 不管施加怎樣的負載條件,結果都會被置換成設定向量。


至於拘束也是一樣,同樣有個類似的指令表、排著漂漂亮亮的指令。其實不管使用這裡面的哪個指令,結果也還是將某個節點的某個方向的變形設為零而已。


圖6 不管是執行怎樣的拘束指令結果還是會置換成把節點變形量設為零


以下是筆者個人的意見,特別是將負載或拘束置換成設計者容易理解的語言時要小心注意。因為在讓設計者容易理解的同時,也可能造成他們會在不明白自己施加的負載與拘束條件之意義下使用這些指令。


本來定義負載或拘束條件的意義,是在於將負載或拘束定要到「哪個節點」的「哪個自由度」上。換句話說,設計者自己必須要理解這些用讓他們容易明白的語言表現出來的負載或拘束是「怎麼作用的?」或「施加在哪個自由度上?」才行。


為了理解自己打算使用的「負載」或「拘束」指令是怎麼樣作用、施加到哪個自由度上,筆者是強烈建議要先拿小的、簡單的模型來檢查確認看看,再正式使用下去會比較好。


然後才是計算出變形量


如果到這裡為止都能理解清楚的話,再來就只剩求出變形量了。也就是計算出「u」這個矩陣。如果像虎克定律一樣只有一條公式一樣,只要在等號的兩邊各除以「K」,就能簡單計算出「u」了;但是在矩陣運算中,並不是這麼簡單。


在矩陣方程式中,除法運算是要先求出「反矩陣」、然後再於等號左右兩側各乘上反矩陣才能成立。這裡是不是看到了很多令人懷念的術語到處亂竄?好比說,「反矩陣」就是一個「乘上原本的矩陣,使其成為『單位矩陣』」的矩陣。而「單位矩陣」則是一個「對角線的元素都是1,其他通通都是零」的矩陣。


圖7 求出剛性矩陣的反矩陣而藉此求出變形量


到這裡應該就能看的出來,其實有限元素法分析模擬的過程中,幾乎絕大部分的時間都花在組合出剛性矩陣與求出其反矩陣上了。


以上就是用來解釋有限元素法的計算原理的一大段長文,不知道大家是否能夠理解?


而這一篇希望請各位讀者務必要記住的內容則是:


「加在結構上的負載」以及「限制結構不能動的拘束」,結果都是在決定節點的自由度。


到這裡為止,已經解釋了要去閱讀有限元素法專門書籍前所必須知道的基本中的基本知識。


不過,各位從事設計工作的朋友要真正下去做分析模擬的話,還是有些需要補充的內容。好比說,到底要怎麼設定加諸於實物的拘束條件或負載條件才對?怎樣的問題才適用於靜態分析模擬的範圍?要如何理解頻譜分析模擬?所以想要寫出來的東西還有很多很多,但是這些實際的分析模擬技巧,就已經超過本文的範圍,如果有機會再另外來介紹這些主題。


===


相關系列文章:

沒有留言:

張貼留言