《有限元程序设计》PPT课件

有 限 元 程 序 设 计 教 学 程 序 ： FEP1 李 明 瑞 编 FEP2 李 明 瑞 编 ZFEP 第 一 章 绪 论 有 限 元 程 序 的 基 本 内 容 有 限 元 求 解 器 1.1 有 限 元 程 序 的 基 本 内 容有 限 元 法 的 解 题 步 骤结 构 离 散 化 为有 限 元 网 格 计 算 单 元 刚 度矩 阵引 入 约 束 条 件 求 应 力 等解 方 程 组 建 立 总 刚 度 矩阵 及 外 载 有 限 元 程 序 的 基 本 内 容 ：数 据 输 入阶 段 有 限 元 矩 阵 的 计算 、 组 集 和 求 解 数 据 输 出阶 段前 处 理 处 理 器 后 处 理 一 、 数 据 输 入 阶 段 前 处 理 读 入 和 生 成 数 据 ， 形 成 有 限 元 网 格 ， 为 有 限 元矩 阵 计 算 作 准 备 。 （ 数 据 或 图 形 输 入 ）1. 标 题 和 控 制 信 息2. 计 算 存 贮 分 配3. 节 点 坐 标 和 约 束 信 息4. 单 元 信 息5. 材 料 信 息6. 载 荷 信 息 二 、 有 限 元 矩 阵 的 计 算 、 组 集 和 求 解 处 理 器 计 算 插 值 函 数 矩 阵 N 、 几 何 矩 阵 B 、 Jacob矩 阵 J 在 高 斯 点 进 行 数 值 积 分 ， 求 得 单 刚 、 单 元 载 荷 组 集 成 总 刚 、 总 载 求 解 三 、 数 据 输 出 阶 段 后 处 理输 出 结 果 （ 场 变 量 和 场 变 量 导 数 等 ） ， 打 印 结 果 。场 变 量 包 括 ： 位 移 、 温 度 、 流 场 的 速 度 势 等 ；场 变 量 导 数 包 括 ： 应 力 、 应 变 、 热 流 、 流 场 速 度 势 等输 出 形 式 ： 数 据 输 出 和 图 形 输 出 1-2 有 限 元 求 解 器带 宽 法 ： 只 存 储 总 刚 矩 阵 的 半 带 宽 内 的 元 素 。等 带 宽 法 ： 每 行 的 带 宽 相 等 ， 常 采 用 二 维 数 组 存 储变 带 宽 法 ： 每 行 的 带 宽 不 等 ， 常 采 用 一 维 数 组 存 储一 、 带 宽 法变 带 宽 分 块 求 解 思 想 ： 一 个 变 量 的 消 元 只 影 响 刚 度 矩 阵 的 一 部 分 元素 ， 较 大 的 节 点 分 量 方 程 组 可 以 分 成 较 小 的 局 部系 数 矩 阵 按 节 点 逐 步 求 解 。 波 前 法 是 另 一 种 分 块 储 存 的 变 带 宽 法 ， 其 消 元 次序 是 按 单 元 编 号 进 行 ， 边 组 集 边 消 元 。 调 入 内 存 的 单 元 所 保 留 的 节 点 称 为 波 前 节 点 ， 所消 去 的 节 点 称 为 消 元 节 点 ， 消 元 节 点 是 与 以 后 调 入 的单 元 没 有 联 系 的 节 点 ， 即 该 节 点 的 方 程 已 组 集 完 。 波 前 法 的 回 代 按 消 元 节 点 的 反 序 进 行 。 对 内 存 的需 求 取 决 与 最 大 波 前 宽 。二 、 波 前 法 V-Fortran 1. Fortran语 言 的 基 本 格 式2. 变 量3. 基 本 语 句4. 子 例 子 程 序5. 函 数 子 程 序 6. 其 它 功 能 、 模 块 第 二 章 FEP2 程 序 的 总 体 设 计 与 输 入数 据 FEP2 的 设 计 任 务 FEP2 的 结 构 设 计 FEP2 的 主 程 序 FEP2 的 主 控 程 序 FEP2 的 数 据 输 入 格 式 与 程 序 实 现 2-1 FEP2 的 设 计 任 务1. FEP2的 简 要 题 目 说 明 FEP2是 一 个 具 有 通 用 性 的 教 学 程 序 ， 可 用 于 计 算一 般 的 线 性 静 力 学 问 题 。 已 设 计 了 平 面 梁 单 元 与 平 面 3-9节 点 元 两 种 单 元 ， 但 留 下 接 口 。2. 支 持 软 件 与 硬 件 FORTRAN77以 上 编 译 器 、 各 种 微 机3. FEP2的 功 能 3. FEP2的 功 能1） 输 入 文 件 名 由 用 户 自 由 定 义 ， 但 限 制 为 4个 字 符 ， 输入 文 件 扩 展 名 为 “ DAT”， 输 出 文 件 扩 展 名 为 “ OUT”。2） 节 点 坐 标 、 单 元 信 息 等 具 有 线 性 自 动 生 成 功 能 。3） 可 以 处 理 多 种 工 况 、 多 种 类 型 单 元 组 合 结 构 问 题 。4） 可 以 处 理 固 定 约 束 和 指 定 位 移 。5） 采 用 变 带 宽 存 储 、 三 角 分 解 法 求 解 刚 度 方 程 。6） 为 多 种 单 元 留 下 接 口 。 2-2 FEP2 的 结 构 设 计FEP2： 由 形 式 主 程 序 、 主 控 程 序 、 功 能 模 块 组 成 。主 程 序 的 主 要 功 能 ： 定 义 输 入 、 输 出 文 件 名 ， 调 用 主 控程 序 PCONTR主 控 程 序 PCONTR的 主 要 功 能 模 块 ：1） 内 存 分 配 2） 网 格 生 成3） 变 带 宽 存 储 设 置 4） 单 刚 的 形 成 与 组 集5） 刚 度 矩 阵 的 分 解 6） 形 成 节 点 载 荷7） 形 成 与 组 集 单 元 载 荷 8） 求 解 位 移9） 求 单 元 应 力 10） 求 结 构 反 力 FEP2 程 序 框 图FEP2 （ 主 程 序 ） PCONTR(主 控 程 序 ) SETMEM 检 查 内 存 PROFIL 形 成 变 带 宽 刚 度 矩 阵 地 址 PFORM (3) 形 成 单 刚 并 组 集 LDLT 总 刚 的 三 角 分 解 GENVEC 形 成 节 点 载 荷 PFOM (3) 构 造 并 组 集 单 元 载 荷 FORBACK 前 消 回 代 求 出 位 移 PFORM (4) 计 算 单 元 应 力 PFORM (6) 计 算 结 构 反 力 2-3 FEP2 的 主 程 序一 、 文 件 名 FEP2 的 文 件 名 可 由 用 户 自 由 定 义 （ 限 制 为 4个 字符 ） ， 通 过 人 机 交 互 方 式 确 定 。 设 计 中 引 入 以 下 技 巧 ：1) 规 定 输 入 文 件 名 FI与 输 出 文 件 名 FO为 8个 字 符 2) 规 定 两 个 字 符 数 组 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), 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 二 、 定 义 数 组 FEP2 中 设 置 了 两 个 数 组 M(1600) 和 A(16380)， 数组 M是 作 为 存 放 坐 标 、 单 元 信 息 、 载 荷 、 位 移 等 ， 数组 A 则 是 存 放 刚 度 矩 阵 用 的 ， 其 上 界 与 FORTRAN（ 编 译 器 ） 版 本 和 计 算 机 内 存 有 关 ， 确 定 了 程 序 的 计算 规 模 。 这 两 个 数 组 也 可 合 并 成 一 个 。 2-4 FEP2 的 主 控 程 序 PCONTR ： 实 质 上 的 主 程 序 。 分 以 下 几 大 段 ：一 、 程 序 头 （ 前 12语 句 ） SUBROUTINE PCONTR(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： 逻 辑 变 量 ， 作 出 错 信 息 控 制 （ .TRUE. 为 错 ） ASEM ： 逻 辑 变 量 ， 作 为 判 断 是 否 要 组 集 总 刚 数 组 A： 存 放 总 刚 的 元 素 数 组 M： 存 放 单 刚 、 单 元 信 息 、 节 点 坐 标 等 NDF 最 大 自 由 度 数 NDM 空 间 维 数 NEN 单 元 的 最 多 节 点 数 NJ 节 点 总 数 NE 总 单 元 数 NEQ 总 方 程 数NUMMAT 材 料 类 型 数 FD , FORC ： 文 字 变 量 ， 作 为 输 出 的 标 题 用 P ： 文 字 变 量 ， 判 断 是 否 要 输 出 单 刚 二 、 输 入 控 制 信 息 、 确 定 主 要 数 组 的 地 址 （ 1434语 句 ） NEN1=NEN+1 每 个 单 元 信 息 的 存 储 量 （ 节 点 号 + 材 料 号 ） NNEQ=NJ*NDF 总 节 点 节 点 自 由 度 数 NST=NDF*NEN 单 刚 阶 数 NXC=1 节 点 坐 标 数 组 XC的 首 址 NIX =NXC+NJ*NDM 单 元 信 息 数 组 IX的 首 址 NID =NIX+NEN1*NE 结 构 方 程 号 编 码 数 组 ID的 首 址 ND =NID+NNEQ 材 料 常 数 数 组 的 首 址 NXL =ND +20*NUMMAT 单 元 节 点 坐 标 数 组 XL的 首 址 NLD =NXL+NEN*NDM 单 元 节 点 自 由 度 数 组 LD的 首 址 NUL=NLD+NST 单 元 节 点 位 移 数 组 UL的 首 址 NS =NUL+NST*2 单 刚 数 组 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 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的 数 据 输 入 格 式 与 程 序 实 现FEP2的 数 据 信 息 （ 自 由 格 式 ） ：1） 控 制 信 息 ： NJ， NE， NUMMAT， NDM， NEN2） 坐 标 信 息 （ 一 组 ） N， NG， （ XL(I) ， I=1， NDM） 节 点 号 节 点 号 增 量 坐 标 或 载 荷 由 GENVEC生 成 节 点 坐 标 或 载 荷 向 量 ， 通 过 NDM, CD， FC 进 行 识 别 生 成 的 是 节 点 坐 标 ， 还 是 载 荷 向 量 。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.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) GO 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）单 元 号 单 元 号 增 量 材 料 号 节 点 号 0， 0， 0， 结 束由 ELDA子 程 序 生 成 单 元 信 息 ， IX(NEN1, NE) 二 维 数 组 。 NEN = NEN + 1 NE 单 元 总 数 4） 材 料 信 息 （ 一 组 或 多 组 ） 对 于 每 一 类 材 料 ， 要 求 输 入 两 个 记 录 ：a） 材 料 类 型 号 ， 单 元 类 型 号 b） 一 批 与 单 元 类 型 相 关 的 常 数 （ 不 超 过 20个 ） 对 于 梁 单 元 ， 输 入 两 个 如 下 12个 常 数 ： NP, E, G, A或 Ix , Iy或 Iz， Rh0， As， N1， N2， Ni， W, a类 型 （ NP=0： 平 面 梁 作 用 平 面 载 荷 ； 1： 平 面 梁 作 用 空 间 载 荷 ）Rh0 ： 单 位 长 度 的 质 量 密 度N 1 ： 梁 左 端 铰 为 1， 非 铰 为 0 N2： 右 端 铰 为 1， 非 铰 为 0 Ni ： 单 元 载 荷 类 型 （ 0， 1， 2， 3， 4） （ 0为 无 载 荷 ）单 元 载 荷 由 EQLOAD 处 理 成 等 效 节 点 载 荷 对 于 平 面 元 ， 输 入 6个 常 数 ： E, XNU, TH ， L， K， I TH： 单 元 厚 度 L： 沿 一 个 方 向 的 高 斯 积 分 点 数 （ 2， 3， 4） K： 沿 一 个 方 向 应 力 输 出 的 高 斯 点 数 I： =0 平 面 应 变 0 平 面 应 力 PMESH 调 用 ELEMLIB， 再 调 用 BEAM or PLANE SUBROUTINE PMESH(XC,IX,ID,D,NEN1) CHARACTER COOR*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) 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） 约 束 信 息 （ 包 括 指 定 位 移 ） （ 一 组 ） PMESH 调 用 BOUN 输 入 ： N， NG， （ IDL(I) ， I=1， NDF） 节 点 号 节 点 号 增 量 坐 标 或 载 荷 0， 0， 0， （ 平 面 元 4个 0， 梁 元 5个 0）IDL(I)： 约 束 信 息 ， 被 约 束 的 自 由 度 给 非 零 值 （ 1或 -1） ，其 它 为 0。1： 该 自 由 度 被 约 束 ， 但 不 继 续 生 成-1： 该 自 由 度 被 约 束 ， 要 求 继 续 生 成 6） 载 荷 信 息 载 荷 输 入 仍 由 GENVEC 完 成 ， 维 数 ： NDF 对 于 指 定 位 移 ， 生 成 方 式 与 节 点 载 荷 相 同 ， 5） 输 入约 束 信 息 ， 6） 输 入 其 值 。 第 三 章 总 刚 度 矩 阵 的 变 带 宽 存 贮 与 求 解 总 刚 度 矩 阵 的 变 带 宽 存 贮 与 对 角 线 地 址 LDLT 三 角 分 解 自 由 度 指 示 矩 阵 ID与 对 角 线 地 址 矩 阵 3-1 总 刚 度 矩 阵 的 变 带 宽 存 贮 与 对 角 线 地 址总 刚 度 矩 阵 的 性 质 ： 对 称 、 正 定 、 稀 疏只 存 贮 半 带 宽 内 的 元 素 （ 下 三 角 、 变 带 宽 ） 稀 疏 对 称 矩 阵 A中 的 元 素 可 用 两 个 一 维 数 组 B和 JP完 全 确 定 。 B的 上 界 为 A中 的 下 三 角 、 变 带 宽 内 元 素 总数 ， 存 贮 A内 的 元 素 。 JP的 上 界 为 A的 阶 数 ， 存 贮 A的 各行 对 角 线 元 素 的 序 号 ， 称 为 对 角 线 元 素 地 址 数 组 。对 应 关 系 ： A ( I , J ) = B ( JP (I) I + J )每 一 行 第 一 个 非 零 元 素 的 列 号 M i : M1=1, Mi= I JP(I) + JP(I-1) + 1 SUBROUTINE PROFIL(JP,ID,IX,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=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 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, 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)=A(M)+S(J,I)100 CONTINUE200 CONTINUE RETURN END 3-2 LDLT 三 角 分 解求 解 方 程 ： Ax=b A 对 称 、 正 定 矩 阵 = Ux=y U 单 位 上 三 角 矩 阵对 A进 行 三 角 分 解 = A=LLTL 下 三 角 矩 阵 ， 通 过 比 较 系 数 法 确 定 ， 需 要 开 方 运 算LDLT三 角 分 解 = A=LDLTL 下 三 角 矩 阵 ， D 对 角 矩 阵 1 iiii LD11 ii AL 11jk kkjkikijij LLLAL i=1, , n; j=2, , iLDLT(A, JP, N)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为 活 化 自 由 度 。 在 PROFIL 中 ， 改 造 数 组 ID， 原 来 的 零 元 素 （ 活 化自 由 度 ） 用 正 数 计 ， 原 来 的 非 零 元 素 （ 非 活 化 自 由 度 ）用 负 数 计 ， 全 部 零 元 素 （ 活 化 ） 的 个 数 记 为 NEQ（ 总刚 矩 阵 的 阶 数 ） 。 计 算 刚 度 矩 阵 的 带 宽 的 方 案 比 较 ：1） 以 各 节 点 自 由 度 为 研 究 对 象 ， 利 用 单 元 信 息 IX， 找出 与 该 点 联 系 的 节 点 ， 然 后 比 较 节 点 自 由 度 之 最 大 者 。运 算 量 ： NJ*NDF*NE*NEN*NDF2） 逐 一 研 究 各 单 元 ， 分 别 取 出 各 节 点 的 各 个 自 由 度 与该 单 元 中 其 它 节 点 的 自 由 度 比 较 ， 差 值 最 大 者 再 与 所 研究 的 自 由 度 之 带 宽 值 （ 初 值 为 零 ） 比 较 ， 取 其 中 最 大 者代 替 原 有 的 带 宽 值 。运 算 量 ： NE*NEN*NDF*NEN*NDF两 者 运 算 量 比 ： NJ : NE 选 方 案 2 带 宽 =自 由 度 之 最 大 差 值 + 1 暂 存 放 于 JP 中 若 已 知 对 角 线 地 址 数 组 JP ， 则 第 I 行 的 带 宽 ： JP ( 1 ) = 1, JP ( I ) JP ( I 1 ) 若 JP 存 放 的 是 带 宽 , 则 对 角 线 地 址 数 组 ： JP ( 1 ) = JP ( 1 ) = 1 JP ( 2 ) = JP ( 2 ) + JP ( 1 ) JP ( I +1) = JP ( I + 1 ) + JP ( I ) 总 刚 矩 阵 ( 按 活 化 自 由 度 存 贮 ) 全 部 存 贮 量 : NK = JP ( NEQ ) 第 四 章 单 刚 与 载 荷 的 组 集 ， 指 定 约 束 的处 理 ， 单 元 分 析 的 预 处 理 子 程 序 PFORM 的 功 能 子 程 序 MODIFY 完 成 将 指 定 位 移 转 化 成 等 价 内力 将 单 刚 组 集 成 总 刚 子 程 序 ADDSTF 子 程 序 BASBLY 组 集 载 荷 与 节 点 反 力 向 量 PRTREAC 输 出 节 点 反 力 PRTDIS 输 出 节 点 位 移 ELEMLIB 单 元 库 管 理 程 序 4-1 子 程 序 PFORM 的 功 能 子 程 序 PFORM 起 到 承 上 启 下 的 作 用 ， 为 调 用 与 单元 运 算 有 关 的 各 种 子 程 序 做 必 要 的 前 后 处 理 。 功 能 参 数 ISW 、 ASEM ：ISW=3， ASEM=.TRUE. 构 造 单 刚 ， 并 组 集ISW=3， ASEM=.FALSE. 构 造 单 元 载 荷 ， 组 集 ， 将 指 定 位 移 转 化 为 成 等 效 载 荷ISW=4 计 算 单 元 内 力 或 应 力 ， 并 输 出ISW=6 构 造 单 元 内 力 ， 并 组 集 成 结 构 反 力ISW=1 读 入 并 构 造 有 关 的 材 料 常 数 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,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 ！ 设 置 打 印 输 出 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) GO TO 105 DO 102 J=1,NDF102 LD(J,I)=0 GO TO 108105 IID=II*NDF-NDF NEL=I ！ 确 定 单 元 的 实 际 节 点 数 DO 106 J1=1,NDM106 XL(J1,I)=X(J1,II) ！ 确 定 单 元 的 节 点 坐 标DO 107 J=1,NDF K=ID(J,II) ！ 确 定 自 由 度 标 志 数 组 IF(K.GT.0) UL(J,I)=U(K) ！ 确 定 单 元 的 活 化 自 由 度 数 INEN=I+NEN IF(K.LT.0) UL(J,I)=U(NEQ-K) ！ 确 定 单 元 的 约 束 位 移 IF(K.LT.0) UL(J,INEN)=F(J,II) ！ 确 定 单 元 的 指 定 位 移 IF(ISW.EQ.6) K= IID+J107 LD(J,I)=K ！ 确 定 自 由 度 标 志 数 组108 CONTINUE IE20=D(20,MA) IF(IE20.NE.IEL) MCT=0 ！ 设 置 打 印 输 出 IEL=IE20 ！ 材 料 号CALL ELEMLIB(D(1,MA),UL,XL,IX(1,N),S,P,NDF,NDM,NST, ISW) ！ 确 定 联 系 信 息 IX IF(ASEM.AND.ISW.EQ.3) CALL ADDSTF(LD,JP,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： 单 元 内 力 4-2 子 程 序 MODIFY 完 成 将 指 定 位 移 转 化 成等 价 内 力按 活 化 自 由 度 和 被 约 束 自 由 度 ， 将 总 刚 方 程 分 块 ： 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 CONTINUE RETURN END 4-3 将 单 刚 组 集 成 总 刚 子 程 序 ADDSTFCALL ELEMLIB(D(1,MA), UL, XL, IX(1,N), S, P, NDF, NDM, NST,ISW) PFOR 35构 造 了 单 刚 ， 存 放 于 S 中 ， ADDSTF 完 成 组 集PFORM 中 ， LD ( NDF, NEN )ADDSTF 中 , LD ( NST )一 维 数 组 与 二 维 满 存 数 组 的 元 素 对 应 关 系 ：IJ = JP( I ) I + J I J SUBROUTINE ADDSTF(LD, JP, S, NST) 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 CONTINUE END 4-4子 程 序 BASBLY 组 集 载 荷 与 节 点 反 力 向 量单 元 自 由 度 指 示 数 组 LDISW=3 组 集 载 荷 向 量 ， 按 节 点 活 化 自 由 度 次 序 排 列ISW=6 节 点 反 力 向 量 ， 按 节 点 排 列 次 序 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(*),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 输 出 节 点 反 力节 点 反 力 两 种 计 算 方 式 ：1 ） 由 总 刚 方 程 求 得 2） 分 别 计 算 每 个 单 元 内 力 （ 单 刚 乘 以 单 元 节 点 位 移 ） ，然 后 组 集 成 总 体 内 力 向 量因 总 刚 分 解 后 ， 已 不 存 在 了 ， 故 选 2） SUBROUTINE 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-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 输 出 节 点 位 移节 点 位 移 按 节 点 次 序 重 新 排 列注 意 ：1） 解 得 得 位 移 向 量 为 活 化 自 由 度 ， 编 号 小 于 NEQ2） 约 束 (包 括 指 定 位 移 )得 值 存 放 在 二 维 数 组 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(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（ 形 参 ） ， 则 程 序 修 改 如 下 ： 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 BEAM(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 第 五 章 平 面 单 元 设 计 平 面 四 节 点 单 元 的 基 本 公 式 数 值 积 分 形 函 数 及 其 导 数 四 节 点 单 刚 的 构 造 3 9 节 点 平 面 单 元 平 面 单 元 PLANE 子 程 序 5-1平 面 四 节 点 单 元 的 基 本 公 式 在 FEP2 中 ， 提 供 了 3 9节 点 平 面 单 元 ， 其 基 本形 式 是 四 节 点 等 参 元 （ 四 边 形 单 元 ） 。子 程 序 ： 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) 5-2 数 值 积 分 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 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+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 形 函 数 及 其 导 数 SUBROUTINE SHAPE(SS,TT,X,SHP,XSJ,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 TRIANGLE BY ADDING THIRD AND FOURTH TOGETHER DO 110 I=1,3110 SHP(I,3) = SHP(I,3)+SHP(I,4) ！ 3 节 点 平 面 单 元 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) = 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 SHP(1,I) = TP END 5-4 四 节 点 单 刚 的 构 造比 较 三 种 不 同 的 算 法 ， 选 运 算 量 最 少 的 一 种 。1） 按 矩 阵 直 接 运 算2） 预 先 算 出 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*WG(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 节 点 平 面 单 元 3 节 点 平 面 单 元 SHAPE 121359 节 点 平 面 单 元 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.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).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. CORRECT 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