實(shí)現(xiàn)原理：

1.記錄格點(diǎn)位置

對三個(gè)方向嵌套三層循環(huán)，生成27節(jié)點(diǎn)的網(wǎng)格，使用計(jì)數(shù)器判斷處于單元的第幾個(gè)點(diǎn)，使用oklst和計(jì)數(shù)器比較，判斷該節(jié)點(diǎn)是否屬于20節(jié)點(diǎn)單元的節(jié)點(diǎn)。如果是，則將該點(diǎn)的坐標(biāo)存儲于coordinates矩陣中，如果不是，則不記錄該點(diǎn)位置，并將該點(diǎn)坐標(biāo)編號記錄在fakeindexlst中，用于后續(xù)修正格點(diǎn)編號。

2.將格點(diǎn)編號與單元關(guān)聯(lián)

因?yàn)?-20格點(diǎn)在每個(gè)單元的位置是固定的，因此其編號在全局格點(diǎn)中有周期性。實(shí)現(xiàn)方法是生成27節(jié)點(diǎn)的單元網(wǎng)格，然后利用編號的周期性取需要的20節(jié)點(diǎn)。值得注意的是，這里的單元內(nèi)格點(diǎn)編號和書上的20節(jié)點(diǎn)編號順序并不相同，但我的程序與這里的格點(diǎn)編號是對應(yīng)的。

3.修正格點(diǎn)編號

這樣雖然形成了對應(yīng)關(guān)系，但格點(diǎn)編號卻是27節(jié)點(diǎn)的編號。因此最后一步是把格點(diǎn)編號修正為20節(jié)點(diǎn)，具體而言，是減去小于該節(jié)點(diǎn)編號的假節(jié)點(diǎn)個(gè)數(shù)，例如20號節(jié)點(diǎn)在修正前是第27號節(jié)點(diǎn)，減去27之前的7個(gè)假節(jié)點(diǎn)，修正為第20號。

測試：

1單元20節(jié)點(diǎn)的格點(diǎn)

?

Y方向有兩個(gè)單元，XZ方向有一個(gè)單元生成的格點(diǎn)。

二、單元分析

處理數(shù)據(jù)的方式是先單元分析逐個(gè)計(jì)算剛度矩陣，再整體處理約束方程。

2.1 問題描述：

三維彈性介質(zhì)的本構(gòu)方程

三維彈性介質(zhì)的幾何方程

無重力情況三維彈性介質(zhì)的平衡方程

所以

2.2 基于弱形式的有限元單元分析解

形函數(shù)選取20節(jié)點(diǎn)的形函數(shù)。

2.3 變換到自然坐標(biāo)系計(jì)算剛度矩陣

導(dǎo)數(shù)的坐標(biāo)變換：

是單元內(nèi)20個(gè)格點(diǎn)在的全局坐標(biāo)向量。

因此某單元的剛度矩陣為：

2.4 數(shù)值積分

使用三個(gè)高斯點(diǎn)，高斯點(diǎn)及對應(yīng)權(quán)重。

三、組集

3.1 組集剛度矩陣

組集時(shí)建立節(jié)點(diǎn)局部自由度索引向量和節(jié)點(diǎn)總體自由度索引向量的關(guān)系，提取單元的20節(jié)點(diǎn)全局格點(diǎn)編號，計(jì)算每個(gè)單元的節(jié)點(diǎn)總體自由度索引向量。對m號節(jié)點(diǎn)，其總體自由度索引是。對每個(gè)單元的節(jié)點(diǎn)總體自由度索引向量，使用雙重循環(huán)將單元的剛度矩陣的值賦給總體剛度矩陣的對應(yīng)位置。

四、邊界約束的考慮

對剛度矩陣和廣義力的修正

記錄有邊界約束的格點(diǎn)全局自由度索引和位移值，然后將總體剛度矩陣對應(yīng)位置的索引處改為1，該索引所在行其他元素設(shè)為0，將力的對應(yīng)位置改為位移值。即

五、求解方程

一些技術(shù)細(xì)節(jié)：

避免直接求解逆矩陣，而是使用高斯消去法求解。MATLAB提供了這一功能。

使用稀疏矩陣存儲總體剛度矩陣。MATLAB提供了這一功能。

結(jié)果展示我怕被人拿去交自己的作業(yè)，就不貼上來了。隨便show幾張圖

結(jié)束。