《有限元程序設(shè)計》PPT課件

《《有限元程序設(shè)計》PPT課件》由會員分享，可在線閱讀，更多相關(guān)《《有限元程序設(shè)計》PPT課件（85頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、有 限 元 程 序 設(shè) 計 教 學(xué) 程 序 ： FEP1 李 明 瑞 編 FEP2 李 明 瑞 編 ZFEP 第 一 章 緒 論 有 限 元 程 序 的 基 本 內(nèi) 容 有 限 元 求 解 器 1.1 有 限 元 程 序 的 基 本 內(nèi) 容有 限 元 法 的 解 題 步 驟結(jié) 構(gòu) 離 散 化 為有 限 元 網(wǎng) 格 計 算 單 元 剛 度矩 陣引 入 約 束 條 件 求 應(yīng) 力 等解 方 程 組 建 立 總 剛 度 矩陣 及 外 載 有 限 元 程 序 的 基 本 內(nèi) 容 ：數(shù) 據(jù) 輸 入階 段 有 限 元 矩 陣 的 計算 、 組 集 和 求 解 數(shù) 據(jù) 輸 出階 段前 處 理 處 理 器 后

2、 處 理 一 、 數(shù) 據(jù) 輸 入 階 段 前 處 理 讀 入 和 生 成 數(shù) 據(jù) ， 形 成 有 限 元 網(wǎng) 格 ， 為 有 限 元矩 陣 計 算 作 準(zhǔn) 備 。 （ 數(shù) 據(jù) 或 圖 形 輸 入 ）1. 標(biāo) 題 和 控 制 信 息2. 計 算 存 貯 分 配3. 節(jié) 點 坐 標(biāo) 和 約 束 信 息4. 單 元 信 息5. 材 料 信 息6. 載 荷 信 息 二 、 有 限 元 矩 陣 的 計 算 、 組 集 和 求 解 處 理 器 計 算 插 值 函 數(shù) 矩 陣 N 、 幾 何 矩 陣 B 、 Jacob矩 陣 J 在 高 斯 點 進(jìn) 行 數(shù) 值 積 分 ， 求 得 單 剛 、 單 元 載 荷

3、 組 集 成 總 剛 、 總 載 求 解 三 、 數(shù) 據(jù) 輸 出 階 段 后 處 理輸 出 結(jié) 果 （ 場 變 量 和 場 變 量 導(dǎo) 數(shù) 等 ） ， 打 印 結(jié) 果 。場 變 量 包 括 ： 位 移 、 溫 度 、 流 場 的 速 度 勢 等 ；場 變 量 導(dǎo) 數(shù) 包 括 ： 應(yīng) 力 、 應(yīng) 變 、 熱 流 、 流 場 速 度 勢 等輸 出 形 式 ： 數(shù) 據(jù) 輸 出 和 圖 形 輸 出 1-2 有 限 元 求 解 器帶 寬 法 ： 只 存 儲 總 剛 矩 陣 的 半 帶 寬 內(nèi) 的 元 素 。等 帶 寬 法 ： 每 行 的 帶 寬 相 等 ， 常 采 用 二 維 數(shù) 組 存 儲變 帶 寬

4、法 ： 每 行 的 帶 寬 不 等 ， 常 采 用 一 維 數(shù) 組 存 儲一 、 帶 寬 法變 帶 寬 分 塊 求 解 思 想 ： 一 個 變 量 的 消 元 只 影 響 剛 度 矩 陣 的 一 部 分 元素 ， 較 大 的 節(jié) 點 分 量 方 程 組 可 以 分 成 較 小 的 局 部系 數(shù) 矩 陣 按 節(jié) 點 逐 步 求 解 。 波 前 法 是 另 一 種 分 塊 儲 存 的 變 帶 寬 法 ， 其 消 元 次序 是 按 單 元 編 號 進(jìn) 行 ， 邊 組 集 邊 消 元 。 調(diào) 入 內(nèi) 存 的 單 元 所 保 留 的 節(jié) 點 稱 為 波 前 節(jié) 點 ， 所消 去 的 節(jié) 點 稱 為 消

5、元 節(jié) 點 ， 消 元 節(jié) 點 是 與 以 后 調(diào) 入 的單 元 沒 有 聯(lián) 系 的 節(jié) 點 ， 即 該 節(jié) 點 的 方 程 已 組 集 完 。 波 前 法 的 回 代 按 消 元 節(jié) 點 的 反 序 進(jìn) 行 。 對 內(nèi) 存 的需 求 取 決 與 最 大 波 前 寬 。二 、 波 前 法 V-Fortran 1. Fortran語 言 的 基 本 格 式2. 變 量3. 基 本 語 句4. 子 例 子 程 序5. 函 數(shù) 子 程 序 6. 其 它 功 能 、 模 塊 第 二 章 FEP2 程 序 的 總 體 設(shè) 計 與 輸 入數(shù) 據(jù) FEP2 的 設(shè) 計 任 務(wù) FEP2 的 結(jié) 構(gòu) 設(shè) 計

6、FEP2 的 主 程 序 FEP2 的 主 控 程 序 FEP2 的 數(shù) 據(jù) 輸 入 格 式 與 程 序 實 現(xiàn) 2-1 FEP2 的 設(shè) 計 任 務(wù)1. FEP2的 簡 要 題 目 說 明 FEP2是 一 個 具 有 通 用 性 的 教 學(xué) 程 序 ， 可 用 于 計 算一 般 的 線 性 靜 力 學(xué) 問 題 。 已 設(shè) 計 了 平 面 梁 單 元 與 平 面 3-9節(jié) 點 元 兩 種 單 元 ， 但 留 下 接 口 。2. 支 持 軟 件 與 硬 件 FORTRAN77以 上 編 譯 器 、 各 種 微 機(jī)3. FEP2的 功 能 3. FEP2的 功 能1） 輸 入 文 件 名 由 用

