当前位置:文档之家› 弹塑性力学讲义程序设计1

弹塑性力学讲义程序设计1


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
线性方程组求解 特点: 稀疏 对称 正定
d j 1 -U a j 1 j d j
U a j 1 .U a j 1 - U a j U a j .U a j
j
2. 预条件共轭梯度法 条件数:最大特征值与最小特征值之比 条件数越大,迭代收敛越困难 预条件目的:减小刚度矩阵的条件数 BKa =Br P = BK c = Br Pa = c
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 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
有限元程序
• 数据输入(前处理) • 求解计算 • 结果整理与输出(后处理)
数据输入(前处理)
• 单元几何: 单元编号; 节点坐标
• 材料性质
• 边界条件(约束)
• 荷载(体积力, 面积力)
求解计算
• 单元刚度矩阵(形函数矩阵,B矩阵,D矩阵,高斯积分)
• 单元荷载列阵
• 组集整体刚度矩阵和整体荷载列阵
• 方程组求解
结果整理(后处理)
• 等值线
• 矢量图
• 变形网格图
编程基本原则
• 结构化,模块化
• 可读性强
变量命名简单直观,适当注释
程序条理清楚,逻辑性强
• 节约内存,计算效率优化
FEAP程序
输入基本变量 D(10,NUMAT) F(NDF,NUMNP) ID(NDF,NUMNP) IE(NUMMAT) IX(NEN1,NUMEL) T(NUMNP) X(NDM, NUMNP) 材料性质数组 节点力数组 约束信息数组 各组材料对应的单元类型 单元编号及材料组号 节点温度 节点坐标
D矩阵
D11 D 12 0 D12 D11 0 0 0 D33
单元刚度矩阵
k ij BiT 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
B矩阵
N i, x Bi 0 N i, y
0
N i, y N i, x
D矩阵
D11 D 12 0 D12 D11 0 0 0 D33
SUBROUTINE SHAPEF(S, T, XL, XSJ, SHP) C--- shape function routine for 4-node isoparametric quadrilateral 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 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 C--- IN NATURAL COORDINATES C
STRE
输出位移
输出应力
TANG
单元数据定位
计算整体刚度矩阵
XL(NDM, NEN) 单元内节点坐标数组 UL (NDM, NEN) LD(NDF*NEN) 单元内节点位移数组 方程编号
单元数组的计算
以平面四边形等参单元为例 形函数
Ni 1 1 i 1 i 4
x N i xi u N i ui
形函数导数
y N i yi v N i vi
N i, x , y , N i, x N i, x , y , N i, y
N i, x 1 y , - y , N i, N i, y J - x , x , N i,
C C---- FOR EACH J NODE COMPUTE DB = D*B C 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---- SG,TG ARE INTERGRATION POINTS IN C---NATRUAL COORDINATES C---- WG IS THE INTERGARION WEIGHT C C---- FOR EACH INTERGRATION COMPUTE 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
1
-1 -1
k Jdd
1 ij
i j i i
k W W f ,
i j
C----ISOPARAMETRIC ELEMENT STIFFNESS COMPUTATION 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 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---- S IS THE ELMENT
SUBROUTINE PGAUSS(L,LINT,R,Z,W) C C---- GAUSS POINTS AND WEIGHTS FOR TWO 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
Q j DB j
N i, x Q11 N i, x Q31 k ij N i, y Q21 N i, x Q31
N i, x Q12 N i, y Q32 N i, y Q22 N i, x Q32
k
ij
e ij
e
t
C C---FOR EACH I NODE COMPUTE S = BT*DB C 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
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
解法: 直接解法 迭代解法
直接解法
1.
一维变带宽储存 带宽优化 对应关系: 每行(每列)的最大宽度(最大高度), 对角线元素位置 PROFIL 子程序完成 2 求解 三角分解
Ka=r K=LU
1 L 12 L . Ln1 0 1 . ... ... ... 0 0 . 1
相关主题