结构有限元刚度方程求解基础Word文档格式.docx
《结构有限元刚度方程求解基础Word文档格式.docx》由会员分享，可在线阅读，更多相关《结构有限元刚度方程求解基础Word文档格式.docx（20页珍藏版）》请在冰豆网上搜索。
fori=n-1,1,-1
b（i）=（b（i）-U（i,i+1~n）b（i+1~n））/U（i,i）
向后消去法的算法实现（Fortran）:
注:
数组du（i）表示向量
whlf（i）表示向量
二维数组whlk（i,j）表示矩阵
变量neqns表示方程数,临时变量temp
du（neqns）=whlf（neqns）/whlk（neqns,neqns）
do100i=neqns-1,1,-1
temp=0.
do200j=i+1,neqns
temp=temp+whlk（i,j）*du（j）
200enddo
du（i）=（whlf（i）-temp）/whlk（i,i）
100enddo
基于列的形式：
交换循环顺序可得到以上算法的列形式，考虑向前消去，一旦x1解出来，该变量可以从第2~n个方程中去掉，我们可只考虑缩小后的方程组
然后我们算出x2，并且从第3~n个方程中去掉x2，依次类推，例如：
我们有x1=3，那么我们处理2x2方程组
算法1.3（向前消去：
列形式）Lx=b的解覆盖b
forj=1~n-1
b（j）=b（j）/L（j,j）
b（j+1~n）=b（j+1~n）-b（j）L（j+1~n,j）
算法1.4（向后消去：
列形式）Ux=b的解覆盖b
forj=n,2,-1
b（j）=b（j）/U（j,j）
b（1~j-1）=b（1~j-1）-b（j）U（1~j-1,j）
b
（1）=b
（1）/U（1,1）
2>
LU分解
如前面所见，三角方程比较容易求解，那么我们可以把刚度阵[K]经过适当的线性变换等价成一个三角方程组—>
Gauss消去法
例：
在方程组
中，将第一个方程乘以2，并且在第二个方程中减去它则：
矩阵形式：
这就是n=2的Gauss消去法
此算法线性变换过程也可写成矩阵形式，则最终求出一个下三角矩阵L和一个上三角矩阵U，使得A=LU，
然后原始问题Ax=b可通过两个三角求解过程求得：
Ly=b,Ux=y=>
Ax=LUx=Ly=b
在n=2的水平上，如果x1≠0且τ=x2/x1,那么
更一般的，设
且
，如果
定义：
则
形如Mk的矩阵的前k个分量都为0时，一般称高斯变换，这样的矩阵是下三角的，元素
称为乘子，向量
称为高斯向量。
上三角化，设
，通常可找到高斯变换
使得，
是上三角矩阵
如果：
则：
同样地
到n-1步后就可以完成上三角化，
注意，分解过程中需要对角线元素不为“0”，有限元刚度阵满足这一条件，如果刚度阵对角线有“0”元素，则说明结构存在刚体位移，无法求解。
算法2.1（高斯消去）U存于A的上三角部分
fork=1~n-1
rows=k+1~n
A（rows,k）=A（rows,k）/A（k,k）
A（rows,rows~n）=A（rows,rows~n）-A（rows,k）A（k,rows~n）
代码实现：
增广矩阵[Kf]的高斯消去
do100k=1,neqns-1
if（whlk（k,k）==0）exit
do200rows=k+1,neqns
v_gauss=whlk（rows,k）/whlk（k,k）
do300j=1,neqns
whlk（rows,j）=whlk（rows,j）-v_gauss*whlk（k,j）
300enddo
whlf（rows）=whlf（rows）-v_gauss*whlf（k）
其他形式：
do400rows=k+1，neqns
whlk（rows,k）=whlk（rows,k）/whlk（k,k）!
先求高斯向量
400enddo
do200i=k+1,neqns
whlk（i,j）=whlk（i,j）-whlk（i,k）*whlk（k,j）
whlf（i）=whlf（i）-whlk（i,k）*whlf（k）
使用后原矩阵被覆盖
[U]存于[K]的上三角.
3>
特殊的LU分解
当问题中出现如对称性,稀疏性等特性时,需要改进一般矩阵算法以提高效率.
1）
分解
可将一般矩阵A分解成
其中[L],[M]为下三角矩阵,[D]为对角矩阵.
假设已知,L的前j-1列,D的前j-1个对角元素,M的前j-1行,为求L的第j列L（j+1~n,j）,M的第j行M（j,1~j-1）和D的对角元d（j）,取出方程
第j列等式,即:
A（1~n,j）=Lv
其中
（表示v是[U]中的第j列），现将v的上半部v（1~j）定义为已知的下三角方程组:
的解,一旦求得v,则可以计算:
d（j）=v（j）
M（j,i）=v（i）/d（i）,i=1~j-1
v的下半部v（j+1~n）有关系式
重新组织后可用于求L的第j列：
算法3.1:
forj=1~n
从L（1~j,1~j）v（1~j）=A（1~j,j）解出v（1~j）
fori=1~j-1
M（j,i）=v（i）/d（i）
L（j+1~n,j）=（A（j+1~n,j）-L（j+1~n,1~j-1）v（1~j-1））/v（j）
2）
分解（LDU）
如果[A]是对称的，则
分解中的[M]=[L].
向量v（1~j）是由
的前j个分量定义的，由于[M]=[L]，得
由L（1~j，1~j）v（1~j）=A（1~j，j）的第j个方程:
L（j，1~j）v（1~j）=A（j，j）
则有关系式
L（j，1~j-1）v（1~j-1）+L（j，j）v（j）=A（j，j）
ð
v（j）=A（j，j）-L（j，1~j-1）*v（1~j-1）
故有:
算法3.2：
v（i）=A（j,j）-L（j,1~j-1）v（1~j-1）
v（j）=A（j,j）-L（j,1~j-1）v（1~j-1）
代码实现，[K]被下三角L，对角D覆盖
programmain
implicitnone
parameterneqns=3
real:
:
whlk（neqns,neqns）=（/10,20,30,20,45,80,30,80,171/）
v（neqns）=0，aa
integer:
i,j,k,l
!
j=1时单独计算，否则出错
aa=1.0/whlk（1,1）
do1000i=2,neqns
whlk（i,1）=whlk（i,1）*aa
1000enddo
do100j=2,neqns
计算v（1~j-1）
do200i=1,j-1
v（i）=whlk（j,i）*whlk（i,i）
计算v（j）
do300k=1,j-1
v（j）=whlk（j,j）-whlk（j,k）*v（k）
！
储存d（j）
whlk（j,j）=v（j）
计算L（j+1~n,j）
do400k=j+1,neqns
do500l=1,j-1
whlk（k,j）=（whlk（k,j）-whlk（k,l）*v（l））/v（j）
500enddo
执行完后
分解后可用Ly=b（向前消去），Dz=y，LTx=Z（向后消去）得到Ax=b的解.
改进，不用v数组,200和300循环可以合并，500循环中的V也可以用上面计算的刚度阵中的元素相乘得到。
顺便调整下循环编号。
programmain！
codeLDU
aa,whlf（neqns）=（/1,2,3/）,u（neqns）
do100i=2,neqns
do200k=1,i-1
whlk（i,i）=whlk（i,i）-whlk（i,k）*whlk（i,k）*whlk（k,k）
do400j=i+1,neqns
do500k=1,i-1
whlk（j,i）=（whlk（j,i）-whlk（j,k）*whlk（i,k）*whlk（k,k））
whlk（j,i）=whlk（j,i）/whlk（i,i）
write（*,110）（（whlk（i,j）,j=1,3）,i=1,3）
110format（3e12.4）
u=whlf把whlf赋值给u
do2000i=1,neqns！
下三角向前消去解Ly=whlf
do2100k=1,i-1
u（i）=u（i）-whlk（i,k）*u（k）
2100enddo
2000enddo
do2200i=1,neqns!
解Dz=y
u（i）=u（i）/whlk（i,i）
2200enddo
do2300i=neqns-1,1,-1!
解上三角U=LTx=z
do2400j=i+1,n
u（i）=u（i）-whlk（j,i）*u（j）
2400enddo
2300enddo
P．S．
以上代码只用到了结构刚度阵的对称特性，关于利用刚度阵对称正定的特性，则存在Cholesky分解。
如果A是对称正定的方阵，则存在唯一一个对角元全部大于0的下三角矩阵G，满足
比较等式
的第j列可得：
也就是说：
如果G的前j-1列已知，则可以计算出v，由
，L是单位上三角，D是对角阵，所以
得出：
（Cholesky）算法3.3
v（j~n）=A（j~n,j）
fork=1~j-1
v（j~n）=v（j~n）-G（j,k）G（j~n,k）
G（j~n,j）=v（j~n）/sqrt（v（j））
parameterneqns=2
whlk（neqns,neqns）=（/2,-2,-2,5/）
v（neqns）=0
do100j=1,neqns
do200i=j,neqns
v（i）=whlk（i,j）
v（i）=v（i）-whlk（j,k）*whlk（i,k）
whlk（i,j）=v（i）/sqrt（v（j））
100enddo
或：
v（neqns）=0,temp
do100i=1,neqns
if（i>
1）then
do200j=i,neqns
do300k=1,i-1
whlk（j,i）=whlk（j,i）-whlk（j,k）*whlk（i,k）
endif
temp=sqrt（whlk（i,i））
do400l=i,neqns
whlk（l.i）=whlk（l,i）/temp
计算得出：
，G覆盖A的下三角。
Cholesky分解和
分解所需的计算量差不多，但是后面只需解
两个三角方程组。
4>
skyline方法
设有如下刚度阵
把每行（或列）第一个非0元素的列号（或行号）记在一个数组中，即index1（i）={1，1，2，2，1，4，3，1}，这组数字，有个特点：
即使对[K]进行LDU分解,
index1（i）对[L],[U]矩阵不变化.
例:
……
可以看到[L]的行或[U]的列中,第一个非0元素的列或行号仍然是{1,1,2,2,1,4,3,1}.而且分解所产生的非0元素如k25,k35,k45,k38,k48,k58,k68都在每行或列第一个非0元素和对角线之间.
利用以上特性,我们可以把二维的刚度阵用一个一维数组存储.
过程如下
我们把按照index1把非0部分的轮廓线划出,这个轮廓线一般被称为skyline,skyline和对角线所夹部分称为profile。
由图可见，skyline方法在非0元素聚集在对角线附近时，计算效率最高。
一维数组储存顺序如下：
由于是对称阵，故只存储上三角或下三角。
为了记住原矩阵中对角线上的元素在一维刚度阵中的位置，需要如下数组：
index2（i）={1，3，5，8，13，16，21，29}
第j列中，从第一个非0元素到对角线元素，一共有j-index1（j）+1个元素，那么已知第j列对角线的位置，下一个对角线的位置为：
其中j-index1（j）+1被称为列高。
原矩阵中Kij在一维刚度阵的位置如图：
综合以上，对照“！
codeLDU”将skyline方法用于改进
分解，
whlk2（neqns,neqns）=（/10,20,30,20,45,80,30,80,171/）
whlk（neqns*neqns）=0.！
方程数多则动态分配
aa
index1（neqns）,index2（neqns）
i,i1,i2,j,j1,j2,k,k1,k2,k3,iii,ncount
do100i=1,neqns!
生成index1
ncount=1
do200j=1,i
if（whlk2（i,j）==0.）then
ncount=ncount+1
else
index1（i）=ncount
exit
index2
（1）=1！
生成index2
do300i=2,neqns
index2（i）=index2（i-1）+i-index1（i）+1
do400i=1,neqns！
将2维刚度阵转为1维
do500j=i,neqns
k=index2（j）-j+i
whlk（k）=whlk（k）+whlk2（i,j）
aa=1.0/whlk
（1）！
修改LDU分解
if（index1（i）==1）then
i1=index2（i-1）+1
whlk（i1）=whlk（i1）*aa
do1100i=2,neqns
i1=index1（i）
i2=index2（i）
do1200k=i1,i-1
k1=index2（i）-i+k
k2=index2（k）
whlk（i2）=whlk（i2）-whlk（k1）*whlk（k1）*whlk（k2）
1200enddo
do1300j=i+1,neqns
j1=index1（j）
if（j1<
=i）then
j2=index2（j）-j+i
if（i1>
j1）iii=i1
if（j1>
i1）iii=j1
do1400k=iii,i-1
k2=index2（j）-j+k
k3=index2（k）
whlk（j2）=whlk（j2）-whlk（k1）*whlk（k2）*whlk（k3）
1400enddo
whlk（j2）=whlk（j2）/whlk（i2）
1300enddo
1100enddo
分解后得到whlk（i）={10，2，5，3，4，1}，储存顺序如图
，
升级会员 