有限元二维热传导波前法MATLAB程序.pdf

本帖最后由 kerbcurb 于 2010-11-24 10:57 编辑

波前法是求解线性方程组的一种方法，广泛用于有限元程序。它最初由英国人（？）B.M. Irons 于1970 在“国际工程计算方法杂志”上发表。30 多年来，波前法有了不少变种。本文所用算法，采于法国人Pascal JOLY 所著《Mise en Oeuvre de la Méthode des Eléments Finis》。


http://space.sznews.com/upload/000/431/33/20070515/1179211789577.MP3

:27bb

资料不错，其中这句话不错
“直到2006 年，我教了一遍有限元后，才弄明白（惟教书才是最好的学习）。
其实简单极了：只需将那二阶偏微分方程，乘上一个任意标量函数，然后
在某个有限单元上积分。”

秘密就在于，这些商业有限元软件，在求解线性方程组时，都尽可能地利
用刚度矩阵的稀疏性。波前法就是这样一种充分考虑了刚度矩阵的稀疏性求解
方程的方法。

这个文献的话有意思，复制出来，呵呵
然而我们没有“任意”地、胡乱地选取，而是居心叵测地，让它恰
恰等于有限元的插值函数！而这些插值函数，恰恰只在给定单元内部非零，在
单元外边一律为零！换句话说，插值函数只是些局部函数。

那么，什么是波前矩阵呢？
就是在某一时刻，彼此发生关系的节点的刚度系数组成的矩阵。这个矩阵
是方的，其中含有最多非零元素的那一行就叫波前。

什么叫某一时刻？就是某一单元！，计算有限
元刚度矩阵有两重循环，最外面那层循环，是对单元循环，即从编号为第一的
单元，到编号最大的单元。在波前法中，当循环到某一单元时，在计算该单元
刚度矩阵以后，还要看看哪一个节点可以消去了，也就是消元。

那么，一个节点，什么时候可以消元了？就是与该节点相连的所有单元都
循环到了的时候。它的含义，就是在整个计算过程中，某一时刻，同时存在于
波前矩阵的节点数，其值最大为11。

但1970
年Irons 发表的那篇关于波前法的文章，大概要算有限元文献里被引用次数最
多的一篇了。波前法的提出，突破了制约有限元应用的瓶颈，使得当时内存量
很小的计算机，可以解算规模很大的题目。可以说，没有波前法，就没有有限
元的大发展，也就没有有限元的今天。

多波前法，就是将网格分成许多子区域，在每个子区域上使用
波前法。这样可以将波前矩阵的阶数降低，以突破内存容许的解题上限并提高
效率。每个子区域间数据的依赖关系与传递，是用所谓消去树来管理的。

波前法是求解线性方程组的一种方法，也是其他数值方法的有利帮手

本帖最后由 kerbcurb 于 2010-11-23 19:48 编辑

我把有限元那一部分用软件转换了一下，发在这里，方便大家研究，代码没有仔细检查错误，估计得有一些错误代码，用的时候自己检查
% Purpose:
%Solution of the heat conduction in a thick cyliner with FEM.
% Calling special functions:
%drawFEM.m for drawing finite element mesh
%cgauss.m for the solution of the system of linear equations.
% Author: Yansheng Jiang
% May 20. 2007

clear
Cdiv = 4% input( 'number of divisions in the circunferencial direction =' ).
Tdiv =  2%input ( 'number of divisions in the thickness direction = ' ):
%tic
a = 100.0; %Inner radius, mm
b = 200.0; %Outer radius, mm

NumeroDeElementos = Cdiv * Tdiv;

NumeroDeNos = ( Cdiv )*( Tdiv + 1);

% Coordonnedades
X = zeros( NumeroDeNos,  1 );
Y = zeros( NumeroDeNos,  1 );
dalpha = 2 * pi / Cdiv;
dT = ( b - a ) / Tdiv;
delta_n = Tdiv + 1;
m = delta_n * ( Cdiv - 1 );
for i = 1:Tdiv + 1

r = a + ( i - 1 ) * dT;

n1 = 1;

nfin = m + n1;

X( n1:delta_n:nfin ) = r * cos(0:dalpha:2 * pi - dalpha );

Y( n1:delta_n:nfin ) = r * sin(0:dalpha:2 * pi - dalpha);
end

% Connectividade de elementos
ICO = zeros(NumeroDeElementos, 4);
for i = 1:Tdiv

ne1 = ( i - 1 ) * Cdiv + 1;

n1 = 1;

n2 = n1 + 1;

n3 = n2 + delta_n;

n4 = n3 - 1;

ICO ( ne1, 1:4 ) = [n1 n2 n3 n4];

ne = ne1;

for j = 2:Cdiv - 1

ne = ne + 1;