7、戶 自 由 定 義 ， 但 限 制 為 4個 字 符 ， 輸入 文 件 擴(kuò) 展 名 為 “ DAT”， 輸 出 文 件 擴(kuò) 展 名 為 “ OUT”。2） 節(jié) 點 坐 標(biāo) 、 單 元 信 息 等 具 有 線 性 自 動 生 成 功 能 。3） 可 以 處 理 多 種 工 況 、 多 種 類 型 單 元 組 合 結(jié) 構(gòu) 問 題 。4） 可 以 處 理 固 定 約 束 和 指 定 位 移 。5） 采 用 變 帶 寬 存 儲 、 三 角 分 解 法 求 解 剛 度 方 程 。6） 為 多 種 單 元 留 下 接 口 。 2-2 FEP2 的 結(jié) 構(gòu) 設(shè) 計FEP2： 由 形 式 主 程 序 、 主

8、控 程 序 、 功 能 模 塊 組 成 。主 程 序 的 主 要 功 能 ： 定 義 輸 入 、 輸 出 文 件 名 ， 調(diào) 用 主 控程 序 PCONTR主 控 程 序 PCONTR的 主 要 功 能 模 塊 ：1） 內(nèi) 存 分 配 2） 網(wǎng) 格 生 成3） 變 帶 寬 存 儲 設(shè) 置 4） 單 剛 的 形 成 與 組 集5） 剛 度 矩 陣 的 分 解 6） 形 成 節(jié) 點 載 荷7） 形 成 與 組 集 單 元 載 荷 8） 求 解 位 移9） 求 單 元 應(yīng) 力 10） 求 結(jié) 構(gòu) 反 力 FEP2 程 序 框 圖FEP2 （ 主 程 序 ） PCONTR(主 控 程 序 ) SETM

9、EM 檢 查 內(nèi) 存 PROFIL 形 成 變 帶 寬 剛 度 矩 陣 地 址 PFORM (3) 形 成 單 剛 并 組 集 LDLT 總 剛 的 三 角 分 解 GENVEC 形 成 節(jié) 點 載 荷 PFOM (3) 構(gòu) 造 并 組 集 單 元 載 荷 FORBACK 前 消 回 代 求 出 位 移 PFORM (4) 計 算 單 元 應(yīng) 力 PFORM (6) 計 算 結(jié) 構(gòu) 反 力 2-3 FEP2 的 主 程 序一 、 文 件 名 FEP2 的 文 件 名 可 由 用 戶 自 由 定 義 （ 限 制 為 4個 字符 ） ， 通 過 人 機(jī) 交 互 方 式 確 定 。 設(shè) 計 中 引

10、入 以 下 技 巧 ：1) 規(guī) 定 輸 入 文 件 名 FI與 輸 出 文 件 名 FO為 8個 字 符 2) 規(guī) 定 兩 個 字 符 數(shù) 組 NAMINP(8)， NAMOUT(8)3) 用 IQUIVALENCE語 句 使 FI與 NAMINP、 FO與NAMOUT 等 價4) 用 DATA語 句 給 FI和 FO賦 初 值 : “ .DAT”、 “ .OUT”5) 輸 入 NAMINP的 前 四 個 字 符 作 為 輸 入 輸 出 文 件 名 COMMON /PSIZE/MAXM,MAXA CHARACTER*8 FI, FO CHARACTER NAMINP(8), NAMOUT(8)

