各向同性线性三维稳态传热有限元计算程序

合集下载

传热问题有限元分析

【问题描述】本例对覆铜板模型进行稳态传热以及热应力分析，图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、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)$ 的解。

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.。