各向同性线性三维稳态传热有限元计算程序
传热问题有限元分析
【问题描述】本例对覆铜板模型进行稳态传热以及热应力分析，图I所示的是铜带以及基板的俯视图，铜带和基板之间由很薄的胶层连接，可以认为二者之间为刚性连接，这样的模型不包含胶层，只有长10mm的铜带（横截面2mm×0.1mm）和同样长10mm的基板（横截面2mm×0.2mm）。
材料性能参数如表1所示，有限元分析模型为实体——实体单元，单元大小0.05mm，边界条件为基板下表面温度为100℃，铜带上表面温度为20℃，通过二者进行传热。
图I 铜带与基板的俯视图表1 材料性能参数名称弹性模量泊松比各向同性导热系数基板 3.5GPa 0.4 300W/(m·℃)铜带110GPa 0.34 401W/(m·℃)【要求】在ANSYS Workbench软件平台上，对该铜板及基板模型进行传热分析以及热应力分析。
1.分析系统选择（1）运行ANSYS Workbench，进入工作界面，首先设置模型单位。
在菜单栏中找到Units下拉菜单，依次选择Units>Metric（kg,m,s,℃，A，N,V）命令。
（2）在左侧工具箱【Toolbox】下方“分析系统”【Analysis Systems】中双击“稳态热分析”【Steady-State Thermal】系统，此时在右侧的“项目流程”【Project Schematic】中会出现该分析系统共7个单元格。
相关界面如图1所示。
图1 Workbench中设置稳态热分析系统（3）拖动左侧工具箱中“分析系统”【Analysis Systems】中的“静力分析”【Static Structural】系统进到稳态热分析系统的【Solution】单元格中，为之后热应力分析做准备。
完成后的相关界面如图2所示。
图2 热应力分析流程图2.输入材料属性（1）在右侧窗口的分析系统A中双击工程材料【Engineering Data】单元格，进入工程数据窗口。
传热学中几种常用的软件及数值解法的介绍
传热学中几种常用软件及数值解法的介绍一、常用软件介绍：1、FLUENT软件简介FLUENT软件是美国FLUENT公司开发的通用CFD流场计算分析软件，囊括了Fluent Dynamic International、比利时Polyflow和Fluent Dynamic International（FDI）的全部技术力量（前者是公认的粘弹性和聚合物流动模拟方面占领先地位的公司，而后者是基于有限元方法CFD软件方面领先的公司）。
FLUENT是用于计算流体流动和传热问题的程序。
由于采用了多种求解方法和多重网格加速收敛技术，因而FLUENT能达到最佳的收敛速度和求解精度。
灵活的非结构化网格和基于解的自适应网格技术及成熟的物理模型，使FLUENT在转捩与湍流、传热与相变、化学反应与燃烧、多相流、旋转机械、动/变形网格、噪声、材料加工、燃料电池等方面有广泛应用。
采用的数值解法有限体积法（Finite Volume Method）程序的结构FLUENT程序软件包由以下几个部分组成：(1)GAMBIT——用于建立几何结构和网格的生成。
(2)FLUENT——用于进行流动模拟计算的求解器。
(3)prePDF——用于模拟PDF燃烧过程。
(4)TGrid——用于从现有的边界网格生成体网格。
(5)Filters(Translators)—转换其他程序生成的网格，用于FLUENT计算。
FLUENT程序可以求解的问题(1)可压缩与不可压缩流动问题。
(2)稳态和瞬态流动问题。
(3)无黏流，层流及湍流问题。
(4)牛顿流体及非牛顿流体。
(5)对流换热问题(包括自然对流和混合对流)。
(6)导热与对流换热耦合问题。
(7)辐射换热。
(8)惯性坐标系和非惯性坐标系下的流动问题模拟。
(9)用Lagrangian轨道模型模拟稀疏相(颗粒，水滴，气泡等)。
(10)一维风扇、热交换器性能计算。
(11)两相流问题。
(12)复杂表面形状下的自由面流动问题。
稳态热传导问题的数值模拟
稳态热传导问题的数值模拟热传导是热能从高温区向低温区传递的过程，在自然界和工程应用中有广泛的应用。
当材料或物体的长度，面积和体积足够大以至于其中的热量可以被视为连续分布时，稳态热传导方程可以用来描述热传导现象。
本文将讨论如何通过数值模拟来解决稳态热传导问题。
1. 稳态热传导方程首先，我们来看一下稳态热传导方程。
稳态热传导方程最常用的形式是二维热传导方程和三维热传导方程。
对于二维情况，可以表示为：$$ \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=0 $$对于三维情况，可以表示为：$$ \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partialy^2}+\frac{\partial^2 T}{\partial z^2}=0 $$其中，T表示温度。
2. 数值模拟方法由于稳态热传导方程在大多数情况下很难用解析方法求解，因此数值模拟方法成为了解决该问题的主要方法之一。
这里我们主要介绍两种数值模拟方法：有限差分法和有限元法。
2.1 有限差分法有限差分法是一种基于迭代计算的数值模拟方法，它将区域离散化为小的网格，并通过有限差分来逼近上述方程。
具体来说，它将偏微分方程近似为差分方程，然后用迭代方法来逼近和求解问题。
在应用有限差分法时，需要将连续的区域离散化为小的网格。
然后，用相邻两个网格点的温度差来逼近该点处的温度。
具体来说，对于二维情况，可以用以下公式来表示：$$ \frac{T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)-4T(i,j)}{h^2}=0 $$其中，h表示网格尺寸，i和j分别表示网格的横向和纵向坐标。
通过递归求解该方程，可以得到整个区域内的温度分布。
2.2 有限元法有限元法是一种更通用的数值模拟方法，可以用于解决各种类型的偏微分方程。
LS-DYNA使用指南中文版本
ANSYS/LS-DYNA将显式有限元程序LS-DYNA和ANSYS程序强大的前后处理结合起来。
用LS-DYNA的显式算法能快速求解瞬时大变形动力学、大变形和多重非线性准静态问题以及复杂的接触碰撞问题。
使用本程序，可以用ANSYS建立模型，用LS-DYNA做显式求解，然后用标准的ANSYS后处理来观看结果。
也可以在ANSYS和ANSYS-LS-DYNA之间传递几何信息和结果信息以执行连续的隐式－显式/显式－隐式分析，如坠落实验、回弹、及其它需要此类分析的应用。
显式动态分析求解步骤概述显式动态分析求解过程与ANSYS程序中其他分析过程类似，主要由三个步骤组成：1：建立模型（用PREP7前处理器）2：加载并求解（用SOLUTION处理器）3：查看结果（用POST1和POST26后处理器）本手册主要讲述了ANSYS/LS-DYNA显式动态分析过程的独特过程和概念。
没有详细论述上面的三个步骤。
如果熟悉ANSYS程序，已经知道怎样执行这些步骤，那么本手册将提供执行显式动态分析所需的其他信息。
如果从未用过ANSYS,就需通过以下两本手册了解基本的分析求解过程：·ANSYS Basic Analysis Guide·ANSYS Modeling and Meshing Guide使用ANSYS/LS-DYNA时，我们建议用户使用程序提供的缺省设置。
多数情况下，这些设置适合于所要求解的问题。
显式动态分析采用的命令在显式动态分析中，可以使用与其它ANSYS分析相同的命令来建立模型、执行求解。
同样，也可以采用ANSYS图形用户界面（GUI）中类似的选项来建模和求解。
然而，在显式动态分析中有一些独特的命令，如下：EDADAPT：激活自适应网格EDASMP：创建部件集合EDBOUND：定义一个滑移或循环对称界面EDBVIS：指定体积粘性系数EDBX：创建接触定义中使用的箱形体EDCADAPT：指定自适应网格控制EDCGEN：指定接触参数EDCLIST：列出接触实体定义EDCMORE：为给定的接触指定附加接触参数EDCNSTR：定义各种约束EDCONTACT：指定接触面控制EDCPU：指定CPU时间限制EDCRB：合并两个刚体EDCSC：定义是否使用子循环EDCTS：定义质量缩放因子EDCURVE：定义数据曲线EDDAMP：定义系统阻尼EDDC：删除或杀死/重激活接触实体定义EDDRELAX：进行有预载荷几何模型的初始化或显式分析的动力松弛EDDUMP：指定重启动文件的输出频率（d3dump）EDENERGY：定义能耗控制EDFPLOT：指定载荷标记绘图EDHGLS：定义沙漏系数EDHIST：定义时间历程输出EDHTIME：定义时间历程输出间隔EDINT：定义输出积分点的数目EDIS：定义完全重启动分析的应力初始化EDIPART：定义刚体惯性EDLCS：定义局部坐标系EDLOAD：定义载荷EDMP：定义材料特性EDNB：定义无反射边界EDNDTSD：清除噪声数据提供数据的图形化表示EDNROT：应用旋转坐标节点约束EDOPT：定义输出类型，ANSYS或LS-DYNA EDOUT：定义LS-DYNA ASCII输出文件EDPART：创建，更新，列出部件EDPC：选择、显示接触实体EDPL：绘制时间载荷曲线EDPVEL：在部件或部件集合上施加初始速度EDRC：指定刚体/变形体转换开关控制EDRD：刚体和变形体之间的相互转换EDREAD：把LS-DYNA的ASCII输出文件读入到POST26的变量中EDRI：为变形体转换成刚体时产生的刚体定义惯性特性EDRST：定义输出RST文件的时间间隔EDSHELL：定义壳单元的计算控制EDSOLV：把“显式动态分析”作为下一个状态主题EDSP：定义接触实体的小穿透检查EDSTART：定义分析状态（新分析或是重启动分析）EDTERM：定义中断标准EDTP：按照时间步长大小绘制单元EDVEL：给节点或节点组元施加初始速度EDWELD：定义无质量焊点或一般焊点EDWRITE：将显式动态输入写成LS-DYNA输入文件PARTSEL：选择部件集合RIMPORT：把一个显式分析得到的初始应力输入到ANSYSREXPORT：把一个隐式分析得到的位移输出到ANSYS/LS-DYNAUPGEOM：相加以前分析得到的位移，更新几何模型为变形构型关于ANSYS命令按字母顺序排列的详细资料(包括每条命令的特定路径)，请参阅《ANSYS Commands Reference》。
第7章 稳态热传导问题的有限元法
)dΒιβλιοθήκη 0(8-18)14
采度用分布Ga函ler数ki和n方换法热，边选界择条权件函代数入为(8，-w181 )式N，i 单将元单的元加内权的积温
分公式为
e
[ Ni x
(x
[N ]) Ni x y
( y
[N ])]{T}e d y
e
e
NiQ d 2 Ni qs d
(8-19)
e 3
Ni h[N ]{T}e d
一点上都满足边界条件（8-11）。对于复杂的工程问
题，这样的精确解往往很难找到，需要设法寻找近似
解。所选取的近似解是一族带有待定参数的已知函数
，一般表示为:
n
u u Ni ai Na
（8-12）
i 1
其中 ai为待定系数，为 Ni已知函数，称为试探函数。试探
函数要取完全的函数序列，是线性独立的。由于试探函数
T
0
t
5
这类问题称为稳态（Steady state）热传导问题。 稳态热传导问题并不是温度场不随时间变化，而是指 温度分布稳定后的状态。
若我们不关心物体内部的温度场如何从初始状态 过渡到最后的稳定温度场，那么随时间变化的瞬态（ Transient）热传导方程就退化为稳态热传导方程，三 维问题的稳态热传导方程为
，取： W j N j W j N j
下面用求解二阶常微分方程为例，说明Galerkin 法（参见，王勖成编著“有限元法基本原理和数值 方法”的1.2.3节）。
12
以二维问题为例，说明用Galerkin法建立稳态温度场 的一般有限元格式的过程。二维问题的稳态热传导方程：
x
x
T x
y
y
1 x j
三维有限元模型
三维有限元模型一、引言三维有限元模型是一种数学计算方法，用于分析和解决复杂的结构问题。
它可以将实际结构转化为由许多小单元组成的离散化模型，并通过数学方程求解每个单元的应力、应变等物理量，最终得出整个结构的响应。
本文将介绍三维有限元模型的基本原理、建模方法和求解过程。
二、三维有限元模型基本原理1. 有限元法基本思想有限元法是一种数值计算方法，它将一个连续的物理问题转化为由许多小单元组成的离散化问题，在每个小单元上建立数学模型，并通过求解代数方程组来得到整个系统的响应。
在三维有限元模型中，通常采用四面体或六面体等简单形状的单元进行离散化。
2. 三维有限元模型建立过程（1）几何建模：根据实际结构进行几何建模，包括确定结构尺寸、形状等。
（2）网格划分：将几何模型划分为许多小单元，并确定每个单元节点坐标。
（3）材料参数：根据实际材料性质确定每个单元的杨氏模量、泊松比等物理参数。
（4）载荷边界条件：根据实际工况确定结构所受载荷和边界条件。
（5）约束边界条件：根据实际结构确定约束边界条件，如支座、铰链等。
（6）求解：将以上信息输入计算机中，通过数学方法求解每个单元的应力、应变等物理量，并得出整个结构的响应。
三、三维有限元模型建模方法1. 网格划分方法三维有限元模型的网格划分可以采用手动或自动方式进行。
手动划分需要经验丰富的工程师进行，通常用于简单结构；自动划分则是利用计算机软件进行，可以快速生成复杂结构的网格。
2. 材料模型在三维有限元模型中，通常采用线性弹性模型来描述材料行为。
这种模型假设材料是各向同性的，并且满足胡克定律。
如果需要考虑非线性效应，则需要采用非线性材料模型。
3. 载荷和边界条件在三维有限元模型中，载荷和边界条件是建模的重要组成部分。
载荷可以是静载荷、动载荷或温度载荷等，边界条件可以是支座、铰链等。
四、三维有限元模型求解过程1. 单元刚度矩阵单元刚度矩阵是计算每个单元应力和应变的关键。
它由每个单元的杨氏模量、泊松比和几何信息确定。
三维有限元法计算过程
三维有限元法计算过程三维有限元法的计算过程：1）网格单元剖分；2）线性插值；3）单元分析；4）总体刚度矩阵合成；5）求解线性方程组等部分组成。
一、偏微分方程对应泛函的极值问题矿井稳恒电流场分布示意图主要任务是分析在给定边界条件下，求解稳定电流场的Laplace 方程或Poisson方程的数值解，即三维椭圆型微分方程的边值问题：)()((0)(0)()()(000z z y y x x I F u n un u F z u z y u y x u x Lu w D ---=⎪⎪⎪⎩⎪⎪⎪⎨⎧=+∂∂=∂∂=∂∂∂∂+∂∂∂∂+∂∂∂∂≡ΓΓ+Γδδδγσσσ 上述微分方程边值问题等价于下面泛函的极小值问题：dS U dxdydz fU z U y U x U U J w D ⎰⎰⎰⎰⎰Γ+Γ+ΓΩ+-∂∂+∂∂+∂∂=222221}])()()[(2{][γσσ二、网格剖分∞1ρi i h ρ............图2-4 单巷道地质模型1、网格单元的类型图2-5 网格单元类型2、网格单元剖分原则及其步长选择 因此，网格内的单元剖分应按以下剖分原则1）、各单元节点（顶点）只能与相邻单元节点（顶点）重合，而不能成为其它单元内点；2）、如果求解区域对称，那么单元剖分也应该对称；3）、在场变化剧烈的区域网格剖分单元要密一些，在场变化平缓的区域单元密度应小。
4）、网格单元体的大小变化应逐步过渡。
根据上述剖分原则，以x 、y 、z 坐标轴原点o 为中心，分别向x 、y 、z 方向的两侧作对称变步长剖分，距o 越远，步长应越大。
常用的变步长方法有：c i x x i i )1(1+=∆-∆+ c x x i i =∆∆+/1(i ≠0)c x x i i =∆-∆+111(i ≠0) 以上各式中c 为常数，1+∆i x 、i x ∆为同一坐标轴上相邻步长值。
以x 方向为例，可知，x 正方向与负方向对称，只相差一负号。
三维热传导方程的解法
三维热传导方程的解法热传导方程是热力学中的一个重要方程，用于描述物质内部温度随时间和位置的变化关系，常用来研究热传导现象和热工艺过程。
三维热传导方程是热传导方程的一种特殊形式，适用于描述三维体积内的热传导行为。
本文将介绍三维热传导方程的解法。
一、三维热传导方程的基本形式三维热传导方程的基本形式如下所示：$$\frac{\partial u}{\partial t} = \alpha \nabla^2 u$$其中，$u$ 表示温度场，$t$ 表示时间，$\alpha$ 为热扩散系数，$\nabla^2$ 是拉普拉斯算子，表示温度场的二阶空间导数之和。
二、三维热传导方程是一个偏微分方程，求解它的方法有很多种，以下将介绍其中的两种方法。
1. 分离变量法分离变量法是求解偏微分方程的常用方法之一，其基本思路是假设方程的解可以表示为若干个函数的乘积形式，然后通过代数推导得到这些函数的形式。
对于三维热传导方程，可以采用以下步骤进行求解：假设温度场 $u$ 可以表示为以下形式：$$u(x,y,z,t) = X(x)Y(y)Z(z)T(t)$$将上式代入三维热传导方程中，得到：$$\frac{X}{\alpha}\cdot\frac{d^2T}{dt^2} =\frac{T}{\alpha}\left(\frac{d^2X}{dx^2}+\frac{d^2Y}{dy^2}+\frac{d ^2Z}{dz^2}\right)$$假设方程的解为 $T(t)=e^{-\lambda\alpha t}$，其中$\lambda$ 为常数，则得到以下形式：$$\frac{X}{\alpha}\cdot\frac{d^2T}{dt^2} + \lambda T = 0$$通过求解上式可以得到 $T(t)$ 的形式。
进而，可以得到 $X(x)$、$Y(y)$ 和 $Z(z)$ 的形式。
将它们代入 $u$ 中，便可以得到温度场$u(x,y,z,t)$ 的解。
三维有限元非线性导热计算程序FemHC
三维有限元非线性导热计算程序FemHC 作者：刘瑜，朱可可，万秋里，王成恩来源：《计算机辅助工程》2022年第02期摘要：介紹采用有限元法求解稳态和瞬态非线性导热方程的计算程序FemHC，明确求解导热问题所依据的控制方程和有限元离散方法，描述FemHC的基本框架和核心类，并通过典型导热问题验证FemHC计算的准确性。
采用FemHC计算三维涡轮叶片的热传导，结果表明FemHC可以用于实际涡轮叶片等复杂结构导热问题的计算。
关键词：导热; 数值模拟; 非线性; 计算程序; 有限元中图分类号： TP391.99; TB115.1文献标志码： BThree dimensional finite element program FemHC fornonlinear heat conduction calculationLIU Yu1， ZHU Keke2， WAN Qiuli1，WANG Cheng’en1（1.School of Mechanical Engineering， Shanghai Jiao Tong University， Shanghai 200240，China;2.AECC Shenyang Engine Research Institute， Shenyang 110015， China）Abstract： The calculation program FemHC for solving steady state and transient nonlinear heat conduction equations by finite element method is introduced. The governing equation and finite element discretization method for solving the heat conduction problem are defined. The basic framework and core classes of FemHC is described. The calculation accuracy of FemHC is verified by typical heat conduction problems. The heat transfer of three-dimensional turbine blades is calculated by FemHC， and the results show that FemHC can be used to calculate the heat conduction problems of complex structures such as actual turbine blades.Key words： heat conduction; numerical simulation; nonlinear; calculation program; finite element0引言当在更严格和更复杂的条件下进行结构导热研究时，非线性导热问题分析变得更加重要，如涡轮叶片导热计算、功能梯度材料[1]导热计算等。
传热问题有限元分析
【问题描述】本例对覆铜板模型进行稳态传热以及热应力分析，图I所示的是铜带以及基板的俯视图，铜带和基板之间由很薄的胶层连接，可以认为二者之间为刚性连接，这样的模型不包含胶层，只有长10mm的铜带（横截面2mm×0.1mm）和同样长10mm的基板（横截面2mm×0.2mm）。
材料性能参数如表1所示，有限元分析模型为实体——实体单元，单元大小0.05mm，边界条件为基板下表面温度为100℃，铜带上表面温度为20℃，通过二者进行传热。
图I 铜带与基板的俯视图表1 材料性能参数名称弹性模量泊松比各向同性导热系数基板 3.5GPa 0.4 300W/(m·℃)铜带110GPa 0.34 401W/(m·℃)【要求】在ANSYS Workbench软件平台上，对该铜板及基板模型进行传热分析以及热应力分析。
1.分析系统选择（1）运行ANSYS Workbench，进入工作界面，首先设置模型单位。
在菜单栏中找到Units下拉菜单，依次选择Units>Metric（kg,m,s,℃，A，N,V）命令。
（2）在左侧工具箱【Toolbox】下方“分析系统”【Analysis Systems】中双击“稳态热分析”【Steady-State Thermal】系统，此时在右侧的“项目流程”【Project Schematic】中会出现该分析系统共7个单元格。
相关界面如图1所示。
图1 Workbench中设置稳态热分析系统（3）拖动左侧工具箱中“分析系统”【Analysis Systems】中的“静力分析”【Static Structural】系统进到稳态热分析系统的【Solution】单元格中，为之后热应力分析做准备。
完成后的相关界面如图2所示。
图2 热应力分析流程图。
四类问题有限元分析的MARC操作指南
2
1
图 1.4 方框拾取
★多边形拾取，如图 1.5 所示，将鼠标移至合适位置，按下键盘上的 Crtl 键，移动鼠标 在多边形每个角点处，按下<ML>，然后松开，直至形成多边形，拾取即告结束。
★闭环线拾取，如图 1.6 所示，将鼠标移至合适位置，按下键盘上的 Crtl 键，按下<ML>， 然后松开，移动鼠标直至闭环形成，拾取即告结束。
MSC.MARC/MARC 模块是功能齐全的高级非线性有限元软件求解器，体现了 30 年来 有限元分析的理论方法和软件实践的完美结合。它具有极强的结构分析能力：
(1) 可以处理各种线性和非线性结构分析包括：线性/非线性静力分析、模态分析、简谐 响应分析、频谱分析、随机振动分析、动力响应分析、自动的静/动力接触、屈曲/失稳、失 效和破坏分析等；
MENTAT 通过对话区的提问来要求用户输入数据或命令，根据提问结尾的标点符号不
同，共有 3 种不同的输入类型：
：
输入数据
>
输入字符串，通常是一个命令、文件名或集名
？
输入 YES 或 NO
用户可以通过鼠标、键盘或二者合用来输入数据、命令等、鼠标器用于菜单区进行选择
菜单以及在图形区检取点、线、面、单元、节点等。值得注意的是，鼠标键在菜单区和图形
1.3 应用 MSC.MARC 进行有限元分析的步骤
应用通用有限元软件进行分析，总体上分为前处理、求解和后处理三步。在 MSC.MARC 上，具体包括如下步骤：
(1) 生成网格（MESH GENERATION）；
6
四类问题有限元分析的MARC操作指南
清华大学机械工程系 曾攀
(2) 定义边界条件（BOUNDARY CONDITIONS）； (3) 定义初始条件（INITIAL CONDITIONS）； (4) 定义材料特性（MATERIAL PROPERTIES）； (5) 定义几何特性（GEOMETRIC PROPERTIES）； (6) 定义载荷工况（LOADCASES）； (7) 定义任务（JOBS）； (8) 后处理（RESULTS）。 下面通过一个简单的例子(图 1.7)，来说明运用 MSC.MARC 进行有限元分析的步骤。
各向同性线性三维动态传热有限元计算程序(有限元计算源码,有限元分析软件)
元计算有限元自动生成系统所开发源代码系列各向同性线性三维动态传热有限元计算程序1.简介元计算所开发的并行有限元程序自动生成系统（pFEPG）可根据用户需要开发出各种有限元计算程序源代码。
该源代码系列即为pFEPG所开发出来的求解各学科典型问题的有限元计算程序。
该组程序为各向同性线性三维动态传热有限元计算程序。
2.starta.for，对温度场的数据进行初始化；implicit real*8 (a-h,o-z)character*12 fname,filename(20)common /aa/ ia(250000000)common /bb/ ib(125000000)maxt=250000000/2c.... open disp0 file to get the numbers of nodes and degree of freedomc.... knode .... number of nodes, kdgof .... number ofd.o.f.open(1,file=' ',form='unformatted')read(1) knode,kdgofclose(1)kvar=knode*kdgofwrite(*,*) 'knode,kdgof,kvar ='write(*,'(1x,4i7)') knode,kdgof,kvarkvar1=kvar+1kcoor=3kelem=31250000knb1=kdgof*knode*1if (knb1/2*2 .lt. knb1) knb1=knb1+1kna4=kcoor*knode*2kna1=kdgof*knode*2kna2=kdgof*knode*2kna3=kdgof*knode*2kna5=knode*1if (kna5/2*2 .lt. kna5) kna5=kna5+1knb4=kelem*1if (knb4/2*2 .lt. knb4) knb4=knb4+1knb3=kvar1*1if (knb3/2*2 .lt. knb3) knb3=knb3+1knc1=kvar1*1if (knc1/2*2 .lt. knc1) knc1=knc1+1knb2=kvar1*1if (knb2/2*2 .lt. knb2) knb2=knb2+1knc3=maxt*1if (knc3/2*2 .lt. knc3) knc3=knc3+1knc2=kvar1*1if (knc2/2*2 .lt. knc2) knc2=knc2+1kna0=1kna1=kna1+kna0kna2=kna2+kna1kna3=kna3+kna2kna4=kna4+kna3kna5=kna5+kna4if (kna5-1.gt.250000000) thenwrite(*,*) 'exceed memory of array ia'write(*,*) 'memory of ia = 250000000'write(*,*) 'memory needed = ',kna5,' in prgram start' stop 55555endifknb0=1knb1=knb1+knb0knb2=knb2+knb1knb3=knb3+knb2knb4=knb4+knb3if (knb4-1.gt.125000000) thenwrite(*,*) 'exceed memory of array ib'write(*,*) 'memory of ib = 125000000'write(*,*) 'memory needed = ',knb4,' in prgram start' stop 55555endifknc0=1knc1=knc1+knc0knc2=knc2+knc1knc3=knc3+knc2if (knc3-1.gt.250000000) thenwrite(*,*) 'exceed memory of array ia'write(*,*) 'memory of ia = 250000000'write(*,*) 'memory needed = ',knc3,' in prgram start' stop 55555endifcall start(knode,kdgof,kcoor,kvar,*kelem,maxt,kvar1,ia(kna0),ia(kna1),ia(kna2),*ia(kna3),ia(kna4),ib(knb0),ib(knb1),ib(knb2),*ib(knb3),ia(knc0),ia(knc1),ia(knc2),*filename)endsubroutine start(knode,kdgof,kcoor,kvar,*kelem,maxt,kvar1,u0,u1,u2,*coor,inodvar,nodvar,numcol,lm,node,*jnz,jdiag,na,*filename)implicit real*8 (a-h,o-z)character*12 filename(20)DIMENSION NODV AR(KDGOF,KNODE),COOR(KCOOR,KNODE),R(3),* U0(KDGOF,KNODE),U1(KDGOF,KNODE),U2(KDGOF,KNODE),* INODV AR(KNODE),node(kelem)dimension lm(kvar1),jnz(kvar1),numcol(kvar1),na(maxt),& jdiag(kvar1)CHARACTER*1 MATERIALlogical filflgC .................................................................C ..... KDGOF NUMBER OF D.O.FC ..... KNODE NUMBER OF NODESC ..... INODV AR ID DATAC ..... NODV AR DENOTE THE EQUA TION NUMBER CORRESPONDING THE D.O.FC ..... U0 U1 U2 INITIAL V ALUEC ..... COOR COORDINATESC ..... NODE ELEMENT NODAL CONNECTIONC .................................................................6 FORMAT (1X, 15I4)7 FORMAT (1X,8F9.3)C.......OPEN ID fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) NUMNOD,NODDOF,((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE (1)call chms(kdgof,knode,NODV AR)c WRITE(*,*) 'NUMNOD =',NUMNOD,' NODDOF =',NODDOFc WRITE (*,*) 'ID ='c WRITE (*,6) ((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)C..... GET THE NATURAL NODAL ORDERDO 12 N=1,KNODEINODV AR(N)=N12 CONTINUEC..... OPEN ORDER.NOD FILE AND READ THE NODAL ORDER IF THE FILE EXISTinquire(file='ORDER.NOD',exist=filflg)if (filflg) thenOPEN (1,FILE='ORDER.NOD',FORM='UNFORMA TTED',STATUS='OLD')READ (1) (INODV AR(I),I=1,NUMNOD)CLOSE(1)WRITE(*,*) 'NODORDER ='WRITE(*,6) (INODV AR(I),I=1,NUMNOD)endifC..... GET NV BY IDNEQ=0DO 20 JNOD=1,NUMNODJ=INODV AR(JNOD)DO 18 I=1,NODDOFIF (NODV AR(I,J).NE.1) GOTO 18NEQ = NEQ + 1NODV AR(I,J) = NEQ18 CONTINUE20 CONTINUEDO 30 JNOD=1,NUMNODJ=INODV AR(JNOD)DO 28 I=1,NODDOFIF (NODVAR(I,J).GE.-1) GOTO 28N = -NODV AR(I,J)-1NODV AR(I,J) = NODV AR(I,N)28 CONTINUE30 CONTINUEC..... OPEN AND WRITE THE NV FILEOPEN(8,STATUS='unknown',FILE=' ' ,FORM='UNFORMA TTED')WRITE(8) ((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE(8)c WRITE(*,*) 'NUMNOD =',NUMNOD,' NODDOF =',NODDOFc WRITE(*,6) ((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)C.... WRITE THE BOUNDAY CONDITION FILE BFD ACCORDING TO THE DISP0 FILE C....OPEN DISP0 FILEOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ(1) NUMNOD,NODDOF,((U0(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE(1)C....OPEN BFD FILEOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')WRITE(1) ((U0(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE(1)C...... GET THE INITIAL TIME FROM TIME0 FILEC.......OPEN TIME0 FileOPEN(1,FILE=' ',FORM='FORMA TTED')READ(1,*) T0,TMAX,DTTIME = T0IT = 0WRITE(*,*) ' TMAX,DT,TIME,IT =',TMAX,DT,TIME,ITCLOSE(1)C.......OPEN TIME FileOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')WRITE(1) TMAX,DT,TIME,ITCLOSE(1)C.......OPEN COOR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)CLOSE(1)c WRITE(*,*) 'COOR ='c WRITE(*,7) ((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)C...... GET THE INITIAL V ALUE FROM THE DATA FILES BY PREPROCESSORinquire(file='disp1',exist=filflg)if (filflg) thenopen(16,file='disp1',form='unformatted',status='old')read(16) numnod,noddof,((U0(J,N),J=1,NODDOF),N=1,NUMNOD)close(16)endifinquire(file='disp2',exist=filflg)if (filflg) thenopen(16,file='disp2',form='unformatted',status='old')read(16) numnod,noddof,((U1(J,N),J=1,NODDOF),N=1,NUMNOD)close(16)endifinquire(file='disp3',exist=filflg)if (filflg) thenopen(16,file='disp3',form='unformatted',status='old')read(16) numnod,noddof,((U2(J,N),J=1,NODDOF),N=1,NUMNOD)close(16)endifc WRITE(*,*) ' U0 = 'c WRITE(*,'(6F13.3)') ((U0(J,N),J=1,NODDOF),N=1,NUMNOD)C WRITE(*,*) ' U1 = 'C WRITE(*,'(6F13.3)') ((U1(J,N),J=1,NODDOF),N=1,NUMNOD)C...... COMPUTE THE INITIAL V ALUE BY BOUND.FORzo = 0.0d0c DO 321 N=1,NUMNODc DO 100 J=1,NCOORc100 R(J) = COOR(J,N)c DO 200 J=1,NODDOFc U0(J,N) = BOUND(R,zo,J)c U1(J,N) = BOUND1(R,zo,J)c U2(J,N) = BOUND2(R,zo,J)c200 CONTINUEc321 CONTINUEC.......OPEN AND WRITE THE INITIAL V ALUE FILE UNODOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown') WRITE(1) ((U0(I,J),J=1,NUMNOD),I=1,NODDOF),* ((U1(I,J),J=1,NUMNOD),I=1,NODDOF),* ((U2(I,J),J=1,NUMNOD),I=1,NODDOF),* ((U0(I,J),J=1,NUMNOD),I=1,NODDOF)CLOSE (1)C.... Open IO fileopen(21,file=' ',form='formatted',status='old')read(21,'(1a)') materialread(21,*) numtypclose(21)DO I=1,NEQNUMCOL(I)=0jnz(I)=0ENDDOC.......OPEN ELEM0 fileOPEN (3,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')do jj=1,2NUMEL=0KELEM=0KEMATE=0DO 2000 ITYP=1,NUMTYPC.......INPUT ENODEREAD (3) NUM,NNODE,* ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODEcc WRITE(*,*) 'NODE ='cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) IF (KELEM.LT.NUM*NNODE) KELEM = NUM*NNODENNE = NNODEIF (MATERIAL.EQ.'Y' .OR. MATERIAL.EQ.'y') THENREAD (3) MMATE,NMATEIF (KEMATE.LT.MMATE*NMATE) KEMATE = MMATE*NMATENNE = NNE-1ENDIFWRITE(*,*) 'MMATE =',MMATE,' NMATE =',NMATEcc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODEcc WRITE(*,*) 'NODE ='cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) DO 1000 NE=1,NUML=0DO 700 INOD=1,NNENODI=NODE((NE-1)*NNODE+INOD)DO 600 IDGF=1,KDGOFINV=NODV AR(IDGF,NODI)IF (INV.LE.0) GOTO 600L=L+1LM(L)=INV600 CONTINUE700 CONTINUENUMEL=NUMEL+1C WRITE (*,*) 'L,LM =',LC WRITE (*,'(1X,15I5)') (LM(I),I=1,L)if (jj.eq.1) thenif (l.gt.0) call ACLn(L,jnz,lm)elseif (l.gt.0) call ACLH(NEQ,jnz,NUMCOL,NA,l,lm,MAXT)endif1000 continue2000 CONTINUEif (jj.eq.1) thenDO 800 N=1,NEQ-1800 jnz(N+1) = jnz(N+1)+jnz(N)DO 850 N=1,NEQ850 jnz(NEQ-N+2) = jnz(NEQ-N+1)jnz(1) = 0IF (JNZ(NEQ+1).GT.MAXT) THENWRITE(*,*) 'EXCEET CORE MEMORY ....',JNZ(NEQ+1),' > ',MAXTSTOP 1111ENDIFrewind(3)endifenddoCLOSE(3)call BCLH(NEQ,jnz,NUMCOL,NA,jdiag,LM,MAXT)MAXA=NUMCOL(NEQ+1)C.......OPEN SYS FileOPEN (2,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown') WRITE(2) NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATECLOSE (2)OPEN(2,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')write(2) (numcol(I),I=1,NEQ+1)write(2) (na(I),I=1,MAXA)write(2) (jdiag(i),i=1,neq)CLOSE(2)C write(*,*) 'NEQ,numcol,na=',NEQC write(*,6) (numcol(I),I=1,NEQ+1)C write(*,6) (na(i),i=1,MAXA)ENDsubroutine chms(kdgof,knode,id)dimension id(kdgof,knode),ms(1000),is(1000)do 1000 k=1,kdgofm = 0do 800 n=1,knodeif (id(k,n).le.-1) id(k,n)=-1if (id(k,n).le.1) goto 800j=id(k,n)j0=0if (m.gt.0) thendo i=1,mif (j.eq.ms(i)) j0=is(i)enddoendifif (j0.eq.0) thenm=m+1ms(m)=jis(m)=nid(k,n)=1elseid(k,n)=-j0-1endif800 continue1000 continuereturnendSUBROUTINE ACLn(ND,NUMCOL,LM)implicit real*8 (a-h,o-z)DIMENSION NUMCOL(*),LM(*)6 FORMAT(1X,10I7)C WRITE (*,*) 'ND= ',ND, (LM(I),I=1,ND)DO 400 I=1,NDNI = LM(I)DO 300 J=1,NDNJ = LM(J)cc IF (NJ.EQ.NI) GOTO 300NUMCOL(NI) = NUMCOL(NI)+1300 CONTINUE400 CONTINUEC WRITE(*,*) 'NUMCOL ='C WRITE(*,6) (NUMCOL(N),N=1,NEQ)1000 RETURNENDSUBROUTINE ACLH(NEQ,JNZ,NUMCOL,MHT,ND,LM,MAXT) implicit real*8 (a-h,o-z)DIMENSION MHT(MAXT),JNZ(*),NUMCOL(*),LM(*)6 FORMA T(1X,10I7)C WRITE (*,*) 'ND= ',ND, (LM(I),I=1,ND)DO 400 I=1,NDNI = LM(I)DO 300 J=1,NDC IF (J.EQ.I) GOTO 300NJ = LM(J)DO 200 K=1,NUMCOL(NI)IF (NJ.EQ.MHT(K+JNZ(NI))) GOTO 300200 CONTINUENUMCOL(NI) = NUMCOL(NI)+1IF (NI.GT.NEQ) THENWRITE(*,*) 'NI ====',NISTOP 1234ENDIFMHT(NUMCOL(NI)+JNZ(NI)) = NJ300 CONTINUE400 CONTINUEC WRITE(*,*) 'NUMCOL ='C WRITE(*,6) (NUMCOL(N),N=1,NEQ)C WRITE(*,*) 'MHT ='C DO 2000 N=1,NEQC2000 WRITE(*,6) (MHT(I+JNZ(N)),I=1,NUMCOL(N))1000 RETURNENDSUBROUTINE ORDER(ND,LM)implicit real*8 (a-h,o-z)DIMENSION LM(*)C WRITE(*,*) '**** ORDER ****'C WRITE(*,*) (LM(I),I=1,ND)DO 200 I=1,NDLS=LM(I)+1DO 100 J=I,NDIF (LM(J).GT.LS) GOTO 100LS=LM(J)J0=J100 CONTINUELM(J0)=LM(I)LM(I)=LS200 CONTINUEC WRITE(*,*) (LM(I),I=1,ND)C WRITE(*,*) '-----------------'RETURNENDSUBROUTINE BCLH(NEQ,JNZ,NUMCOL,NA,jdiag,LMI,MAXT) implicit real*8 (a-h,o-z)DIMENSION JNZ(*),NUMCOL(*),NA(*),jdiag(neq),LMI(*)DO 600 N=1,NEQCC NN = (N-1)*MAXCOLNN=JNZ(N)LI = NUMCOL(N)DO 500 I=1,LI500 LMI(I) = NA(NN+I)CALL ORDER(LI,LMI)DO 550 I=1,LI550 NA(NN+I) = LMI(I)600 CONTINUENSUM = 0DO 700 N=1,NEQCC NN=(N-1)*MAXCOLNN=JNZ(N)DO 700 I=1,NUMCOL(N)NSUM = NSUM+1NA(NSUM) = NA(NN+I)700 CONTINUEDO 800 N=1,NEQ-1800 NUMCOL(N+1) = NUMCOL(N+1)+NUMCOL(N)DO 850 N=1,NEQ850 NUMCOL(NEQ-N+2) = NUMCOL(NEQ-N+1)NUMCOL(1) = 0C WRITE(*,*) 'NSUM,NUMCOL(NEQ+1) =',NSUM,NUMCOL(NEQ+1)C WRITE(*,*) 'NUMCOL ='C WRITE(*,6) (NUMCOL(N),N=1,NEQ+1)C WRITE(*,*) 'NA ='C WRITE(*,6) (NA(I),I=1,NSUM)do 1400 ir=1, neq-1nir0 = numcol(ir)+1nir1 = numcol(ir+1)do 1380 id=nir0,nir1if(na(id).eq.ir) thenjdiag(ir)=idgoto 1381endif1380 continue1381 continue1400 continuejdiag(neq)=numcol(neq+1)1000 RETURN6 FORMA T(1X,5I15)END3.bft.for，更新时间并处理动态边界implicit real*8 (a-h,o-z)character*12 fname,filename(20)logical filflgcharacter*12 fname1,fname2common /aa/ ia(250000000)c.... open disp0 file to get the numbers of nodes and degree of freedom c.... knode .... number of nodes, kdgof .... number ofd.o.f.open(1,file=' ',form='unformatted')read(1) knode,kdgofclose(1)kvar=knode*kdgofKCOOR=3kna3=kdgof*knode*1if (kna3/2*2 .lt. kna3) kna3=kna3+1kna1=kcoor*knode*2kna2=kdgof*knode*2kna0=1kna1=kna1+kna0kna2=kna2+kna1kna3=kna3+kna2if (kna3-1.gt.250000000) thenwrite(*,*) 'exceed memory of array ia'write(*,*) 'memory of ia = 250000000'write(*,*) 'memory needed = ',kna3,' in prgram bft'stop 55555endifcall bft(knode,kdgof,kcoor,tmax,*time,dt,it,ia(kna0),ia(kna1),ia(kna2),*filename)C ...... OPEN THE FILE TO OBTAIN GRAPH FILE NAMESinquire(file='plotname',exist=filflg)if (.NOT. filflg) thenfname1 = 'unod'open(26,file='plotname',form='formatted',status='unknown')write(26,'(8a)') fname1close(26)endifopen(26,file='plotname',form='formatted',status='old')C ...... OPEN THE BATCH FILE FOR STORING THE RESULT TO GRAPHICopen(27,file='post.bat',form='formatted',status='unknown')C ...... STORE THE RESULT FOR EACH nstep TIME STEPcccc nstep = 1rstep = tmax/dt/100.0if (rstep.lt.1.5) thennstep = 1elsenstep = rstep+0.5endifik = it/nstepkk = it-ik*nstepcc write(*,*) 'nstep,it,ik,kk =',nstep,it,ik,kkif (kk.gt.0) goto 99999998 CONTINUEC ...... GET THE GRAPHIC FILE NAMEread(26,'(8a)',END=9999) fname1C WRITE(*,*) 'fname1 =',fname1fname2 = fname1call getname(fname2,ik)C ...... WRITE COPY COMMAND TO post.bat FILE FOR STORING THE RESULTwrite(27,*) 'copy ',fname1,' ',fname2GOTO 99989999 CONTINUEclose(26)write(27,*)close(27)c write(*,*) 'it =',itc write(*,*) fname1c write(*,*) fname2C.......OPEN STOP file IF THE LAST TIME IS ARRIVEDIF (TIME-TMAX.GT.-1.0d-20) THENcc OPEN(1,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown') open(16,file='plotno',form='formatted',status='unknown')itime = 1write(16,*) itime,ikclose(16)ENDIFendsubroutine bft(knode,kdgof,kcoor,tmax,*time,dt,it,coor,bf,nodvar,*filename)implicit real*8 (a-h,o-z)character*12 filename(20)DIMENSION NODV AR(KDGOF,KNODE),COOR(KCOOR,KNODE),* BF(KDGOF,KNODE),R(3)C.......OPEN TIME File AND UPDA TE THE TIMEOPEN(1,FILE=' ',FORM='UNFORMATTED')READ(1) TMAX,DT,TIME,ITT = TIME+DTTIME = TIME+DTIT = IT+1cc WRITE(*,*) ' TMAX,DT,TIME,IT =',TMAX,DT,TIME,IT REWIND(1)WRITE(1) TMAX,DT,TIME,ITCLOSE(1)C.......OPEN COOR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) KNODE,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,KNODE) CLOSE(1)C.......OPEN NODVAR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) ((NODV AR(I,J),I=1,KDGOF),J=1,KNODE)CLOSE (1)C WRITE(*,*) 'KNODE =',KNODE,' KDGOF =',KDGOFC WRITE (*,*) 'NODV AR ='C WRITE (*,6) ((NODV AR(I,J),I=1,KDGOF),J=1,KNODE)PUTE BOUNDARY CONDITIONDO 333 N=1,KNODEDO 100 J=1,NCOOR100 R(J) = COOR(J,N)DO 200 J=1,KDGOFID = NODV AR(J,N)BF(J,N) = 0.0IF (ID.LT.0) BF(J,N) = BOUND(R,T,J)C IF (ID.GT.0) BF(J,N) = FORCE(R,T,J)200 CONTINUE333 CONTINUECC WRITE(*,*) ' BF = 'CC WRITE(*,'(6F13.3)') ((BF(J,N),J=1,KDGOF),N=1,KNODE)C.......OPEN BFD file and WRITE BOUNDARY CONDITIONOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown')WRITE(1) ((BF(I,J),I=1,KDGOF),J=1,KNODE)CLOSE (1)C.......OPEN STOP file IF THE LAST TIME IS ARRIVEDIF (TIME-TMAX.GT.-1.0d-20) THENOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='unknown') ENDIFENDsubroutine getname(name,IT)implicit real*8 (a-h,o-z)character name*12,ch3*3c IF (IT.LT.10) WRITE(UNIT=CH3,FMT='(I1)') ITc IF (IT.GE.10) WRITE(UNIT=CH3,FMT='(I2)') ITc IF (IT.GE.100) WRITE(UNIT=CH3,FMT='(I3)') ITcall getext(it,ch3)c write(*,*) 'name =',name,'++++ CH3 =',CH3do 10 i=1,12if (name(i:i).eq.' ') thenj=igoto 20endif10 continue20 continueif (j.gt.9) thenwrite(*,*) 'Error, plot filename too long .......',namewrite(*,*) ' the length of name must be less or equal 8 character'stop 111endifc read(*,'(a3)') ch3name(j:j)='.'name(j+1:j+4)=ch3c write(*,*) 'plot filename = ',namereturnendsubroutine getext(ii,ch3)implicit real*8 (a-h,o-z)character ch3*3it = iich3 = ' 'k = 0if (ii.ge.100) thenn = it/100k = k+1call getchar(n,k,ch3)it = it - n*100endifif (ii.ge.10) thenn = it/10k = k+1call getchar(n,k,ch3)it = it - n*10endifn = itk = k+1call getchar(n,k,ch3)returnendsubroutine getchar(n,k,ch3)implicit real*8 (a-h,o-z)character ch3*3if (n.eq.0) ch3(k:k) = '0'if (n.eq.1) ch3(k:k) = '1'if (n.eq.2) ch3(k:k) = '2'if (n.eq.3) ch3(k:k) = '3'if (n.eq.4) ch3(k:k) = '4'if (n.eq.5) ch3(k:k) = '5'if (n.eq.6) ch3(k:k) = '6'if (n.eq.7) ch3(k:k) = '7'if (n.eq.8) ch3(k:k) = '8'if (n.eq.9) ch3(k:k) = '9'returnendreal*8 function bound(r,t,j)implicit real*8 (a-h,o-z)c implicit (a-h,o-z)dimension r(2)bound=t*t+r(1)*r(1)+r(2)*r(2) c write(*,*) 'bound =',boundreturnendreal*8 function bound1(r,t,j)implicit real*8 (a-h,o-z)c implicit (a-h,o-z)dimension r(2)bound1=2.*treturnendreal*8 function bound2(r,t,j)implicit real*8 (a-h,o-z)c implicit (a-h,o-z)dimension r(2)bound2=2.returnend4.eshunre3da.for，Galerkin法求解温度场的主程序implicit real*8 (a-h,o-z)character*12 fname,filename(20)common /aa/ ia(250000000)common /bb/ ib(125000000)common /cc/ ic(62500000)open(1,file=' ',form='unformatted',status='old')read(1) knode,kdgofclose(1)MAXT=250000000/(1+2*2)C.......OPEN SYS FileOPEN (2,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')read(2) NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATECLOSE (2)neq1=neq+1IF (MAXA.GT.MAXT) THENWRITE(*,*) 'MATRIX A EXCEED CORE MEMERY .... ',MAXAWRITE(*,*) 'REQUIRED CORE MEMERY ........... ',MAXTSTOP 0000ENDIFKV AR=KNODE*KDGOFKCOOR=3C KELEM=31250000WRITE(*,*) 'KNODE,KDGOF,KV AR,KCOOR,KELEM ='WRITE(*,'(1X,6I7)') KNODE,KDGOF,KVAR,KCOOR,KELEMkna3=kdgof*knode*1if (kna3/2*2 .lt. kna3) kna3=kna3+1knc1=kdgof*knode*2knc2=kcoor*knode*2knc11=kdgof*knode*2knc4=neq*2knb1=maxa*2knb3=maxa*1if (knb3/2*2 .lt. knb3) knb3=knb3+1kna1=neq1*1if (kna1/2*2 .lt. kna1) kna1=kna1+1kna2=neq*1if (kna2/2*2 .lt. kna2) kna2=kna2+1knc10=kemate*2kna4=kelem*1if (kna4/2*2 .lt. kna4) kna4=kna4+1knc12=100000*2knc3=neq*2knb2=maxa*2knc5=neq*2knc6=neq*2knc7=neq*2knc8=neq*2knc9=neq*2kna0=1kna1=kna1+kna0kna2=kna2+kna1kna3=kna3+kna2kna4=kna4+kna3if (kna4-1.gt.125000000) thenwrite(*,*) 'exceed memory of array ib'write(*,*) 'memory of ib = 125000000'write(*,*) 'memory needed = ',kna4,' in prgram eshunre3da' stop 55555endifknb0=1knb1=knb1+knb0knb2=knb2+knb1knb3=knb3+knb2if (knb3-1.gt.250000000) thenwrite(*,*) 'exceed memory of array ia'write(*,*) 'memory of ia = 250000000'write(*,*) 'memory needed = ',knb3,' in prgram eshunre3da' stop 55555endifknc0=1knc1=knc1+knc0knc2=knc2+knc1knc3=knc3+knc2knc4=knc4+knc3knc5=knc5+knc4knc6=knc6+knc5knc7=knc7+knc6knc8=knc8+knc7knc9=knc9+knc8knc10=knc10+knc9knc11=knc11+knc10knc12=knc12+knc11if (knc12-1.gt.62500000) thenwrite(*,*) 'exceed memory of array ic'write(*,*) 'memory of ic = 62500000'write(*,*) 'memory needed = ',knc12,' in prgram eshunre3da'stop 55555endifcall eshunre3da(knode,kdgof,kvar,kcoor,*numtyp,numel,neq,kelem,kemate,maxa,*maxt,neq1,ib(kna0),ib(kna1),ib(kna2),*ib(kna3),ia(knb0),ia(knb1),ia(knb2),ic(knc0),*ic(knc1),ic(knc2),ic(knc3),ic(knc4),ic(knc5),*ic(knc6),ic(knc7),ic(knc8),ic(knc9),ic(knc10),*ic(knc11),*filename)endsubroutine eshunre3da(knode,kdgof,kvar,kcoor,*numtyp,numel,neq,kelem,kemate,maxa,*maxt,neq1,numcol,jdiag,nodvar,node,*a,an,na,u,coor,u1,*f,x0,r0,p,w,hp,*emate,eu,sml,*filename)implicit real*8 (a-h,o-z)character*12 filename(20)DIMENSION NODV AR(KDGOF,KNODE),U(KDGOF,KNODE),COOR(KCOOR,KNODE), *eu(kdgof,knode),& F(NEQ),A(MAXA),na(maxa),numcol(neq1),jdiag(neq),EMA TE(KEMA TE),& NODE(KELEM),SML(100000),u1(neq),an(maxa),& X0(NEQ),R0(NEQ),p(NEQ),w(NEQ),hp(neq)6 FORMA T (1X,15I5)7 FORMA T (1X,5e15.5)1001 FORMA T(1X,9I7)C.......OPEN TIME FileOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ(1) TMAX,DT,TIME,ITWRITE(*,*) ' TMAX,DT,TIME,IT =',TMAX,DT,TIME,ITCLOSE(1)C.......OPEN NODVAR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) ((NODV AR(I,J),I=1,KDGOF),J=1,KNODE)CLOSE (1)cc WRITE(*,*) 'KDGOF =',KDGOF,' KNODE =',KNODEcc WRITE (*,*) 'NODV AR ='cc WRITE (*,6) ((NODV AR(I,J),I=1,KDGOF),J=1,KNODE)C.......OPEN COOR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD) CLOSE(1)cc WRITE(*,*) 'NUMNOD,NCOOR=',NUMNOD,NCOORC.......OPEN BF fileOPEN (32,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')cc READ (32) ((U(J,I),J=1,KDGOF),I=1,KNODE)cc CLOSE (32)cc WRITE (*,*) 'BF ='cc WRITE(*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)numtyp = 2open(11,file='unod',form='unformatted',status='old')read(11) ((eu(j,i),i=1,knode),j=1,kdgof)C.......OPEN DIAG fileOPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ(2) (numcol(I),I=1,NEQ1)READ(2) (na(I),I=1,maxa)READ(2) (jdiag(I),I=1,NEQ)CLOSE(2)C.......OPEN ELEM0 fileOPEN (3,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')c write(*,*) 'Input errilu'c read(*,*) rerrilurerrilu = 0.001itime=01 continueitime=itime+1errilu = rerriluwrite(*,*) 'errilu =',erriluif (itime.gt.1) thenwrite(*,*) 'Nonlinear Iteration Times ========',itimerewind(3)rewind(32)endifREAD (32) ((U(J,I),J=1,KDGOF),I=1,KNODE)cc WRITE (*,*) 'BF ='cc WRITE(*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)DO 111 I=1,MAXAA(I) = 0.0111 CONTINUEDO 2300 I=1,NEQ2300 CONTINUENUMEL=0C.......OPEN EMATE+ENODE+ELOAD fileDO 2000 ITYP=1,NUMTYPC.......INPUT ENODEREAD (3) NUM,NNODE,* ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODEcc WRITE(*,*) 'NODE ='cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) NNE = NNODEnne = nne-1K=0DO 115 J=1,NNEJNOD = NODE(J)DO 115 L=1,KDGOFIF (NODV AR(L,JNOD).NE.0) K=K+1115 CONTINUEWRITE(*,*) 'K =',Kkk=k*kk0=1k1=k0+k*kk2=k1+kk3=k2+kk4=k3+k*kk5=k4+k*kCALL ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE,* ITYP,NCOOR,NUM,TIME,DT,neq,neq1,maxt,MAXA,NODV AR,COOR,NODE,EMATE, & A,na,numcol,jdiag,&sml(k0),sml(k1),sml(k2),sml(k3),sml(k4),&eu,*U)2000 CONTINUECLOSE(1)CLOSE(2)close(11)DO 2050 IJ=1,NEQif (itime.le.1) u1(IJ) = 0.0F(IJ)=0.0D02050 CONTINUEDO 2200 I=1,KNODEDO 2100 J=1,KDGOFIJ=NODV AR(J,I)IF (IJ.LE.0) GOTO 2100F(IJ)=F(IJ)+U(J,I)2100 CONTINUE2200 CONTINUEcc WRITE (*,*) 'U ='cc WRITE (*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)cc WRITE (*,*) 'NEQ =',NEQ,' F ='cc WRITE(*,7) (F(I),I=1,NEQ)if (itime.le.1) thenC.......OPEN LMATRIX FILEOPEN (2,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown')CLOSE (2)endifMAXA = JDIAG(NEQ)C WRITE(*,*) 'JDIAG ='WRITE(*,*) 'ILUCG_SOLVER MEMORY REQUIRED .... ',MAXAIF (MAXA.GT.MAXT) THENWRITE(*,*) 'WARNING MA TRIX A EXCEED CORE MEMORY .... ',MAXTc STOP 0000ENDIFcall ilu(na, a, an, numcol, neq, maxa,jdiag)call redu(na, a,an, u1, f, numcol, neq, maxa, kkk,& x0,r0,p,w,hp,jdiag,errilu)C WRITE(*,*) ' U1 = 'C WRITE(*,7) (A(I),I,MAXA)C WRITE(*,7) (F(I),I=1,NEQ)NOUT = 20OPEN(NOUT,FILE=' ',FORM='FORMATTED',STA TUS='unknown')DO 3200 INOD=1,KNODEDO 3100 IDFG=1,KDGOFN=NODV AR(IDFG,INOD)C WRITE (*,*) 'N =',Nif(n.le.0) theneu(IDFG,INOD)=u(IDFG,INOD)elseeu(IDFG,INOD)=u1(N)endif3100 CONTINUE3200 CONTINUEDO 3400 N=1,KNODEWRITE (NOUT,3600) N,(eu(I,N),I=1,KDGOF)3400 CONTINUE3600 FORMA T (1X,I5,1X,6E11.4,9(/6X,6E11.4))CLOSE (NOUT)open(10,file='unod',form='unformatted',status='unknown')write(10) ((eu(j,i),i=1,knode),j=1,kdgof)close(10)close (3)close (32)RETURNENDSUBROUTINE ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE,& ITYP,NCOOR,NUM,TIME,DT,neq,neq1,maxt,MAXA,NODV AR,COOR,NODE,EMATE,& A,na,numcol,jdiag,*es,em,ef,Estifn,Estifv,eu,*U)implicit real*8 (a-h,o-z)DIMENSION NODV AR(KDGOF,KNODE),COOR(KCOOR,KNODE),NODE(KELEM), & U(KDGOF,KNODE),EMATE(300),& A(MAXa),na(MAXa),numcol(neq1),jdiag(neq),*es(k,k),em(k),ef(k),eu(kdgof,knode),*Estifn(k,k),Estifv(kk),*R(500),PRMT(500),COEF(500),LM(500)17 FORMA T (1X,15I5)18 FORMA T (1X,8e9.2)READ (3) MMATE,NMATE,((EMATE((I-1)*NMATE+J),J=1,NMATE),* I=1,MMATE)WRITE(*,*) 'MMATE =',MMATE,' NMATE =',NMATEWRITE (*,*) 'EMATE ='WRITE (*,18) ((EMATE((I-1)*NMATE+J),J=1,NMATE),* I=1,MMATE)DO 1000 NE=1,NUMNR=0DO 130 J=1,NNEJNOD = NODE((NE-1)*NNODE+J)IF (JNOD.LT.0) JNOD = -JNODPRMT(NMA TE+7+J) = JNODDO 120 I=1,NCOORNR=NR+1120 R(NR) = COOR(I,JNOD)130 CONTINUEIMA TE = NODE(NNODE*NE)DO 140 J=1,NMATE140 PRMT(J) = EMA TE((IMATE-1)*NMATE+J)PRMT(NMATE+1)=TIMEPRMT(NMATE+2)=DTPRMT(NMA TE+3)=IMATEprmt(NMA TE+4)=NEprmt(NMA TE+5)=NUMprmt(NMA TE+6)=ITprmt(NMA TE+7)=NMATEprmt(NMA TE+8)=ITIMEprmt(NMA TE+9)=ITYPgoto (1,2), ityp。
6 稳态热传导问题的有限元法
6. 穩態熱傳導問題的有限元法本章的內容如下：6.1熱傳導方程與換熱邊界6.2穩態溫度場分析的一般有限元列式 6.3三角形單元的有限元列式 6.4溫度場分析舉例6.1熱傳導方程與換熱邊界在分析工程問題時，經常要瞭解工件內部的溫度分佈情況，例如發動機的工作溫度、金屬工件在熱處理過程中的溫度變化、流體溫度分佈等。
物體內部的溫度分佈取決於物體內部的熱量交換，以及物體與外部介質之間的熱量交換，一般認為是與時間相關的。
物體內部的熱交換採用以下的熱傳導方程（Fourier 方程）來描述，Q z T z y T y x T x tT c+⎪⎭⎫⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂z y x λλλρ （6-1）式中ρ為密度，kg/m 3； c 為比熱容，K)J/(kg ⋅；z y x λλλ,,為導熱係數，)k m w ⋅；T 為溫度，℃；t 為時間，s ；Q 為內熱源密度，w/m 3。
對於各向同性材料，不同方向上的導熱係數相同，熱傳導方程可寫為以下形式，Q zT yT xT tT c222222+∂∂+∂∂+∂∂=∂∂λλλρ （6-2）除了熱傳導方程，計算物體內部的溫度分佈，還需要指定初始條件和邊界條件。
初始條件是指物體最初的溫度分佈情況，() z y,x,T T00t ==（6-3）邊界條件是指物體外表面與周圍環境的熱交換情況。
在傳熱學中一般把邊界條件分為三類。
1）給定物體邊界上的溫度，稱為第一類邊界條件。
物體表面上的溫度或溫度函數為已知，s sT T=或 ),,,(t z y x T Ts s=（6-4）2）給定物體邊界上的熱量輸入或輸出，稱為第二類邊界條件。
已知物體表面上熱流密度，s sz zy yx xq n zT n yT n xT =∂∂+∂∂+∂∂)(λλλ或),,,()(t z y x q n zT n yT n xT s sz zy yx x=∂∂+∂∂+∂∂λλλ（6-5）3）給定對流換熱條件，稱為第三類邊界條件。
稳态热传导问题有限元法
6. 稳态热传导问题的有限元法本章的内容如下：6.1热传导方程与换热边界6.2稳态温度场分析的一般有限元列式 6.3三角形单元的有限元列式 6.4温度场分析举例6.1热传导方程与换热边界在分析工程问题时，经常要了解工件内部的温度分布情况，例如发动机的工作温度、金属工件在热处理过程中的温度变化、流体温度分布等。
物体内部的温度分布取决于物体内部的热量交换，以及物体与外部介质之间的热量交换，一般认为是与时间相关的。
物体内部的热交换采用以下的热传导方程（Fourier 方程）来描述，Q z T z y T y x T x t T c+⎪⎭⎫⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂z y x λλλρ （6-1）式中ρ为密度，kg/m 3； c 为比热容，K)J/(kg ⋅；z y x λλλ,,为导热系数，)k m w ⋅；T 为温度，℃；t 为时间，s ；Q 为内热源密度，w/m 3。
对于各向同性材料，不同方向上的导热系数相同，热传导方程可写为以下形式，Q zTy T x T t T c 222222+∂∂+∂∂+∂∂=∂∂λλλρ（6-2）除了热传导方程，计算物体内部的温度分布，还需要指定初始条件和边界条件。
初始条件是指物体最初的温度分布情况，() z y,x,T T 00t ==（6-3）边界条件是指物体外表面与周围环境的热交换情况。
在传热学中一般把边界条件分为三类。
1）给定物体边界上的温度，称为第一类边界条件。
物体表面上的温度或温度函数为已知，s s T T =或),,,(t z y x T T s s =（6-4）2）给定物体边界上的热量输入或输出，称为第二类边界条件。
已知物体表面上热流密度，s sz z y y x xq n z T n y T n x T =∂∂+∂∂+∂∂)(λλλ或),,,()(t z y x q n zT n y T n x T s sz z y y x x=∂∂+∂∂+∂∂λλλ（6-5）3）给定对流换热条件，称为第三类边界条件。
传热问题的基本方程有限元分析
温度场u分布云图
热流场x方向分布云图
热流场y方向分布云图
➢有限元语言描述文件
为生成该问题有限元计算的所有程序源代码，针对之前的ELAB有限元分析得到的微分方程弱 形式，ELAB软件提供简洁的有限元语言描述文件，包括微分方程描述文件、多物理场描述文件以 及求解命令流控制文件。
针对该问题的有限元描述文件包括heatxy.fde(温度场fde文件)， hfxy.fde(热流场fde文件)， heat.mdi， heat.gcn ✓微分方程描述文件heatxy.fde(温度场fde文件)
V
(nxkx
u x
nyky
u ) ud y
将边界条件代入上式(注意，对于已知温度边界条件，虚位移δu为0，可得 :
V
(c
u t
u
kx
u x
u
x
ky
u y
u )dV
y
Q udV
V
q q0 ud
瞬态热传导有限元分析 ➢工程背景
u t
u
kx
u x
u x
ky
u y
u )dV y
Q udV
V
q q0 ud
未知变量：
DISP u u
未知变量定义微分方程弱形式中 的变量
材料参数：
MATE ek ec q 1.0 1.0 0.0 kx(ky) ρc q
材料参数行对应微分方程弱形式 中的变量(考虑各向同性材料，各
在heatxy.fde给出单元的待求未知量，涉及到的材料参数，单元的形函数表达式，刚度 矩阵表达式和载荷表达式，以及为描述刚度矩阵和载荷向量而自定义的函数。 以下给出微分方程描述文件中与微分方程弱形式对应的部分(详细的解析见《有限元分析基础 和应用》中相关章节)：
有限元第12章 热传导问题
第12章热传导问题1. 引言2. 稳态热传导问题33. 瞬态热传导问题一般格式直接积分法模态叠加法解的稳定性与时间步长选择44. 热应力的计算1.1 典型加工方法中的传热问题焊接汽车各个典型部件的加工方法注塑冲压铸造1.1典型加工方法中的传热问题焊接注塑铸造锻压1.1 典型加工方法中的传热问题注塑1.1 典型加工方法中的传热问题焊接1.1 典型加工方法中的传热问题铸造1.1 典型加工方法中的传热问题锻压冷冲热冲1.1 典型加工方法中的传热问题⏹传热问题广泛出现在材料加工领域⏹温度场与宏观力学性能和微观组织变化关系密切1.2 温度场基本方程微分方程边界条件初始条件1.2 温度场基本方程退化为二维问题1.2 温度场基本方程退化为稳态问题稳态热传导问题以前各章所讨论的弹性静力学问题相同，采用C0型插值函数的有限单元进行离散以后，可以直接得到有限元求解方程。
瞬态热传导问题，在空间域有限元离散后，得到的是一阶常微分方程组，不能对它直接求解。
如何进行求解，原则上和下—章将讨论的动力学问题类同，可以采用模态叠加法或直接积分法。
热能传递的三种基本方式：1.2 温度场基本方程热能传递的三种基本方式：热对流：是指由于流体的宏观运动而引起的流体各部分之间发生相对位移，冷、热流体相互掺混所导致的热量传递过程。
热对流仅能发生在流体中。
包括自然对流与强制对流，前者是由于流体冷、热各部分的密度不同而引起的；括自然对流与强制对流前者是于流体冷热各部分的密度不同而引起的；后者是由于水泵、风机或其他压差作用所造成的。
Th q ∆=牛顿冷却公式为表面换热系数，不仅取决于流体物性，以及表面形状等，还与流体速度有密切关系。
h1.2 温度场基本方程1.2 温度场基本方程热能传递的三种基本方式：热辐射：物体通过电磁波来传递能量的方式称为辐射。
物体会因各种原因发出辐射能，其中因热的原因而发出辐射能的现象称为热辐射。
Stefan-Boltzmann定理其中为热力学温度（K），为环境温度，为Stefan-Boltzmann常量i i it)理想黑体其值等于般量。
三维传热方程求解python
三维传热方程求解python传热是热力学中的一个重要概念，它描述了热量如何通过传导、对流和辐射等方式从一个物体传递到另一个物体。
在工程领域中，我们经常需要对传热过程进行分析和计算，以设计高效的热交换设备或优化能源利用。
在三维传热问题中，我们需要求解三维传热方程来描述热量在三维空间中的传递。
三维传热方程是一个偏微分方程，它的求解需要借助数值方法和计算工具。
Python作为一种强大的科学计算工具，可以用来求解三维传热方程。
三维传热方程可以写作：∂T/∂t = α(∂²T/∂x² + ∂²T/∂y² + ∂²T/∂z²)其中，T是温度，t是时间，x、y和z分别是空间坐标，α是热扩散系数。
为了求解这个方程，我们需要确定边界条件和初始条件。
边界条件描述了热量在物体表面的传递，初始条件描述了初始温度分布。
现在我们来看一个具体的例子，假设我们有一个长方体物体，边界上的温度分布为固定值，初始温度分布为均匀的初值。
我们的目标是求解物体内部的温度分布随时间的演化。
我们需要将三维空间离散化为网格。
我们可以使用numpy库来创建一个三维数组来表示网格。
假设我们将空间分为Nx、Ny和Nz 个网格点，那么我们可以创建一个形状为(Nx, Ny, Nz)的数组来表示温度分布。
然后，我们可以使用数值方法来逐步求解三维传热方程。
一个常用的数值方法是显式差分法，它基于时间步进和空间离散化的思想。
具体来说，我们可以使用下面的迭代公式来更新温度分布：T_new[i, j, k] = T_old[i, j, k] + αΔt(Δ²xT[i+1, j, k] - 2Δ²xT[i, j, k] + Δ²xT[i-1, j, k] + Δ²yT[i, j+1, k] - 2Δ²yT[i, j, k] + Δ²y T[i, j-1, k] + Δ²zT[i, j, k+1] - 2Δ²zT[i, j, k] + Δ²zT[i, j, k-1])其中，T_new是更新后的温度分布，T_old是上一步的温度分布，Δt、Δx、Δy和Δz分别是时间步长和空间步长。
三维稳态导热问题数值求解 实验内容的ppt讲解
第一类边界条件
tAB=ti,M=200℃； tBC=tL,j=100℃ ； tCD=ti,0=50℃； tDA=t0,j=50℃；
上边界200℃
其余边界50℃
导热系数为常数、稳态、无内热源时的导 热微分方程式（控制方程）为
t t t 2 2 0 0℃； tbob=50℃； tlb=50℃； trb=50℃； tfb=50℃； tbab=50℃；
划分为10×10×10的三维网格后，Δx=Δy=Δz
实验一、三维稳态导热问题数值求解
一正方体金属块，其长宽 高均为0.1m,上边界面温 度为200℃，其他5个面 温度均为50℃，利用C语 言在10×10×10三维网 格上编写该三维稳态导热 问题计算程序，并求出图 2中所示中间面的温度分 布。
实验一、三维稳态导热问题数值求解
一正方体金属块，其长宽高均为0.1m,上边界面温度 为200℃，其他5个面温度均为50℃，利用C语言在 10×10×10三维网格上编写该三维稳态导热问题计 算程序，并求出图中所示中间面的温度分布。
a /( c) ，称为热扩散率。 式中，
对于二维问题，导热微分方程式为
2 2 t t t a x 2 y 2
初始条件（时间条件）
t 0 t0 50 C
第一类边界条件
t
tAB=200℃； tBC=100℃； tCD=50℃； tDA=50℃；
划分为10×5的二维网格后，Δx=Δy 内节点离散方程（显式差分格式）
t i , j Fo t i 1, j t i 1, j t i, j 1 t i , j 1 1 4Fot i , j
k 1 k k k k k
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
元计算有限元自动生成系统所开发源代码系列各向同性线性三维稳态传热有限元计算程序1.简介元计算（）公司所开发的并行有限元程序自动生成系统（pFEPG）可根据用户需要开发出各种有限元计算程序源代码。
该源代码系列即为pFEPG所开发出来的求解各学科典型问题的有限元计算程序。
该组程序为各向同性线性三维稳态传热有限元计算程序。
2.starta.for，对温度场的数据进行初始化；implicit real*8 (a-h,o-z)character*12 fname,filename(20)common /aa/ ia(250000000)common /bb/ ib(125000000)c.... open disp0 file to get the numbers of nodes and degree of freedomc.... knode .... number of nodes, kdgof .... number ofd.o.f.open(1,file=' ',form='unformatted')read(1) knode,kdgofclose(1)kvar=knode*kdgofwrite(*,*) 'knode,kdgof,kvar ='write(*,'(1x,4i7)') knode,kdgof,kvarkvar1=kvar+1kcoor=3kelem=31250000knb1=kdgof*knode*1if (knb1/2*2 .lt. knb1) knb1=knb1+1kna4=kcoor*knode*2kna1=kdgof*knode*2kna2=kdgof*knode*2kna3=kdgof*knode*2kna5=knode*1if (kna5/2*2 .lt. kna5) kna5=kna5+1knb4=kelem*1if (knb4/2*2 .lt. knb4) knb4=knb4+1knb2=kvar1*1if (knb2/2*2 .lt. knb2) knb2=knb2+1knb3=kvar1*1if (knb3/2*2 .lt. knb3) knb3=knb3+1kna0=1kna1=kna1+kna0kna2=kna2+kna1kna3=kna3+kna2kna4=kna4+kna3kna5=kna5+kna4if (kna5-1.gt.250000000) thenwrite(*,*) 'exceed memory of array ia'write(*,*) 'memory of ia = 250000000'write(*,*) 'memory needed = ',kna5,' in prgram start'stop 55555endifknb0=1knb1=knb1+knb0knb2=knb2+knb1knb3=knb3+knb2knb4=knb4+knb3if (knb4-1.gt.125000000) thenwrite(*,*) 'exceed memory of array ib'write(*,*) 'memory of ib = 125000000'write(*,*) 'memory needed = ',knb4,' in prgram start'stop 55555endifcall start(knode,kdgof,kcoor,kvar,*kelem,maxt,kvar1,ia(kna0),ia(kna1),ia(kna2),*ia(kna3),ia(kna4),ib(knb0),ib(knb1),ib(knb2),*ib(knb3),*filename)endsubroutine start(knode,kdgof,kcoor,kvar,*kelem,maxt,kvar1,u0,u1,u2,*coor,inodvar,nodvar,numcol,lm,node,*filename)implicit real*8 (a-h,o-z)character*12 filename(20)DIMENSION NODV AR(KDGOF,KNODE),COOR(KCOOR,KNODE),R(3),* U0(KDGOF,KNODE),U1(KDGOF,KNODE),U2(KDGOF,KNODE),* INODV AR(KNODE),node(kelem)DIMENSION NUMCOL(KV AR1),LM(KV AR1)CHARACTER*1 MATERIALlogical filflgC .................................................................C ..... KDGOF NUMBER OF D.O.FC ..... KNODE NUMBER OF NODESC ..... INODV AR ID DATAC ..... NODV AR DENOTE THE EQUA TION NUMBER CORRESPONDING THE D.O.F C ..... U0 U1 U2 INITIAL V ALUEC ..... COOR COORDINA TESC ..... NODE ELEMENT NODAL CONNECTIONC .................................................................6 FORMAT (1X, 15I4)7 FORMAT (1X,8F9.3)C.......OPEN ID fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ (1) NUMNOD,NODDOF,((NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE (1)call chms(kdgof,knode,NODV AR)c WRITE(*,*) 'NUMNOD =',NUMNOD,' NODDOF =',NODDOFc WRITE (*,*) 'ID ='c WRITE (*,6) ((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)C..... GET THE NATURAL NODAL ORDERDO 12 N=1,KNODEINODV AR(N)=N12 CONTINUEC..... OPEN ORDER.NOD FILE AND READ THE NODAL ORDER IF THE FILE EXISTinquire(file='ORDER.NOD',exist=filflg)if (filflg) thenOPEN (1,FILE='ORDER.NOD',FORM='UNFORMA TTED',STATUS='OLD')READ (1) (INODV AR(I),I=1,NUMNOD)CLOSE(1)WRITE(*,*) 'NODORDER ='WRITE(*,6) (INODV AR(I),I=1,NUMNOD)endifC..... GET NV BY IDNEQ=0DO 20 JNOD=1,NUMNODJ=INODV AR(JNOD)DO 18 I=1,NODDOFIF (NODV AR(I,J).NE.1) GOTO 18NEQ = NEQ + 1NODV AR(I,J) = NEQ18 CONTINUE20 CONTINUEDO 30 JNOD=1,NUMNODJ=INODVAR(JNOD)DO 28 I=1,NODDOFIF (NODVAR(I,J).GE.-1) GOTO 28N = -NODV AR(I,J)-1NODV AR(I,J) = NODV AR(I,N)28 CONTINUE30 CONTINUEC..... OPEN AND WRITE THE NV FILEOPEN(8,STATUS='unknown',FILE=' ' ,FORM='UNFORMA TTED')WRITE(8) ((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE(8)c WRITE(*,*) 'NUMNOD =',NUMNOD,' NODDOF =',NODDOFc WRITE(*,6) ((NODV AR(I,J),I=1,NODDOF),J=1,NUMNOD)C.... WRITE THE BOUNDAY CONDITION FILE BFD ACCORDING TO THE DISP0 FILE C....OPEN DISP0 FILEOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ(1) NUMNOD,NODDOF,((U0(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE(1)C....OPEN BFD FILEOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')WRITE(1) ((U0(I,J),I=1,NODDOF),J=1,NUMNOD)CLOSE(1)C...... GET THE INITIAL TIME FROM TIME0 FILEC.......OPEN TIME0 FileOPEN(1,FILE=' ',FORM='FORMA TTED')READ(1,*) T0,TMAX,DTTIME = T0IT = 0WRITE(*,*) ' TMAX,DT,TIME,IT =',TMAX,DT,TIME,ITCLOSE(1)C.......OPEN TIME FileOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')WRITE(1) TMAX,DT,TIME,ITCLOSE(1)C.......OPEN COOR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ (1) NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)CLOSE(1)c WRITE(*,*) 'COOR ='c WRITE(*,7) ((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)C...... GET THE INITIAL V ALUE FROM THE DATA FILES BY PREPROCESSORinquire(file='disp1',exist=filflg)if (filflg) thenopen(16,file='disp1',form='unformatted',status='old')read(16) numnod,noddof,((U0(J,N),J=1,NODDOF),N=1,NUMNOD) close(16)endifinquire(file='disp2',exist=filflg)if (filflg) thenopen(16,file='disp2',form='unformatted',status='old')read(16) numnod,noddof,((U1(J,N),J=1,NODDOF),N=1,NUMNOD) close(16)endifinquire(file='disp3',exist=filflg)if (filflg) thenopen(16,file='disp3',form='unformatted',status='old')read(16) numnod,noddof,((U2(J,N),J=1,NODDOF),N=1,NUMNOD) close(16)endifc WRITE(*,*) ' U0 = 'c WRITE(*,'(6F13.3)') ((U0(J,N),J=1,NODDOF),N=1,NUMNOD) C WRITE(*,*) ' U1 = 'C WRITE(*,'(6F13.3)') ((U1(J,N),J=1,NODDOF),N=1,NUMNOD)C...... COMPUTE THE INITIAL V ALUE BY BOUND.FORzo = 0.0d0c DO 321 N=1,NUMNODc DO 100 J=1,NCOORc100 R(J) = COOR(J,N)c DO 200 J=1,NODDOFc U0(J,N) = BOUND(R,zo,J)c U1(J,N) = BOUND1(R,zo,J)c U2(J,N) = BOUND2(R,zo,J)c200 CONTINUEc321 CONTINUEC.......OPEN AND WRITE THE INITIAL VALUE FILE UNODOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown')WRITE(1) ((U0(I,J),J=1,NUMNOD),I=1,NODDOF),* ((U1(I,J),J=1,NUMNOD),I=1,NODDOF),* ((U2(I,J),J=1,NUMNOD),I=1,NODDOF),* ((U0(I,J),J=1,NUMNOD),I=1,NODDOF)CLOSE (1)c.... open IO fileopen(21,file=' ',form='formatted',status='old')read(21, '(1a)') materialread(21,*) numtypclose(21)DO I=1,NEQNUMCOL(i)=1ENDDOC.......OPEN ELEM0 fileOPEN (3,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')NUMEL=0KELEM=0KEMATE=0DO 2000 ITYP=1,NUMTYPC.......INPUT ENODEREAD (3) NUM,NNODE,* ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODEcc WRITE(*,*) 'NODE ='cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) IF (KELEM.LT.NUM*NNODE) KELEM = NUM*NNODENNE = NNODEIF (MATERIAL.EQ.'Y' .OR. MATERIAL.EQ.'y') THENREAD (3) MMATE,NMATEIF (KEMATE.LT.MMATE*NMATE) KEMATE = MMATE*NMATENNE = NNE-1ENDIFWRITE(*,*) 'MMATE =',MMATE,' NMATE =',NMATEcc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODEcc WRITE(*,*) 'NODE ='cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) DO 1000 NE=1,NUML=0DO 700 INOD=1,NNENODI=NODE((NE-1)*NNODE+INOD)DO 600 IDGF=1,KDGOFINV=NODV AR(IDGF,NODI)IF (INV.LE.0) GOTO 600L=L+1LM(L)=INV600 CONTINUE700 CONTINUENUMEL=NUMEL+1C WRITE (*,*) 'L,LM =',LC WRITE (*,'(1X,15I5)') (LM(I),I=1,L)if (l.gt.0) call ACLH(NEQ,NUMCOL,l,lm)1000 continue2000 CONTINUEc CLOSE(1)CLOSE(3)call BCLH(NEQ,NUMCOL)MAXA=NUMCOL(NEQ)C.......OPEN SYS FileOPEN (2,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown') WRITE(2) NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATECLOSE (2)OPEN(2,FILE=' ',FORM='UNFORMATTED',STATUS='unknown')write(2) (NUMCOL(I),I=1,NEQ)CLOSE(2)c write(*,*) 'NEQ,NUMCOL=',NEQc write(*,6) (NUMCOL(i),i=1,NEQ)ENDsubroutine chms(kdgof,knode,id)dimension id(kdgof,knode),ms(1000),is(1000)do 1000 k=1,kdgofm = 0do 800 n=1,knodeif (id(k,n).le.-1) id(k,n)=-1if (id(k,n).le.1) goto 800j=id(k,n)j0=0if (m.gt.0) thendo i=1,mif (j.eq.ms(i)) j0=is(i)enddoendifif (j0.eq.0) thenm=m+1ms(m)=jis(m)=nid(k,n)=1elseid(k,n)=-j0-1endif800 continue1000 continuereturnendSUBROUTINE ACLH(NEQ,NUMCOL,ND,LM)implicit real*8 (a-h,o-z)DIMENSION LM(ND),NUMCOL(NEQ)LS=LM(1)+1DO 100 I=1,ND110 IF(LM(I)-LS) 120,100,100120 LS=LM(I)100 CONTINUEDO 200 I=1,NDII=LM(I)ME=II-LSIF(ME.GT.NUMCOL(II)) NUMCOL(II)=ME200 CONTINUERETURNENDSUBROUTINE BCLH(NEQ,NUMCOL)implicit real*8 (a-h,o-z)DIMENSION NUMCOL(NEQ)C NUMCOL(1) = 1DO 490 I=2,NEQ490 NUMCOL(I) = NUMCOL(I) + NUMCOL(I-1) + 1RETURNEND3.ewenre3da.for，Galerkin法求解温度场的主程序implicit real*8 (a-h,o-z)character*12 fname,filename(20)common /aa/ ia(250000000)common /bb/ ib(125000000)common /cc/ ic(62500000)open(1,file=' ',form='unformatted',status='old')read(1) knode,kdgofclose(1)MAXT=250000000/2/2C.......OPEN SYS FileOPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')read(2) NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATECLOSE (2)IF (MAXA.GT.MAXT) THENWRITE(*,*) 'MATRIX A EXCEED CORE MEMERY .... ',MAXA WRITE(*,*) 'REQUIRED CORE MEMERY ........... ',MAXTSTOP 0000ENDIFKV AR=KNODE*KDGOFKCOOR=3C KELEM=31250000WRITE(*,*) 'KNODE,KDGOF,KV AR,KCOOR,KELEM ='WRITE(*,'(1X,6I7)') KNODE,KDGOF,KV AR,KCOOR,KELEMkna1=kdgof*knode*1if (kna1/2*2 .lt. kna1) kna1=kna1+1knc1=kdgof*knode*2knc2=kcoor*knode*2knc7=kdgof*knode*2knc3=neq*2knb1=maxa*2knb2=maxa*2kna2=neq*1if (kna2/2*2 .lt. kna2) kna2=kna2+1knc6=kemate*2kna3=kelem*1if (kna3/2*2 .lt. kna3) kna3=kna3+1knc8=100000*2knc5=neq*2knc4=kdgof*knode*2kna0=1kna1=kna1+kna0kna2=kna2+kna1kna3=kna3+kna2if (kna3-1.gt.125000000) thenwrite(*,*) 'exceed memory of array ib'write(*,*) 'memory of ib = 125000000'write(*,*) 'memory needed = ',kna3,' in prgram ewenre3da'stop 55555endifknb0=1knb1=knb1+knb0knb2=knb2+knb1if (knb2-1.gt.250000000) thenwrite(*,*) 'exceed memory of array ia'write(*,*) 'memory of ia = 250000000'write(*,*) 'memory needed = ',knb2,' in prgram ewenre3da'stop 55555endifknc0=1knc1=knc1+knc0knc2=knc2+knc1knc3=knc3+knc2knc4=knc4+knc3knc5=knc5+knc4knc6=knc6+knc5knc7=knc7+knc6knc8=knc8+knc7if (knc8-1.gt.62500000) thenwrite(*,*) 'exceed memory of array ic'write(*,*) 'memory of ic = 62500000'write(*,*) 'memory needed = ',knc8,' in prgram ewenre3da'stop 55555endifcall ewenre3da(knode,kdgof,kvar,kcoor,*numtyp,numel,neq,kelem,kemate,maxa,*maxt,neq1,ib(kna0),ib(kna1),ib(kna2),*ia(knb0),ia(knb1),ic(knc0),ic(knc1),ic(knc2),*ic(knc3),ic(knc4),ic(knc5),ic(knc6),ic(knc7),*filename)endsubroutine ewenre3da(knode,kdgof,kvar,kcoor,*numtyp,numel,neq,kelem,kemate,maxa,*maxt,neq1,nodvar,jdiag,node,a,*b,u,coor,f,ubf,u1,*emate,eu,sml,*filename)implicit real*8 (a-h,o-z)character*12 filename(20)DIMENSION NODV AR(KDGOF,KNODE),U(KDGOF,KNODE),COOR(KCOOR,KNODE), *eu(kdgof,knode),& F(NEQ),A(MAXA),B(MAXA),JDIAG(NEQ),EMATE(KEMATE),& NODE(KELEM),SML(100000),u1(neq),UBF(KDGOF,KNODE)6 FORMA T (1X,15I5)7 FORMA T (1X,5e15.5)1001 FORMA T(1X,9I7)C.......OPEN TIME FileOPEN(1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ(1) TMAX,DT,TIME,ITWRITE(*,*) ' TMAX,DT,TIME,IT =',TMAX,DT,TIME,ITCLOSE(1)C.......OPEN NODVAR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ (1) ((NODV AR(I,J),I=1,KDGOF),J=1,KNODE)CLOSE (1)cc WRITE(*,*) 'KDGOF =',KDGOF,' KNODE =',KNODEcc WRITE (*,*) 'NODV AR ='cc WRITE (*,6) ((NODV AR(I,J),I=1,KDGOF),J=1,KNODE)C.......OPEN COOR fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ (1) NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD) CLOSE(1)cc WRITE(*,*) 'NUMNOD,NCOOR=',NUMNOD,NCOORC.......OPEN BFD fileOPEN (1,FILE=' ',FORM='UNFORMATTED',STA TUS='OLD')READ (1) ((UBF(J,I),J=1,KDGOF),I=1,KNODE)CLOSE (1)cc WRITE (*,*) 'BF ='cc WRITE(*,7) ((UBF(J,I),J=1,KDGOF),I=1,KNODE)numtyp = 2C.......OPEN DIAG fileOPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')READ(2) (JDIAG(I),I=1,NEQ)CLOSE(2)C.......OPEN ELEM0 fileOPEN (3,FILE=' ',FORM='UNFORMATTED',STATUS='OLD')itime=01 continueitime=itime+1if (itime.gt.1) thenwrite(*,*) 'Nonlinear Iteration Times ========',itimerewind(3)endifDO 111 I=1,KNODEDO 111 J=1,KDGOFU(J,I) = UBF(J,I)111 CONTINUEcc WRITE (*,*) 'BF ='cc WRITE(*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)DO 112 I=1,MAXAA(I) = 0.0B(I) = 0.0112 CONTINUEDO 2300 I=1,NEQ2300 CONTINUENUMEL=0C.......OPEN EMATE+ENODE+ELOAD fileC OPEN (3,FILE=' ',FORM='UNFORMA TTED',STATUS='OLD')DO 2000 ITYP=1,NUMTYPC.......INPUT ENODEREAD (3) NUM,NNODE,* ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODEcc WRITE(*,*) 'NODE ='cc WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) NNE = NNODEnne = nne-1K=0DO 115 J=1,NNEJNOD = NODE(J)DO 115 L=1,KDGOFIF (NODV AR(L,JNOD).NE.0) K=K+1115 CONTINUEWRITE(*,*) 'K =',Kkk=k*kk0=1k1=k0+k*kk2=k1+kk3=k2+kk4=k3+k*kk5=k4+k*kCALL ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE, * ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODV AR,COOR,NODE,EMATE, & A,B,JDIAG,&sml(k0),sml(k1),sml(k2),sml(k3),sml(k4),&eu,*U)2000 CONTINUEDO 2050 IJ=1,NEQif (itime.le.1) u1(IJ) = 0.0F(IJ)=0.0D02050 CONTINUEDO 2200 I=1,KNODEDO 2100 J=1,KDGOFIJ=NODV AR(J,I)IF (IJ.LE.0) GOTO 2100F(IJ)=F(IJ)+U(J,I)U1(IJ)=F(IJ)2100 CONTINUE2200 CONTINUECC IF (IT.GT.0) THENcc WRITE (*,*) 'U ='cc WRITE (*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE)cc WRITE (*,*) 'NEQ =',NEQ,' F ='cc WRITE(*,7) (F(I),I=1,NEQ)if (itime.le.1) thenC.......OPEN LMATRIX FILEOPEN (2,FILE=' ',FORM='UNFORMATTED',STA TUS='unknown')CLOSE (2)endifWRITE(*,*) 'NIN_SOLVER MEMORY REQUIRED .... ',MAXAIF (MAXA.GT.MAXT) THENWRITE(*,*) 'WARNING MA TRIX A EXCEED CORE MEMORY .... ',MAXT c STOP 0000ENDIFCALL REDU(A,B,U1,JDIAG,NEQ,MAXA,1)C WRITE(*,*) ' U1 = 'C WRITE(*,7) (A(I),I,MAXA)C WRITE(*,7) (F(I),I=1,NEQ)NOUT = 20OPEN(NOUT,FILE=' ',FORM='FORMATTED',STA TUS='unknown')DO 3200 INOD=1,KNODEDO 3100 IDFG=1,KDGOFN=NODV AR(IDFG,INOD)C WRITE (*,*) 'N =',Nif(n.le.0) theneu(IDFG,INOD)=u(IDFG,INOD)elseeu(IDFG,INOD)=u1(N)endif3100 CONTINUE3200 CONTINUEDO 3400 N=1,KNODEWRITE (NOUT,3600) N,(eu(I,N),I=1,KDGOF)3400 CONTINUE3600 FORMA T (1X,I5,1X,6E11.4,9(/6X,6E11.4))CLOSE (NOUT)open(10,file='unod',form='unformatted',status='unknown')write(10) ((eu(j,i),i=1,knode),j=1,kdgof)close(10)close (3)RETURNENDSUBROUTINE ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE, *ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODV AR,COOR,NODE,EMATE,&A,B,JDIAG,*es,em,ef,Estifn,Estifv,eu,*U)implicit real*8 (a-h,o-z)DIMENSION NODV AR(KDGOF,KNODE),COOR(KCOOR,KNODE),NODE(KELEM), *U(KDGOF,KNODE),EMA TE(300),&A(MAXa),B(MAXa),JDIAG(neq),*es(k,k),em(k),ef(k),eu(kdgof,knode),*Estifn(k,k),Estifv(kk),*R(500),PRMT(500),COEF(500),LM(500)17 FORMA T (1X,15I5)18 FORMA T (1X,8e9.2)READ (3) MMATE,NMATE,((EMATE((I-1)*NMATE+J),J=1,NMATE),* I=1,MMATE)WRITE(*,*) 'MMATE =',MMATE,' NMATE =',NMATEWRITE (*,*) 'EMATE ='WRITE (*,18) ((EMATE((I-1)*NMATE+J),J=1,NMATE),* I=1,MMATE)DO 1000 NE=1,NUMNR=0DO 130 J=1,NNEJNOD = NODE((NE-1)*NNODE+J)IF (JNOD.LT.0) JNOD = -JNODPRMT(NMA TE+7+J) = JNODDO 120 I=1,NCOORNR=NR+1120 R(NR) = COOR(I,JNOD)130 CONTINUEIMA TE = NODE(NNODE*NE)DO 140 J=1,NMATE140 PRMT(J) = EMA TE((IMATE-1)*NMATE+J) PRMT(NMATE+1)=TIMEPRMT(NMATE+2)=DTPRMT(NMA TE+3)=IMATEprmt(NMA TE+4)=NEprmt(NMA TE+5)=NUMprmt(NMA TE+6)=ITprmt(NMA TE+7)=NMATEprmt(NMA TE+8)=ITIMEprmt(NMA TE+9)=ITYPgoto (1,2), ityp1 call aec8g2(r,coef,prmt,es,em,ec,ef,ne)goto 32 call agq4g2(r,coef,prmt,es,em,ec,ef,ne)goto 33 continueC WRITE(*,*) 'ES EM EF ='C DO 555 I=1,KC555 WRITE(*,18) (ES(I,J),J=1,K)C WRITE(*,18) (EM(I),I=1,K)C WRITE(*,18) (EF(I),I=1,K)CC IF (IT.GT.0) THENdo 201 i=1,kdo 201 j=1,kEstifn(i,j)=0.0201 continuedo 202 i=1,kEstifn(i,i)=Estifn(i,i)do 202 j=1,kEstifn(i,j)=Estifn(i,j)+es(i,j)202 continueL=0M=0I=0DO 700 INOD=1,NNENODI=NODE((NE-1)*NNODE+INOD)DO 600 IDGF=1,KDGOFINV=NODV AR(IDGF,NODI)IF (INV.EQ.0) GOTO 600I=I+1IF (INV.LT.0) GOTO 305L=L+1LM(L)=INVU(IDGF,NODI)=U(IDGF,NODI)*+ef(i)305 J=0DO 500 JNOD=1,NNENODJ=NODE((NE-1)*NNODE+JNOD)DO 400 JDGF=1,KDGOFJNV=NODV AR(JDGF,NODJ)IF (JNV.EQ.0) GOTO 400J=J+1IF (JNV.LT.0) GOTO 400IF (INV.LT.0) GOTO 310M=M+1Estifv(m)=Estifn(i,j)310 CONTINUEIF (INV.LT.0)* U(JDGF,NODJ)=U(JDGF,NODJ)-ESTIFN(I,J)*U(IDGF,NODI) 400 CONTINUE500 CONTINUE600 CONTINUE700 CONTINUEC WRITE (*,*) 'U ='C WRITE (*,18) ((U(J,I),J=1,KDGOF),I=1,KNODE)LRD=MNER=NUMEL+NEC WRITE(*,*) '**************************'C WRITE(*,*) (ESTIFV(I),I=1,LRD)C WRITE (*,*) 'Einform ............'C WRITE (*,'(1X,15I5)') L,LRD,(LM(I),I=1,L)DO 800 I=1,LJ=LM(I)800 CONTINUEcall ADDA(A,B,JDIAG,L,LM,ESTIFV,NEQ,MAXA)1000 CONTINUERETURNENDSUBROUTINE ADDA(A,B,JDIAG,ND,LM,ESTIF,NEQ,MAXA)implicit real*8 (a-h,o-z)DIMENSION A(MAXA),B(MAXA),JDIAG(NEQ),LM(ND),ESTIF(ND,ND) C WRITE (*,*) ND, (LM(I),I=1,ND)C WRITE (*,*) ((ESTIF(I,J),J=1,ND),I=1,ND)DO 300 I=1,NDII = LM(I)DO 280 J=1,IJJ = LM(J)IF (II.LT.JJ) GOTO 240K = JDIAG(II) - II + JJA(K) = A(K) + ESTIF(I,J)B(K) = B(K) + ESTIF(J,I)GOTO 280240 K = JDIAG(JJ) - JJ + IIB(K) = B(K) + ESTIF(I,J)A(K) = A(K) + ESTIF(J,I)280 CONTINUE300 CONTINUERETURNENDSUBROUTINE REDU(A,B,U,JDIAG,NEQ,MAXA,KKK)implicit real*8 (a-h,o-z)DIMENSION A(MAXA),B(MAXA),JDIAG(NEQ),U(NEQ)DOUBLE PRECISION CPUTE L U & u, L*U(T) = A, A*u = fDO 500 I=2,NEQI1 = I-1NI = JDIAG(I)LI = I-NI+JDIAG(I-1)+1IF (KKK.GT.1) GOTO 333DO 200 J=LI,INJ = JDIAG(J)LJ = J-NJ+1IF (J.GT.1) LJ = LJ+JDIAG(J-1)LIJ = MAX0(LI,LJ)J1 = J-1IF (J.EQ.I) GOTO 130C = 0.0D0DO 100 L=LIJ,J1100 C = C + A(NI-I+L)*B(NJ-J+L)A(NI-I+J) = (A(NI-I+J) - C)/B(NJ)130 C = 0.0D0DO 150 L=LIJ,J1150 C = C + B(NI-I+L)*A(NJ-J+L)200 B(NI-I+J) = B(NI-I+J) - C333 C = 0.0D0DO 400 L=LI,I1400 C = C + A(NI-I+L)*U(L)U(I) = U(I)-C500 CONTINUEJ = MAXA + 1N = NEQ + 1700 N = N - 1J = J - 1U(N) = U(N)/B(J)IF (N.EQ.1) GOTO 999M = J-JDIAG(N-1)-1DO 800 I=1,MJ = J - 1800 U(N-I) = U(N-I) - B(J)*U(N)GOTO 700999 RETURNEND3.1.aec8g2.for，计算区域内单元刚度矩阵和荷载向量的子程序subroutine aec8g2(coorr,coefr,& prmt,estif,emass,edamp,eload,num)c .... coorr ---- nodal coordinate valuec .... coefr ---- nodal coef valueimplicit real*8 (a-h,o-z)dimension estif(8,8),elump(8),emass(8),& eload(8)dimension prmt(*),& egux(8),eguy(8),eguz(8),& coorr(3,8),coor(3)common /raec8g2/ru(8,32),& cu(8,4)c .... store shape functions and their partial derivativesc .... for all integral pointscommon /vaec8g2/rctr(3,3),crtr(3,3)common /daec8g2/ refc(3,8),gaus(8),& nnode,ngaus,ndisp,nrefc,ncoor,nvar,& nvard(1),kdord(1),kvord(8,1)c .... nnode ---- the number of nodesc .... nrefc ---- the number of numerical integral pointsc .... ndisp ---- the number of unknown functionsc .... nrefc ---- the number of reference coordinatesc .... nvar ---- the number of unknown varibles varc .... refc ---- reference coordinates at integral pointsc .... gaus ---- weight number at integral pointsc .... nvard ---- the number of var for each unknownc .... kdord ---- the highest differential order for each unknown c .... kvord ---- var number at integral points for each unknownek=prmt(1)q=prmt(2)time=prmt(3)dt=prmt(4)imate=prmt(5)+0.5ielem=prmt(6)+0.5nelem=prmt(7)+0.5it=prmt(8)+0.5nmate=prmt(9)+0.5itime=prmt(10)+0.5ityp=prmt(11)+0.5if (num.eq.1) call aec8g2ic .... initialize the basic datado 10 i=1,nvareload(i)=0.0do 10 j=1,nvarestif(i,j)=0.010 continuedo 999 igaus=1,ngauscall aec8g2t(nnode,nrefc,ncoor,refc(1,igaus),coor,coorr, & rctr,crtr,det)c .... coordinate transfer from reference to original systemc .... rctr ---- Jacobi's matrixc .... crtr ---- inverse matrix of Jacobi's matrixx=coor(1)y=coor(2)z=coor(3)rx=refc(1,igaus)ry=refc(2,igaus)rz=refc(3,igaus)iu=(igaus-1)*4+1if (num.gt.1) goto 2c .... the following is the shape function caculationcall aec8g21(refc(1,igaus),ru(1,iu),rctr,crtr)2 continuec .... the following is the shape function transformation c .... from reference coordinates to original coordinatescall shapn(nrefc,ncoor,8,ru(1,iu),cu,crtr,1,4,4)weigh=det*gaus(igaus)do 100 i=1, 8egux(i) = 0.0eguy(i) = 0.0eguz(i) = 0.0100 continuevol = 1.0d0do 101 i=1,8iv=kvord(i,1)stif=+cu(i,2)egux(iv)=egux(iv)+stif101 continuedo 102 i=1,8iv=kvord(i,1)stif=+cu(i,3)eguy(iv)=eguy(iv)+stif102 continuedo 103 i=1,8iv=kvord(i,1)stif=+cu(i,4)eguz(iv)=eguz(iv)+stif103 continuec .... the following is the stiffness computationdo 202 iv=1,8do 201 jv=1,8stif=+egux(iv)*egux(jv)*ek*vol& +eguy(iv)*eguy(jv)*ek*vol& +eguz(iv)*eguz(jv)*ek*volestif(iv,jv)=estif(iv,jv)+stif*weigh201 continue202 continuec .... the following is the load vector computationdo 501 i=1,8iv=kvord(i,1)stif=+cu(i,1)*q*voleload(iv)=eload(iv)+stif*weigh501 continue999 continue998 continuereturnendsubroutine aec8g2iimplicit real*8 (a-h,o-z)common /daec8g2/ refc(3,8),gaus(8),& nnode,ngaus,ndisp,nrefc,ncoor,nvar,& nvard(1),kdord(1),kvord(8,1)c .... initial datac .... refc ---- reference coordinates at integral pointsc .... gaus ---- weight number at integral pointsc .... nvard ---- the number of var for each unknownc .... kdord ---- the highest differential order for each unknown c .... kvord ---- var number at integral points for each unknownngaus= 8ndisp= 1nrefc= 3ncoor= 3nvar = 8nnode= 8kdord(1)=1nvard(1)=8kvord(1,1)=1kvord(2,1)=2kvord(3,1)=3kvord(4,1)=4kvord(5,1)=5kvord(6,1)=6kvord(7,1)=7kvord(8,1)=8refc(1,1)=5.773502692e-001refc(2,1)=5.773502692e-001refc(3,1)=5.773502692e-001gaus(1)=1.000000000e+000refc(1,2)=5.773502692e-001refc(2,2)=5.773502692e-001refc(3,2)=-5.773502692e-001gaus(2)=1.000000000e+000refc(1,3)=5.773502692e-001refc(2,3)=-5.773502692e-001refc(3,3)=5.773502692e-001gaus(3)=1.000000000e+000refc(1,4)=5.773502692e-001refc(2,4)=-5.773502692e-001refc(3,4)=-5.773502692e-001gaus(4)=1.000000000e+000refc(1,5)=-5.773502692e-001refc(2,5)=5.773502692e-001refc(3,5)=5.773502692e-001gaus(5)=1.000000000e+000refc(1,6)=-5.773502692e-001refc(2,6)=5.773502692e-001refc(3,6)=-5.773502692e-001gaus(6)=1.000000000e+000refc(1,7)=-5.773502692e-001refc(2,7)=-5.773502692e-001refc(3,7)=5.773502692e-001gaus(7)=1.000000000e+000refc(1,8)=-5.773502692e-001refc(2,8)=-5.773502692e-001refc(3,8)=-5.773502692e-001gaus(8)=1.000000000e+000endsubroutine aec8g2t(nnode,nrefc,ncoor,refc,coor,coorr,& rc,cr,det)implicit real*8 (a-h,o-z)dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10), & coorr(ncoor,nnode),coor(ncoor)c .... compute coordinate tranformationc .... coorr ---- nodal coordinatesc .... rc ---- jacobi's matrixc .... cr ---- inverse matrix of jacobi's matrixcall taec8g2(refc,coor,coorr,rc)c .... coordinate transfer caculationc .... from reference coordinate refc to origenal coordinate coorn=nrefcm=n*2det = 1.0do 10 i=1,ndo 10 j=1,nif (i.le.ncoor) a(i,j) = rc(i,j)if (i.gt.ncoor) a(i,j)=1.0a(i,n+j)=0.0if (i.eq.j) a(i,n+i) = 1.010 continuec write(*,*) 'a ='c do 21 i=1,nc21 write(*,8) (a(i,j),j=1,m)do 400 i=1,namax = 0.0l = 0do 50 j=i,nc = a(j,i)if (c.lt.0.0) c = -cif (c.le.amax) goto 50amax = cl = j50 continuedo 60 k=1,mc = a(l,k)a(l,k) = a(i,k)a(i,k) = c60 continuec = a(i,i)det = c*detdo 100 k=i+1,m100 a(i,k) = a(i,k)/cdo 300 j=1,nif (i.eq.j) goto 300do 200 k=i+1,m200 a(j,k) = a(j,k)-a(i,k)*a(j,i)c write(*,*) 'i =',i,' j =',j,' a =' c do 11 ii=1,nc11 write(*,8) (a(ii,jj),jj=1,m) 300 continue400 continuedo 500 i=1,nrefcdo 500 j=1,ncoor500 cr(i,j) = a(i,n+j)c write(*,*) 'a ='c do 22 i=1,nc22 write(*,8) (a(i,j),j=1,m)c write(*,*) 'rc ='c do 24 i=1,ncoorc24 write(*,8) (rc(i,j),j=1,nrefc)c write(*,*) 'cr ='c do 23 i=1,nrefcc23 write(*,8) (cr(i,j),j=1,ncoor)c write(*,*) 'det =',detif (det.lt.0.0) det=-detc write(*,*) 'det =',det8 format(1x,6f12.3)endsubroutine aec8g21(refc,shpr,rctr,crtr)c .... compute shape functions and their partial derivativesc .... shapr ---- store shape functions and their partial derivativesimplicit real*8 (a-h,o-z)dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)external faec8g21rx=refc(1)ry=refc(2)rz=refc(3)call dshap(faec8g21,refc,shpr,3,8,1)c .... shape function and their derivatives computationc .... compute partial derivatives by centered differencec .... which is in the file ccshap.for of FEPG libraryreturnendreal*8 function faec8g21(refc,n)c .... shape function caculationimplicit real*8 (a-h,o-z)common /ccaec8g2/ xa(8),ya(8),za(8)common /vaec8g2/ rctr(3,3),crtr(3,3)dimension refc(3)common /coord/ coor(3),coora(27,3)x=coor(1)y=coor(2)z=coor(3)rx=refc(1)ry=refc(2)rz=refc(3)goto (1,2,3,4,5,6,7,8) n1 faec8g21=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.。