结构有限元刚度方程求解基础Word文档格式.docx
- 文档编号:17613291
- 上传时间:2022-12-07
- 格式:DOCX
- 页数:20
- 大小:291.75KB
结构有限元刚度方程求解基础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},储存顺序如图
,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 结构 有限元 刚度 方程 求解 基础