11、, HEAD*50 COMMON/HEAD/HEAD1 EQUIVALENCE (FI,NAMINP(1),(FO,NAMOUT(1) DATA FI,FO/ .DAT, .OUT/ WRITE(*, (A) INPUT FILE NAME (4 LETTERS ONLY) : READ (*, (4A1) (NAMINP(I), I=1,4) DO 10 I=1,4 10 NAMOUT(I)=NAMINP(I) OPEN(6, FILE=FO) OPEN(5, FILE=FI) MAXM=16000 MAXA=16380 CALL PCONTR(FO) END 二 、 定 義 數(shù) 組 FEP

12、2 中 設(shè) 置 了 兩 個 數(shù) 組 M(1600) 和 A(16380)， 數(shù)組 M是 作 為 存 放 坐 標(biāo) 、 單 元 信 息 、 載 荷 、 位 移 等 ， 數(shù)組 A 則 是 存 放 剛 度 矩 陣 用 的 ， 其 上 界 與 FORTRAN（ 編 譯 器 ） 版 本 和 計 算 機(jī) 內(nèi) 存 有 關(guān) ， 確 定 了 程 序 的 計算 規(guī) 模 。 這 兩 個 數(shù) 組 也 可 合 并 成 一 個 。 2-4 FEP2 的 主 控 程 序 PCONTR ： 實 質(zhì) 上 的 主 程 序 。 分 以 下 幾 大 段 ：一 、 程 序 頭 （ 前 12語 句 ） SUBROUTINE PCONTR(

13、FO) CHARACTER*8 FO LOGICAL ERR, ASEM COMMON /M/M(6000) COMMON /A/A(16380) COMMON /CDATA/NDF,NDM,NEN,NJ,NE,NEQ COMMON /PSIZE/MAXM,MAXA COMMON /MDATA/NUMMAT CHARACTER FD*25,FORC*6,P COMMON /PRINT/P COMMON /ASEM/ASEM DATA FORC/ FORC /FD/ NODAL FORCE/DISPL/ FO： 輸 出 文 件 名 ERR： 邏 輯 變 量 ， 作 出 錯 信 息 控 制 （ .

14、TRUE. 為 錯 ） ASEM ： 邏 輯 變 量 ， 作 為 判 斷 是 否 要 組 集 總 剛 數(shù) 組 A： 存 放 總 剛 的 元 素 數(shù) 組 M： 存 放 單 剛 、 單 元 信 息 、 節(jié) 點 坐 標(biāo) 等 NDF 最 大 自 由 度 數(shù) NDM 空 間 維 數(shù) NEN 單 元 的 最 多 節(jié) 點 數(shù) NJ 節(jié) 點 總 數(shù) NE 總 單 元 數(shù) NEQ 總 方 程 數(shù)NUMMAT 材 料 類 型 數(shù) FD , FORC ： 文 字 變 量 ， 作 為 輸 出 的 標(biāo) 題 用 P ： 文 字 變 量 ， 判 斷 是 否 要 輸 出 單 剛 二 、 輸 入 控 制 信 息 、 確 定 主

15、 要 數(shù) 組 的 地 址 （ 1434語 句 ） NEN1=NEN+1 每 個 單 元 信 息 的 存 儲 量 （ 節(jié) 點 號 + 材 料 號 ） NNEQ=NJ*NDF 總 節(jié) 點 節(jié) 點 自 由 度 數(shù) NST=NDF*NEN 單 剛 階 數(shù) NXC=1 節(jié) 點 坐 標(biāo) 數(shù) 組 XC的 首 址 NIX =NXC+NJ*NDM 單 元 信 息 數(shù) 組 IX的 首 址 NID =NIX+NEN1*NE 結(jié) 構(gòu) 方 程 號 編 碼 數(shù) 組 ID的 首 址 ND =NID+NNEQ 材 料 常 數(shù) 數(shù) 組 的 首 址 NXL =ND +20*NUMMAT 單 元 節(jié) 點 坐 標(biāo) 數(shù) 組 XL的 首

16、 址 NLD =NXL+NEN*NDM 單 元 節(jié) 點 自 由 度 數(shù) 組 LD的 首 址 NUL=NLD+NST 單 元 節(jié) 點 位 移 數(shù) 組 UL的 首 址 NS =NUL+NST*2 單 剛 數(shù) 組 S的 首 址 NP =NS +NST*NST 總 剛 對 角 線 地 址 向 量 JP的 首 址 READ(5,*) NJ,NE,NUMMAT,NDF,NDM,NEN WRITE(*,2000) NJ,NE,NUMMAT,NDF,NDM,NEN WRITE(6,2000) NJ,NE,NUMMAT,NDF,NDM,NEN NEN1=NEN+1 NNEQ=NJ*NDF NST=NDF*NEN

17、 NXC=1 NIX=NXC+NJ*NDM NS =NUL+NST*2 NP =NS +NST*NST NF =NP +NST NU =NF +NNEQ NJP=NU +NNEQ NEND=NJP+NNEQ 2-5 FEP2的 數(shù) 據(jù) 輸 入 格 式 與 程 序 實 現(xiàn)FEP2的 數(shù) 據(jù) 信 息 （ 自 由 格 式 ） ：1） 控 制 信 息 ： NJ， NE， NUMMAT， NDM， NEN2） 坐 標(biāo) 信 息 （ 一 組 ） N， NG， （ XL(I) ， I=1， NDM） 節(jié) 點 號 節(jié) 點 號 增 量 坐 標(biāo) 或 載 荷 由 GENVEC生 成 節(jié) 點 坐 標(biāo) 或 載 荷 向 量

18、 ， 通 過 NDM, CD， FC 進(jìn) 行 識 別 生 成 的 是 節(jié) 點 坐 標(biāo) ， 還 是 載 荷 向 量 。CD： NODAL COORDINATES FC ： COORCD： NODAL FORCE/DISPL FC ： FORC SUBROUTINE GENVEC(NDM, X, CD, FC, NJ, ERR) CHARACTER CD*18, FC*6 LOGICAL ERR REAL X(NDM,*),XL(6) ERR=.FALSE. N=0 NG=0102 L=N LG=NG READ (5,*,ERR=999) N,NG,(XL(I),I=1,NDM)101 IF (N

19、.LE.0 .OR. N.GT.NJ) GO TO 109 DO 103 I=1,NDM103 X(I,N) = XL(I) IF (LG) 104,102,104 104 LG=ISIGN(LG,N-L) LI= (IABS(N-L+LG)-1)/IABS(LG) DO 105 I=1,NDM 105 XL(I)=(X(I,N)-X(I,L)/LI106 L = L+LG IF(N-L)*LG.LE.0) GO TO 102 IF(L.LE.0.OR.L.GT.NJ) GO TO 108 DO 107 I=1,NDM LMLG=L-LG107 X(I,L)=X(I,LMLG)+XL(I) G

20、O TO 106108 WRITE(6,3000) L,CD WRITE(*,3000) L,CD ERR=.TRUE. GO TO 102 109 CONTINUE 3） 單 元 信 息 （ 一 組 ） L， LX， LK， （ IDL(K) ， K=1， NEN）單 元 號 單 元 號 增 量 材 料 號 節(jié) 點 號 0， 0， 0， 結(jié) 束由 ELDA子 程 序 生 成 單 元 信 息 ， IX(NEN1, NE) 二 維 數(shù) 組 。 NEN = NEN + 1 NE 單 元 總 數(shù) 4） 材 料 信 息 （ 一 組 或 多 組 ） 對 于 每 一 類 材 料 ， 要 求 輸 入 兩 個

21、 記 錄 ：a） 材 料 類 型 號 ， 單 元 類 型 號 b） 一 批 與 單 元 類 型 相 關(guān) 的 常 數(shù) （ 不 超 過 20個 ） 對 于 梁 單 元 ， 輸 入 兩 個 如 下 12個 常 數(shù) ： NP, E, G, A或 Ix , Iy或 Iz， Rh0， As， N1， N2， Ni， W, a類 型 （ NP=0： 平 面 梁 作 用 平 面 載 荷 ； 1： 平 面 梁 作 用 空 間 載 荷 ）Rh0 ： 單 位 長 度 的 質(zhì) 量 密 度N 1 ： 梁 左 端 鉸 為 1， 非 鉸 為 0 N2： 右 端 鉸 為 1， 非 鉸 為 0 Ni ： 單 元 載 荷 類 型

22、 （ 0， 1， 2， 3， 4） （ 0為 無 載 荷 ）單 元 載 荷 由 EQLOAD 處 理 成 等 效 節(jié) 點 載 荷 對 于 平 面 元 ， 輸 入 6個 常 數(shù) ： E, XNU, TH ， L， K， I TH： 單 元 厚 度 L： 沿 一 個 方 向 的 高 斯 積 分 點 數(shù) （ 2， 3， 4） K： 沿 一 個 方 向 應(yīng) 力 輸 出 的 高 斯 點 數(shù) I： =0 平 面 應(yīng) 變 0 平 面 應(yīng) 力 PMESH 調(diào) 用 ELEMLIB， 再 調(diào) 用 BEAM or PLANE SUBROUTINE PMESH(XC,IX,ID,D,NEN1) CHARACTER C

23、OOR*6,XX*18 COMMON /CDATA/NDF, NDM, NEN, NJ, NE, NEQ COMMON /MDATA/NUMMAT COMMON /ELDATA/LM, N, MA, MCT, IEL, NEL DIMENSION XC(NDM,*), IX(NEN1,*), ID(NDF,*), D(20,*) LOGICAL ERR DATA XX/ NODAL COORDINATES/COOR/ COOR / 1 CALL GENVEC(NDM,XC,XX,COOR,NJ,ERR) 1 CALL GENVEC(NDM,XC,XX,COOR,NJ,ERR) IF(ERR)

24、WRITE(*,3000) IF(ERR) STOP 2 CALL ELDAT(IX,NEN1) 3 DO 304 N=1,NUMMAT WRITE(6,2000) MA,IEL WRITE(*,2000) MA,IEL IF(MA.EQ.0) GO TO 4 CALL PZERO(D(1,MA),20) D(20,MA)=IEL CALL LEMLIB(D(1,MA),P,XC,IX,S,P,NDF,NDM,NST,1)304 CONTINUE 2000 FORMAT(/ MATERIAL TYPE ,I4, ELEMENT TYPE ,I4/)4 CALL BOUN(ID) END 5）

25、約 束 信 息 （ 包 括 指 定 位 移 ） （ 一 組 ） PMESH 調(diào) 用 BOUN 輸 入 ： N， NG， （ IDL(I) ， I=1， NDF） 節(jié) 點 號 節(jié) 點 號 增 量 坐 標(biāo) 或 載 荷 0， 0， 0， （ 平 面 元 4個 0， 梁 元 5個 0）IDL(I)： 約 束 信 息 ， 被 約 束 的 自 由 度 給 非 零 值 （ 1或 -1） ，其 它 為 0。1： 該 自 由 度 被 約 束 ， 但 不 繼 續(xù) 生 成-1： 該 自 由 度 被 約 束 ， 要 求 繼 續(xù) 生 成 6） 載 荷 信 息 載 荷 輸 入 仍 由 GENVEC 完 成 ， 維 數(shù) ：

26、 NDF 對 于 指 定 位 移 ， 生 成 方 式 與 節(jié) 點 載 荷 相 同 ， 5） 輸 入約 束 信 息 ， 6） 輸 入 其 值 。 第 三 章 總 剛 度 矩 陣 的 變 帶 寬 存 貯 與 求 解 總 剛 度 矩 陣 的 變 帶 寬 存 貯 與 對 角 線 地 址 LDLT 三 角 分 解 自 由 度 指 示 矩 陣 ID與 對 角 線 地 址 矩 陣 3-1 總 剛 度 矩 陣 的 變 帶 寬 存 貯 與 對 角 線 地 址總 剛 度 矩 陣 的 性 質(zhì) ： 對 稱 、 正 定 、 稀 疏只 存 貯 半 帶 寬 內(nèi) 的 元 素 （ 下 三 角 、 變 帶 寬 ） 稀 疏 對 稱

27、 矩 陣 A中 的 元 素 可 用 兩 個 一 維 數(shù) 組 B和 JP完 全 確 定 。 B的 上 界 為 A中 的 下 三 角 、 變 帶 寬 內(nèi) 元 素 總數(shù) ， 存 貯 A內(nèi) 的 元 素 。 JP的 上 界 為 A的 階 數(shù) ， 存 貯 A的 各行 對 角 線 元 素 的 序 號 ， 稱 為 對 角 線 元 素 地 址 數(shù) 組 。對 應(yīng) 關(guān) 系 ： A ( I , J ) = B ( JP (I) I + J )每 一 行 第 一 個 非 零 元 素 的 列 號 M i : M1=1, Mi= I JP(I) + JP(I-1) + 1 SUBROUTINE PROFIL(JP,ID,I

28、X,NEN1,NK) COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ DIMENSION JP(1),ID(NDF,1),IX(NEN1,*),IDL(18) NEQ=0 NAD=0 DO 50 N=1,NJ DO 40 I=1,NDF J=ID(I,N) IF(J) 30,20,3020 NEQ=NEQ+1 ID(I,N)=NEQ JP(NEQ)=0 GO TO 4030 NAD=NAD-1 ID(I,N)=NAD40 CONTINUE 50 CONTINUE C.COMPUTE ROW BANDWIDTH DO 500 N=1,NE DO 400 I=1,NEN II

29、=IX(I,N) IF(II.EQ.0) GO TO 400 DO 300 K=1,NDF KK=ID(K,II) IF(KK.LE.0) GO TO 300 DO 200 J=1,NEN JJ=IX(J,N) IF(JJ.EQ.0) GO TO 200 DO 100 L=1,NDF LL=ID(L,JJ) IF(LL.LE.0) GO TO 100 M=MAX0(KK,LL) JP(M)=MAX0(JP(M),IABS(KK-LL) 100 CONTINUE200 CONTINUE300 CONTINUE400 CONTINUE500 CONTINUE C.COMPUTE DIAGONAL

30、POINTERS FOR PROFIL NK=1 JP(1)=1 IF(NEQ.EQ.1) RETURN DO 600 N=2,NEQ600 JP(N)=JP(N)+JP(N-1)+1 NK=JP(NEQ) WRITE(*,2000)NK2000 FORMAT( THE MAXIMUM STORAGE OF GLOBAL STIFFNESS =,I6/) RETURN END SUBROUTINE ADDSTF(LD,JP,S,NST) CHARACTER Y COMMON /PRINT/Y COMMON /A/A(16380) COMMON /CDATA/NDF, NDM, NEN, NEJ

31、, NE, NEQ COMMON /ELDATA/LM, N, MA ,MCT, IEL, NEL DIMENSION LD(*), JP(*), S(NST,NST) IF(Y.EQ.Y.OR.Y.EQ.y) WRITE(*,2000) N,S2000 FORMAT( ELEMENT NUMBER ,I4, ELEMENT STIFFNESS/(6E12.5) DO 200 J=1,NST K=LD(J) IF(K.LE.0) GO TO 200 L=JP(K)-K DO 100 I=1,NST M=LD(I) IF(M.GT.K.OR.M.LE.0) GO TO 100 M=L+M A(M

32、)=A(M)+S(J,I)100 CONTINUE200 CONTINUE RETURN END 3-2 LDLT 三 角 分 解求 解 方 程 ： Ax=b A 對 稱 、 正 定 矩 陣 = Ux=y U 單 位 上 三 角 矩 陣對 A進(jìn) 行 三 角 分 解 = A=LLTL 下 三 角 矩 陣 ， 通 過 比 較 系 數(shù) 法 確 定 ， 需 要 開 方 運 算LDLT三 角 分 解 = A=LDLTL 下 三 角 矩 陣 ， D 對 角 矩 陣 1 iiii LD11 ii AL 11jk kkjkikijij LLLAL i=1, , n; j=2, , iLDLT(A, JP, N

33、)DIMENSION A(*), JP(*) Ax=b = LDLT x=b 令 LT x=y = LD y=b 前 消 公 式 ：11 by 11ij jjjijii LyLby i=2, , n回 代 公 式 ： nnnn Lyx kknik kkiii LxLyx )( 1 i=n-1, , 1 SUBROUTINE FORBACK(A, B, JP, N) DIMENSION A(*), B(*), JP(*) 3-3 自 由 度 指 示 矩 陣 ID與 對 角 線 地 址 矩 陣ID(NDF, NJ)： 非 零 值 （ 1或 -1） 為 被 約 束 的 自 由 度 ， 0為 活 化

34、自 由 度 。 在 PROFIL 中 ， 改 造 數(shù) 組 ID， 原 來 的 零 元 素 （ 活 化自 由 度 ） 用 正 數(shù) 計 ， 原 來 的 非 零 元 素 （ 非 活 化 自 由 度 ）用 負(fù) 數(shù) 計 ， 全 部 零 元 素 （ 活 化 ） 的 個 數(shù) 記 為 NEQ（ 總剛 矩 陣 的 階 數(shù) ） 。 計 算 剛 度 矩 陣 的 帶 寬 的 方 案 比 較 ：1） 以 各 節(jié) 點 自 由 度 為 研 究 對 象 ， 利 用 單 元 信 息 IX， 找出 與 該 點 聯(lián) 系 的 節(jié) 點 ， 然 后 比 較 節(jié) 點 自 由 度 之 最 大 者 。運 算 量 ： NJ*NDF*NE*NEN

35、*NDF2） 逐 一 研 究 各 單 元 ， 分 別 取 出 各 節(jié) 點 的 各 個 自 由 度 與該 單 元 中 其 它 節(jié) 點 的 自 由 度 比 較 ， 差 值 最 大 者 再 與 所 研究 的 自 由 度 之 帶 寬 值 （ 初 值 為 零 ） 比 較 ， 取 其 中 最 大 者代 替 原 有 的 帶 寬 值 。運 算 量 ： NE*NEN*NDF*NEN*NDF兩 者 運 算 量 比 ： NJ : NE 選 方 案 2 帶 寬 =自 由 度 之 最 大 差 值 + 1 暫 存 放 于 JP 中 若 已 知 對 角 線 地 址 數(shù) 組 JP ， 則 第 I 行 的 帶 寬 ： JP (

36、 1 ) = 1, JP ( I ) JP ( I 1 ) 若 JP 存 放 的 是 帶 寬 , 則 對 角 線 地 址 數(shù) 組 ： JP ( 1 ) = JP ( 1 ) = 1 JP ( 2 ) = JP ( 2 ) + JP ( 1 ) JP ( I +1) = JP ( I + 1 ) + JP ( I ) 總 剛 矩 陣 ( 按 活 化 自 由 度 存 貯 ) 全 部 存 貯 量 : NK = JP ( NEQ ) 第 四 章 單 剛 與 載 荷 的 組 集 ， 指 定 約 束 的處 理 ， 單 元 分 析 的 預(yù) 處 理 子 程 序 PFORM 的 功 能 子 程 序 MODIFY

37、 完 成 將 指 定 位 移 轉(zhuǎn) 化 成 等 價 內(nèi)力 將 單 剛 組 集 成 總 剛 子 程 序 ADDSTF 子 程 序 BASBLY 組 集 載 荷 與 節(jié) 點 反 力 向 量 PRTREAC 輸 出 節(jié) 點 反 力 PRTDIS 輸 出 節(jié) 點 位 移 ELEMLIB 單 元 庫 管 理 程 序 4-1 子 程 序 PFORM 的 功 能 子 程 序 PFORM 起 到 承 上 啟 下 的 作 用 ， 為 調(diào) 用 與 單元 運 算 有 關(guān) 的 各 種 子 程 序 做 必 要 的 前 后 處 理 。 功 能 參 數(shù) ISW 、 ASEM ：ISW=3， ASEM=.TRUE. 構(gòu) 造 單

38、 剛 ， 并 組 集ISW=3， ASEM=.FALSE. 構(gòu) 造 單 元 載 荷 ， 組 集 ， 將 指 定 位 移 轉(zhuǎn) 化 為 成 等 效 載 荷ISW=4 計 算 單 元 內(nèi) 力 或 應(yīng) 力 ， 并 輸 出ISW=6 構(gòu) 造 單 元 內(nèi) 力 ， 并 組 集 成 結(jié) 構(gòu) 反 力ISW=1 讀 入 并 構(gòu) 造 有 關(guān) 的 材 料 常 數(shù) SUBROUTINE PFORM(F, LD, UL, XL, S, P, U, X, IX, D, ID, JP, NST, NEN1,ISW) LOGICAL ASEM COMMON/ASEM/ASEM COMMON/CDATA/NDF,NDM,NEN,

39、NJ,NE,NEQ COMMON/ELDATA/LM,N,MA,MCT,IEL,NEL DIMENSION UL(NDF,*), U(*), F(NDF,*), S(NST,*), P(*) , JP(*) , + XL(NDM,*), X(NDM,*), D(20,*), IX(NEN1,*), ID(NDF,*), + LD(NDF,*) IEL=0 MCT=0 ！ 設(shè) 置 打 印 輸 出 DO 110 N=1, NE CALL PZERO(UL,NST*(NST+3) MA=IX(NEN1,N) ！ 材 料 號 DO 108 I=1,NEN II=IX(I,N) IF(II.NE.0)

40、GO TO 105 DO 102 J=1,NDF102 LD(J,I)=0 GO TO 108105 IID=II*NDF-NDF NEL=I ！ 確 定 單 元 的 實 際 節(jié) 點 數(shù) DO 106 J1=1,NDM106 XL(J1,I)=X(J1,II) ！ 確 定 單 元 的 節(jié) 點 坐 標(biāo)DO 107 J=1,NDF K=ID(J,II) ！ 確 定 自 由 度 標(biāo) 志 數(shù) 組 IF(K.GT.0) UL(J,I)=U(K) ！ 確 定 單 元 的 活 化 自 由 度 數(shù) INEN=I+NEN IF(K.LT.0) UL(J,I)=U(NEQ-K) ！ 確 定 單 元 的 約 束 位

41、 移 IF(K.LT.0) UL(J,INEN)=F(J,II) ！ 確 定 單 元 的 指 定 位 移 IF(ISW.EQ.6) K= IID+J107 LD(J,I)=K ！ 確 定 自 由 度 標(biāo) 志 數(shù) 組108 CONTINUE IE20=D(20,MA) IF(IE20.NE.IEL) MCT=0 ！ 設(shè) 置 打 印 輸 出 IEL=IE20 ！ 材 料 號CALL ELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST, ISW) ！ 確 定 聯(lián) 系 信 息 IX IF(ASEM.AND.ISW.EQ.3) CALL ADDSTF(LD,JP,

42、S,NST)130 IF(ISW.EQ.6) CALL BASBLY(F,P,LD,NST)120 CONTINUE IF(.NOT.ASEM.AND.ISW.EQ.3) CALL MODIFY(U,LD,S,UL(1,NEN+1),NST) IF(.NOT.ASEM.AND.ISW.EQ.3) CALL BASBLY(U,P,LD,NST)109 CONTINUE110 CONTINUE ENDUL： 單 元 位 移 S： 單 剛 P： 單 元 內(nèi) 力 4-2 子 程 序 MODIFY 完 成 將 指 定 位 移 轉(zhuǎn) 化 成等 價 內(nèi) 力按 活 化 自 由 度 和 被 約 束 自 由 度 ，

43、 將 總 剛 方 程 分 塊 ： 2222121 1212111 FuKuK FuKuKu1 ： 活 化 自 由 度 的 位 移 u2 ： 被 約 束 自 由 度 的 位 移 （ 包 括 指 定 位 移 ）F1 已 知 F2 ： 未 知 SUBROUTINE MODIFY(U,LD,S,DUL,NS) REAL U, DUL, S DIMENSION U(1), LD(1), S(NS, 1), DUL(1) DO 110 I=1, NS K=LD(I) IF(K.LE.0) GO TO 110 DO 100 J=1,NS100 U(K)=U(K) - S(I, J) * DUL(J)110

44、CONTINUE RETURN END 4-3 將 單 剛 組 集 成 總 剛 子 程 序 ADDSTFCALL ELEMLIB(D(1,MA), UL, XL, IX(1,N), S, P, NDF, NDM, NST,ISW) PFOR 35構(gòu) 造 了 單 剛 ， 存 放 于 S 中 ， ADDSTF 完 成 組 集PFORM 中 ， LD ( NDF, NEN )ADDSTF 中 , LD ( NST )一 維 數(shù) 組 與 二 維 滿 存 數(shù) 組 的 元 素 對 應(yīng) 關(guān) 系 ：IJ = JP( I ) I + J I J SUBROUTINE ADDSTF(LD, JP, S, NST)

45、 CHARACTER Y DIMENSION LD(*), JP(*), S(NST, NST) IF(Y.EQ.Y.OR.Y.EQ.y) WRITE(*,2000) N,S2000 FORMAT( ELEMENT NUMBER ,I4, ELEMENT STIFFNESS/(6E12.5) DO 200 J=1,NST K=LD(J) IF(K.LE.0) GO TO 200 L=JP(K)-K DO 100 I=1,NST M=LD(I) IF(M.GT.K.OR.M.LE.0) GO TO 100 M=L+M A(M)=A(M)+S(J,I)100 CONTINUE 200 CONTIN

46、UE END 4-4子 程 序 BASBLY 組 集 載 荷 與 節(jié) 點 反 力 向 量單 元 自 由 度 指 示 數(shù) 組 LDISW=3 組 集 載 荷 向 量 ， 按 節(jié) 點 活 化 自 由 度 次 序 排 列ISW=6 節(jié) 點 反 力 向 量 ， 按 節(jié) 點 排 列 次 序 130 IF(ISW.EQ.6) CALL BASBLY(F,P,LD,NST) PFOR 37IF(.NOT.ASEM.AND.ISW.EQ.3) CALL BASBLY(U,P,LD,NST) PFOR 40 SUBROUTINE BASBLY(B,P,LD,NS) BASB 1 DIMENSION B(*),P

47、(*),LD(*) BASB 2 DO 100 I=1,NS BASB 4 II=LD(I) BASB 5 IF(II.GT.0) B(II)=B(II)+P(I) BASB 6100 CONTINUE BASB 7 RETURN BASB 8 END BASB 9 4-5 PRTREAC 輸 出 節(jié) 點 反 力節(jié) 點 反 力 兩 種 計 算 方 式 ：1 ） 由 總 剛 方 程 求 得 2） 分 別 計 算 每 個 單 元 內(nèi) 力 （ 單 剛 乘 以 單 元 節(jié) 點 位 移 ） ，然 后 組 集 成 總 體 內(nèi) 力 向 量因 總 剛 分 解 后 ， 已 不 存 在 了 ， 故 選 2） SU

48、BROUTINE PRTREAC(R) COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ REAL R(NDF,*),RSUM(6) NNEQ=NDF*NJ DO 50 K=1,NDF RSUM(K)=0.0 50 CONTINUE DO 100 N=1,NJ,50 J=MIN0(NJ,N+49) WRITE(6,2000) (K,K=1,NDF) DO 100 I=N,J DO 75 K=1,NDF R(K,I)=-R(K,I) RSUM(K)=RSUM(K)+R(K,I) 75 CONTINUE DO 76 K=1,NDF76 IF(ABS(R(K,I).GT.1.0E

49、-3) GO TO 77 GO TO 10077 WRITE(6,2001) I,(R(K,I),K=1,NDF)100 CONTINUE WRITE(6,2002)(RSUM(K),K=1,NDF) RETURN2000 FORMAT(/5X,NODAL REACTIONS/2X,NODE,6(I8, DOF)2001 FORMAT(I6,6F12.4)2002 FORMAT(3X,SUM,6F12.4) END 4-6 PRTDIS 輸 出 節(jié) 點 位 移節(jié) 點 位 移 按 節(jié) 點 次 序 重 新 排 列注 意 ：1） 解 得 得 位 移 向 量 為 活 化 自 由 度 ， 編 號 小 于

50、 NEQ2） 約 束 (包 括 指 定 位 移 )得 值 存 放 在 二 維 數(shù) 組 F（ 載 荷 ）中 SUBROUTINE PRTDIS(ID,U,F) COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ DIMENSION ID(NDF,*),U(*),F(NDF,*),UL(6) DO 102 II=1,NJ,50 WRITE(6,2000)(K,K=1,NDF) JJ=MIN0(NJ,II+49) DO 102 N=II,JJ DO 100 I=1,NDF K=ID(I,N) IF(K.LT.0) U(NEQ-K)=F(I,N) IF(K.LT.0) UL(I)=U

51、(NEQ-K)100 IF(K.GT.0) UL(I)=U(K) WRITE(6,2001) N,(UL(I),I=1,NDF)102 CONTINUE2000 FORMAT(/ /1X,NODAL DISPLACEMENTS/2X,NODE,6(I8,DISP)2001 FORMAT(I6,6E12.4) END 4-7 ELEMLIB 單 元 庫 管 理 程 序FEP2 中 已 有 兩 種 單 元 ： BEAM 、 PLANE若 要 增 加 單 元 ， 可 修 改 第 七 句 。 例 如 ， 已 編 好 一 個 板 單元 ， 名 為 ： PLATE（ 形 參 ） ， 則 程 序 修 改 如

52、 下 ： GOTO (1, 2, 3) IEL 3 CALL PLATE（ D , U, X, IX, S, P, I, J, K, ISW) RETURN SUBROUTINE ELEMLIB(D,U,X,IX,S,P,I,J,K,ISW) REAL P(K),S(K,K),U(1) DIMENSION IX(1),D(1),X(1) COMMON/CDATA/NDF,NDM,NEN,NJ,NE,NEQ COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL IF(IEL.LE.0.OR.IEL.GT.2) GO TO 400 GO TO (1,2) IEL1 CALL BE

53、AM(D,U,X,IX,S,P,I,J,K,ISW) GO TO 1002 CALL PLANE(D,U,X,IX,S,P,I,J,K,ISW) GO TO 100100 RETURN400 WRITE(*,4000) IEL STOP4000 FORMAT(5X, *FATAL ERROR* ELEMENT CLASS NUMBER,I5, INPUT) END 第 五 章 平 面 單 元 設(shè) 計 平 面 四 節(jié) 點 單 元 的 基 本 公 式 數(shù) 值 積 分 形 函 數(shù) 及 其 導(dǎo) 數(shù) 四 節(jié) 點 單 剛 的 構(gòu) 造 3 9 節(jié) 點 平 面 單 元 平 面 單 元 PLANE 子 程 序

54、5-1平 面 四 節(jié) 點 單 元 的 基 本 公 式 在 FEP2 中 ， 提 供 了 3 9節(jié) 點 平 面 單 元 ， 其 基 本形 式 是 四 節(jié) 點 等 參 元 （ 四 邊 形 單 元 ） 。子 程 序 ： PLANE(D, UL, XL, IX, S, P, NDF, NDM, NST, ISW) PLANE 中 包 含 的 子 程 序 PGAUSS(L,LINT,SG,TG,WG)SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.) SHAP2(SS,TT,SHP,IX,NEL)PSTRES(SIG,SIG(4),SIG(5),SIG(6)

55、 5-2 數(shù) 值 積 分 SUBROUTINE PGAUSS(L,LINT,R,Z,W) COMMON /ELDATA/ LM,N,MA,MCT,IEL,NEL DIMENSION LR(9),LZ(9),LW(9),R(3),Z(3),W(3),G4(4),H4(4) DATA LR/-1,1,1,-1,0,1,0,-1,0/,LZ/-1,-1,1,1,-1,0,1,0,0/ DATA LW/4*25,4*40,64/ LINT = L*L GO TO (1,2,3,4),L1 R(1) = 0.0E0 Z(1) = 0.0E0 W(1) = 4.0E0 IF(NEL.EQ.3) THEN

56、Z(1) = -1.0E0/3.0E0 W(1) = 3.0E0 END IF RETURN 2 G = 1.0E0/SQRT(3.0E0) DO 21 I = 1,4 R(I) = G*LR(I) Z(I) = G*LZ(I)21 W(I) = 1.0E0 RETURN3 G = SQRT(0.6E0) H = 1.0E0/81.0E0 DO 31 I = 1,9 R(I) = G*LR(I) Z(I) = G*LZ(I)31 W(I) = H*LW(I) RETURN 4 G = SQRT(4.8E0) H = SQRT(30.0E0)/36.0E0 G4(1) = SQRT(3.0E0+

57、G)/7.0E0) G4(4) = - G4(1) G4(2) = SQRT(3.0E0-G)/7.0E0) G4(3) = -G4(2) H4(1) = 0.50E0 H； H4(2) = 0.50E0 + H H4(3) = 0.50E0 + H； H4(4) = 0.50E0 - H I = 0 DO 41 J = 1,4 DO 41 K = 1,4 I = I + 1 R(I) = G4(K) Z(I) = G4(J) W(I) = H4(J)*H4(K)41 CONTINUE END 5-3 形 函 數(shù) 及 其 導(dǎo) 數(shù) SUBROUTINE SHAPE(SS,TT,X,SHP,XSJ

58、,NDM,NEL,IX,FLG)DIMENSION SHP(3,*), X(NDM,1), S(4),T(4),XS(2,2), SX(2,2), IX(*) DATA S/-0.5E0,0.5E0,0.5E0,-0.5E0/,T/-0.5E0,-0.5E0,0.5E0,0.5E0/ DO 100 I = 1,4 SHP(3,I) = (0.5e0+S(I)*SS)*(0.5e0+T(I)*TT) SHP(1,I) = S(I)*(0.5e0+T(I)*TT)100 SHP(2,I) = T(I)*(0.5e0+S(I)*SS) IF(NEL.GE.4) GO TO 120 C. FORM T

59、RIANGLE BY ADDING THIRD AND FOURTH TOGETHER DO 110 I=1,3110 SHP(I,3) = SHP(I,3)+SHP(I,4) ！ 3 節(jié) 點 平 面 單 元 C. ADD QUADRATIC TERMS IF NECESSARY 120 IF(NEL.GT.4) CALL SHAP2(SS,TT,SHP,IX,NEL) C. CONSTRUCT JACOBIAN AND ITS INVERSE DO 130 I = 1,NDM DO 130 J = 1,2 XS(I,J) = 0.0 DO 130 K = 1,NEL130 XS(I,J) =

60、 XS(I,J) + X(I,K)*SHP(J,K) XSJ = XS(1,1)*XS(2,2)-XS(1,2)*XS(2,1) IF(FLG) RETURN SX(1,1) = XS(2,2)/XSJ SX(2,2) = XS(1,1)/XSJ SX(1,2) =-XS(1,2)/XSJ SX(2,1) =-XS(2,1)/XSJC. FORM GLOBAL DERIVATIVES DO 140 I = 1,NEL TP = SHP(1,I)*SX(1,1)+SHP(2,I)*SX(2,1) SHP(2,I) = SHP(1,I)*SX(1,2)+SHP(2,I)*SX(2,2) 140 S

61、HP(1,I) = TP END 5-4 四 節(jié) 點 單 剛 的 構(gòu) 造比 較 三 種 不 同 的 算 法 ， 選 運 算 量 最 少 的 一 種 。1） 按 矩 陣 直 接 運 算2） 預(yù) 先 算 出 D Bj ， 再 與 BiT相 乘 ， 展 開3） 按 矩 陣 運 算 全 部 展 開計 算 ： BiT D Bj 3 L = D(5) IF(L*L.NE.LINT) CALL PGAUSS(L,LINT,SG,TG,WG) DO 320 L = 1,LINT CALL SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX, .FALSE.) XSJ = XSJ*W

62、G(L)*D(7) J1 = 1 DO 320 J=1,NEL W11 = SHP(1,J)*XSJ W12 = SHP(2,J)*XSJ K1 = J1 DO 310 K = J,NEL S(J1 ,K1 ) = S(J1 ,K1 ) + W11*SHP(1,K) S(J1 ,K1+1) = S(J1 ,K1+1) + W11*SHP(2,K) S(J1+1,K1 ) = S(J1+1,K1 ) + W12*SHP(1,K) S(J1+1,K1+1) = S(J1+1,K1+1) + W12*SHP(2,K) 310 K1 = K1 + NDF320 J1 = J1 + NDF 5-5 39

63、 節(jié) 點 平 面 單 元 3 節(jié) 點 平 面 單 元 SHAPE 121359 節(jié) 點 平 面 單 元 SUBROUTINE SHAP2(S,T,SHP,IX,NEL) DIMENSION IX(*),SHP(3,*) S2 = (1.0E0-S*S)/2.0E0 T2 = (1.0E0-T*T)/2.0E0 DO 100 I=5,9 DO 100 J = 1,3 100 SHP(J,I) = 0.0 IF(IX(5).EQ.0) GO TO 101 SHP(1,5) = -S*(1.0E0-T) SHP(2,5) = -S2 SHP(3,5) = S2*(1.0E0-T)101 IF(NEL

64、.LT.6) GO TO 107 IF(IX(6).EQ.0) GO TO 102 SHP(1,6) = T2 SHP(2,6) = -T*(1.0E0+S) SHP(3,6) = T2*(1.0E0+S) 102 IF(NEL.LT.7) GO TO 107 SHAP 19 IF(IX(7).EQ.0) GO TO 103 SHAP 20 SHP(1,7) = -S*(1.0E0+T) SHAP 21 SHP(2,7) = S2 SHAP 22 SHP(3,7) = S2*(1.0E0+T) SHAP 23103 IF(NEL.LT.8) GO TO 107 SHAP 24 IF(IX(8)

65、.EQ.0) GO TO 104 SHAP 25 SHP(1,8) = -T2 SHAP 26 SHP(2,8) = -T*(1.0E0-S) SHAP 27 SHP(3,8) = T2*(1.0E0-S) SHAP 28C. INTERIOR NODE (LAGRANGIAN)104 IF(NEL.LT.9) GO TO 107 SHAP 30 IF(IX(9).EQ.0) GO TO 107 SHAP 31 SHP(1,9) = -4.0E0*S*T2 SHAP 32 SHP(2,9) = -4.0E0*T*S2 SHAP 33 SHP(3,9) = 4.0E0*S2*T2 C. CORR

66、ECT EDGE NODES FOR INTERIOR NODE DO 106 J= 1,3 SHAP 36 DO 105 I = 1,4 SHAP 37105 SHP(J,I) = SHP(J,I) - 0.25E0*SHP(J,9) SHAP 38 DO 106 I = 5,8 SHAP 39106 IF(IX(I).NE.0) SHP(J,I) = SHP(J,I) - .5E0*SHP(J,9) SHAP 40C. CORRECT CORNER NODES FOR PRESENSE OF MIDSIDE NODES107 K = 8 SHAP 42 DO 109 I = 1,4 SHAP 43 L = I + 4 SHAP 44 DO 108 J = 1,3 SHAP 45108 SHP(J,I) = SHP(J,I) - 0.5E0*(SHP(J,K)+SHP(J,L) SHAP 46109 K = L SHAP 47 END 5-6 平 面 單 元 PLANE 子 程 序 SUBROUTINE PLANE(D,UL,XL,IX,S,P,NDF,NDM,NST,ISW) CH