ICO( ne, 1:4 ) =ICO( ne - 1, 1:4  ) + delta_n;

end


nefin = ne + 1;


ICO( nefin, 1:2 ) = [ICO( ne, 4 ) ICO( ne, 3 )];

ICO( nefin, 3:4 ) = [n2 n1];
end

%Draw finite element mesh

%drawFEM( NumeroDeElementos, [X Y], ICO );
%Material data
kconx =1; % Coefficient of heat conduction in direction x
kcony = kconx; % Coefficient of heat conduction in direction y
%

%lnicializar a matriz de rigidez global e Tersa independent

TotalDof = NumeroDeNos;

K = zeros ( TotalDof, TotalDof );

F = zeros ( TotalDof, 1 );
%fundoes de interpolacao
syms ksi eta;
%eta
%n4 o o n3
%--> ksi
%n1 n2
psi = 0.25 * [( 1 - ksi ) * ( 1 - eta )   ( 1 + ksi ) * ( 1 - eta )   ( 1 + ksi ) * ( 1 + eta )  ( 1 - ksi ) * ( 1 + eta )];
%derivatives of shape functions with reference to natural coordinates

num_nodes_per_elem = 4;

num_gausspoints_per_elem = 4;

dpsidksi = diff ( psi, ksi );

dpsideta = diff ( psi, eta );

GaussPointCor = 0.677350269189626;

GaussWeight= 1;

GaussPoints = [ -GaussPointCor, -GaussPointCor
                          GaussPointCor, -GaussPointCor
                          GaussPointCor, GaussPointCor
                          -GaussPointCor, GaussPointCor];
dpsidksi_num  = zeros( num_gausspoints_per_elem, num_nodes_per_elem );
dpsideta_num = zeros( num_gausspoints_per_elem, num_nodes_per_elem );
for i = 1:num_gausspoints_per_elem

dpsidksi_num( i, : )  = double ( subs ( dpsidksi, { ksi, eta }, { GaussPoints( i, 1 ), GaussPoints( i, 2 ) } ) );

dpsideta_num( i, : ) = double ( subs ( dpsideta, { ksi, eta }, { GaussPoints( i, 1 ), GaussPoints( i, 2 ) } ) );
end


for elemento = 1:NumeroDeElementos
% numeros globais dos nos do elemento

Nos =ICO( elemento, : );

xe = X(Nos);

ye = Y(Nos);


Kelem = zeros( num_nodes_per_elem, num_nodes_per_elem );


for igauss = 1:num_gausspoints_per_elem
%Jacobian

J = [dpsidksi_num( igauss, : ) * xe   dpsidksi_num( igauss, : ) * ye
              dpsideta_num( igauss, : ) * xe  dpsideta_num( igauss, : ) * ye ];

invj = inv( J );
%Derivatives of shape function with reference to global coordinates

dpsidx = invj( 1, : ) * [dpsidksi_num( igauss, : )
                                           dpsideta_num( igauss, : ) ];

dpsidy = invj( 2, : ) * [dpsidksi_num( igauss, : )
                                            dpsideta_num( igauss, : ) ];

detJ = det( J );

Kelem = Kelem + ( kconx * dpsidx' * dpsidx + kcony * dpsidy' * dpsidy ) * detJ * GaussWeight;

end

% Assemblage to the global stiffness matrix

K( Nos, Nos ) = K( Nos, Nos ) + Kelem;

clear J invj detJ xe ye Kelem;
end

% Apply essential boundary conditions using Reddy's algorithm, p. 173
% boudary conditions: temperature at the inner surface of cylinder =rc.

Tin = 1;

% normal thermal flux dT/dn = 0 at the outer surface.
% nodes on the inner surface where temperature Tin is applied,
ebcdof = [1:delta_n:m + 1];
for i = 1:length( ebcdof )

F = F - K( :, ebcdof( i ) ) * Tin;
end

F( ebcdof ) = Tin;
K( ebcdof, : ) = 0;
K( :, ebcdof ) = 0;

for i = 1:length( ebcdof )

K( ebcdof( i ), ebcdof( i ) ) = 1;
end
K
% Solution for nodal temperature

%x = egauss ( K, F ) ;
x =K \ F;

%toc

感谢楼主分享，最近在找这方面的资料！

回复 tingo 的帖子

波前法是有限元与矩阵处理结合起来综合考虑的，实际上其他的算法如果综合考虑分析也可以得到不少优化。

http://netmusic.sdgb.cn/UpLoadMtv/200861061694201.mp3

{:7_1234:}

本帖最后由 kerbcurb 于 2010-12-2 12:21 编辑

程序修正了若干小问题，具体可以到本人空间去看最新版本
本人空间http://www.mwtee.com/space-uid-18294.html

谢谢分享。。。。