波前法的基本概念与算法

cem-uestc · 发表于 2010-11-15 18:01:27

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

kerbcurb · 发表于 2010-11-17 15:22:50

本帖最后由 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