程序设计1详细讲解

结果整理（后处理）
• 等值线 • 矢量图 • 变形网格图
编程基本原则
• 结构化，模块化 • 可读性强
变量命名简单直观，适当注释 程序条理清楚，逻辑性强 • 节约内存，计算效率优化
FEAP程序
输入基本变量 D(10,NUMAT) F(NDF,NUMNP) ID(NDF,NUMNP) IE(NUMMAT) IX(NEN1,NUMEL) T(NUMNP) X(NDM, NUMNP)
DO 200 I = 1,2 DO 200 J = 1,2 XS(I,J) = 0.0 D0 200 K = 1,4 200 XS(I,J) = XS(I,J) + XL(I,K)*SHP(J,K)
C C---COMPUTE JACOBIAN DETERMINANT C
XSJ = XS(1,1)*XS(2,2) - XS(1,2)*XS(2,1) C C--- TRANSFORM NATURAL DERIVATIVES TO C--- X,Y DERIVATIVES C
材料性质数组 节点力数组 约束信息数组 各组材料对应的单元类型 单元编号及材料组号 节点温度 节点坐标
NDF NDM NEN NEN1 NUMEL NUMMAT NUMNP
节点自由度数目 问题的维数 单元的节点数 NEN+1 单元总数 材料性质组数 节点总数
COOR ELEM BOUN FORC MATE END
DO 300 I = 1,4 TEMP = (XS(2,2)*SHP(1,I) – XS(2,1)*SHP(2,I))/XSJ SHP(2,I) = (XS(2,2)*SHP(1,I) – XS(2,1)*SHP(2,I))/XSJ 300 SHP(1,I) = TEMP RETURN END
SUBROUTINE PGAUSS(L,LINT,R,Z,W) C C---- GAUSS POINTS AND WEIGHTS FOR TWO
D矩阵 单元刚度矩阵
D11 D12
D12
D11
0
0
0
0
D33
kij BiT DB j
Q j DB j
D11 N j, x
Q j
D12
N
j, x
D33 N j, y
D12 N j, y
D11 N j, y
Q33 N j, x
kij
N i, xQ11 N i, y Q21
有限元程序
• 数据输入(前处理) • 求解计算 • 结果整理与输出(后处理)
数据输入(前处理)
• 单元几何: 单元编号; 节点坐标 • 材料性质 • 边界条件(约束) • 荷载(体积力, 面积力)
求解计算
• 单元刚度矩阵(形函数矩阵，B矩阵，D矩阵，高斯积分) • 单元荷载列阵 • 组集整体刚度矩阵和整体荷载列阵 • 方程组求解
DIMENSIONS C
DATA LR/-1,1,1,-1,0,1,0,-1,0/,LZ/-1,-1,1,1,-1,0,1,0,0/ LW/4/
LINT = L + L GO TO (1,2,3) L C---- 11 INTERGRATION 1 R(1) = 0. Z(1) = 0. W(1) = 4. RETURN
U 22
...
U 2n
. ... .
0
...
U
nn
Ly = r Ua = y
迭代解法 1。共轭梯度法
Min U(a)
a j1 a j j d j
d是搜索方向，是步长
d j1 -U a j1 j d j
j
U a j1 .U a j1 - U a j U a j .U a j
2. 预条件共轭梯度法
Leabharlann Baidu
C---- CONTRIBUTION TO ELEMENT STIFFNESS
C
DO 100 L = 1,LINT
CALL SHAPEF(SG(L),TG(L),XL, XSJ, SHP)
DV = XSJ*WG(L)
D11 = D(1)*DV
D22 = D(2)*DV
D33 = D(3)*DV
C C---- FOR EACH J NODE COMPUTE DB = D*B C
C--- IN NATURAL COORDINATES
C
DO 100 I = 1,4
SHP(3,I) = (0.5+SI(I)*S)*(0.5+TI(I)*T)
SHP(1,I) = SI(I)*(0.5+TI(I)*T) 100 SHP(2,I) = (0.5+SI(I)*S)*TI(I) C C---COMPUTE JACOBIAN TRANSROMATION C---FROM X,Y TO S,T C
C---- 2 2 INTERGRATION 2 G = 1./SQRT(3.)
DO 21 I = 1,4 R(I) = G*LR(I) Z(I) = G*LZ(I) 21 W(I) = 1. RETURN C---- 3 3 INTERGRATION 3 G = 1./SQRT(0.6) H = 1./81. DO 31 I = 1,9 R(I) = G*LR(I) Z(I) = G*LZ(I) 31 W(I) = H*LW(I) RETURN END
Bi 0
N i, y
0
N i, y
N i, x
D11 D12
D12
D11
0
0
0
0
D33
SUBROUTINE SHAPEF(S, T, XL, XSJ, SHP)
C--- shape function routine for 4-node isoparametric quadrilateral
C---- FOR ISOTROPCI LINEAR ELASTICITY
C----PLANE STRESS AND PLANE STRAIN DIFFER ONLY IN
C--- VALUES OF THE CONSTANTS D(1),D(2) ,D(3) SUPPLIED
C
C----FOR PLANE STRESS
以平面四边形等参单元为例 形函数
Ni
1 4
1
i
1
i
x Ni xi
y Ni yi
u Niui
形函数导数
v Nivi
N i, N i,
x ,
x,
y, Ni, x
y
,
N
i,
y
Ni, x
N
i,
y
1 J
y, - x,
- y, Ni,
x,
N
i,
B矩阵 D矩阵
N i, x
C C---COMPUTE LOWER TRIANGULAR PART BY SYMMETRY C
NL = NEL + NEL DO 200 I = 2,NEL DO 200 J = 1, I 200 S(I,J) = S(J,I) RETURN END
线性方程组求解 特点:
稀疏 对称 正定
解法: 直接解法 迭代解法
C--- D(1) = E/(1.-NU*NU)
C---- D(2) = NU*D(1)
C---- D(3) = E/(2.*(1.+NU))
C---- DV IS AN AREA WEIGHTING
C---- LINT IS THE NUMBER FO INTEGRATION POINTS
C---- NEL IS THE NUMBER OF NODES ON THE ELEMENT
C XSJ = Jacobian determinant
C
DIMENSION SHP(3,4), XL(2,4), XS(2,2), SI(4), TI(4)
DATA SI,TI/-0.5, 0.5, 0.5, -0.5, -0.5, -0.5, 0.5, 0.5 /
C
C---COMPUTE SHAPE FUNCTIONS AND DERIVATIVES
DO 100 I = 1,J S(I+I-1,J+J-1) = S(I+I-1,J+J-1) + SHP(1,I)*DB11+SHP(2,I)*DB31 S(I+I-1,J+J-1) = S(I+I-1,J+J ) + SHP(1,I)*DB12+SHP(2,I)*DB32 S(I+I ,J+J-1) = S(I+I ,J+J-1) + SHP(1,I)*DB31+SHP(2,I)*DB21 100 S(I+I , J+J) = S(I+I ,J+J ) + SHP(1,I)*DB32SHP(2,I)*DB22
条件数：最大特征值与最小特征值之比
条件数越大，迭代收敛越困难
预条件目的：减小刚度矩阵的条件数
BKa =Br
P = BK
c = Br
Pa = c
C---- S
IS THE ELMENT
C---- SG,TG ARE INTERGRATION POINTS IN
C----
NATRUAL COORDINATES
C---- WG IS THE INTERGARION WEIGHT
C
C---- FOR EACH INTERGRATION COMPUTE
宏命令
TANG FORM SOLV DISP STRE
计算整体刚度矩阵 形成线性方程组的右端项 求解方程 输出位移 输出应力
TANG
计算整体刚度矩阵
单元数据定位 XL（NDM, NEN) 单元内节点坐标数组 UL (NDM, NEN) 单元内节点位移数组 LD(NDF*NEN) 方程编号
单元数组的计算
直接解法
1. 一维变带宽储存 带宽优化 对应关系: 每行(每列)的最大宽度(最大高度)， 对角线元素位置 PROFIL 子程序完成
2 求解 三角分解
Ka=r K=LU
1 0 ... 0
L
L12
1
...
0
. . ... .
Ln1
Ln2 ....
1
U11
L
0
.
0
U12 ... U1n
N i, N i,
x Q31 x Q31
Ni, xQ12 Ni, yQ32
Ni, yQ22
Ni, xQ32
kij
e
t
1 -1
1 -1
k ij
Jdd
kij e
WiWj f i ,i
ij
C----ISOPARAMETRIC ELEMENT STIFFNESS COMPUTATION
DO 100 J = 1, NEL DB11 = D11*SHP(1,J) DB12 = D12*SHP(2,J) DB21 = D12*SHP(1,J) DB22 = D11*SHP(2,J) DB31 = D33*SHP(2,J) DB32 = D33*SHP(1,J)
C C---FOR EACH I NODE COMPUTE S = BT*DB C
C SHP(1,I) = X-direivative fo I –node shape function
C SHP(2,I) = Y-direivative fo I –node shape function
C SHP(3,I) = shape function for I-node
C XJ
= Jacobian array