有限元FEM求解一维电磁场问题 Rits法 Galerkin法

3.FEM

两块无线大的PEC板位于YOZ平面，一块位于x=0电位为0V，另一块位于x=4电位为20V，使用一维有限元方法求解板间电位。

（1）求解电位的解析表达式

由泊松方程

两边积分2次得到

由边界条件，，得到电位的解析表达式为：

（2）Ritz变分法FEM

利用讲义（21）式

将区域离散化后，在每个子区域上的的表达式为：

代入（21）式得到

然后对求偏导数，令其等于零，即

得到含有N-2个未知量的N-2个线性方程组，求解线性方程组即得到离散点上的解，然后由的表达式得到整个区域上的解。

编写Matlab程序如下：

clc; clear;

N=20;

xstart=0; xend=4; %x range

xx=linspace(xstart, xend, N);

syms x y F1 F2

%construct y

y=ones(1,N);

y=sym(y);

for i=1:N

y(i)=['y', num2str(i)];

end

y(1)=0; y(N)=20; %bound

f=x^2 + 1/2*x + 1/4;

tic

% Ritz method

for i=1:N-1

Fd(i) = 0.5*int( ( ( y(i+1)-y(i) ) / ( xx(i+1)-xx(i) ) )^2, x, xx(i), xx(i+1) )...

+int( f * ( y(i) * ( xx(i+1)-x ) / ( xx(i+1)-xx(i) ) + y(i+1) * ( x-xx(i) ) / ( xx(i+1)-xx(i) ) ) , x, xx(i), xx(i+1) );

end

F=sum(Fd);

for i=2:N-1

Fdiff(i)=collect(diff(F, y(i)));

for j=2:N-1

temp=coeffs(Fdiff(i), y(j)); % Extract the coefficient matrix A, Ax=b

if( size(temp) == 1 )

A(i-1,j-1) = 0;

else

A(i-1,j-1) = temp(2);

end

temp=coeffs(Fdiff(i)); % Extract the right matrix b

b(i-1)=temp(1);

end

end

A=double(A); b=double(b);

yy(2:N-1)=inv(A)*(-b');

t=toc

yy(1)=double( y(1) );

yy(N)=double( y(N) );

plot( xx, yy, 'b*'); hold on

plot( xx, yy,'b--');

%Analytic solution

ax=xstart:0.001:xend;

ay=1/12*ax.^4 + 1/12*ax.^3 + 1/8*ax.^2 - 13/6*ax;

plot(ax, ay, 'r')

title(['Ritz FEM 电位分布，N=', num2str(N)]);

xlabel('距离'); ylabel('电位');

legend('精确数值解','数值解', '解析解');

N=5,10,20,30时的电位分布分别如下图：

（3）Garlerkin法FEM

利用讲义（29）式，加权余量为：

将区域离散化后，在每个子区域上的的表达式为：

检验函数取为：

将和代入（29）式得到：

令

得到含有N-2个未知量的N-2个线性方程组，求解线性方程组即得到离散点上的解，然后由的表达式得到整个区域上的解。

编写Matlab程序如下：

clc; clear;

N=20;

xstart=0; xend=4; %x range

xx=linspace(xstart, xend, N);

syms x F1 F2

%construct y

y=ones(1,N);

y=sym(y);

for i=1:N

y(i)=['y', num2str(i)];

end

y(1)=0; y(N)=20; %bound

f=x^2 + 1/2*x + 1/4;

tic

% Galerkin method

for i=2:N-1

F(i) = -int( ( ( y(i)-y(i-1) ) / ( xx(i)-xx(i-1) ) ) * 1 / ( xx(i)-xx(i-1) ), x, xx(i-1), xx(i) )...

-int( ( ( y(i+1)-y(i) ) / ( xx(i+1)-xx(i) ) ) * (-1) / ( xx(i+1)-xx(i) ), x, xx(i), xx(i+1) )...

-int( f * ( x-xx(i-1) ) / ( xx(i)-xx(i-1) ) , x, xx(i-1), xx(i) )...

-int( f * ( xx(i+1)-x ) / ( xx(i+1)-xx(i) ) , x, xx(i), xx(i+1) );

end

for i=2:N-1

Ff(i)=collect(F(i));

for j=2:N-1

temp=coeffs(Ff(i), y(j)); % Extract the coefficient matrix A, Ax=b

if( size(temp) == 1 )

A(i-1,j-1) = 0;

else

A(i-1,j-1) = temp(2);

end

temp=coeffs(Ff(i)); % Extract the right matrix b

b(i-1)=temp(1);

end

end

A=double(A); b=double(b);

yy(2:N-1)=inv(A)*(-b');

t=toc

yy(1)=double( y(1) );

yy(N)=double( y(N) );

plot( xx, yy, 'b*', xx, yy, 'b--'); hold on

%Analytic solution

ax=xstart:0.001:xend;

ay=1/12*ax.^4 + 1/12*ax.^3 + 1/8*ax.^2 - 13/6*ax;

plot(ax, ay, 'r')

title(['Galerkin FEM 电位分布，N=', num2str(N)]);

xlabel('距离'); ylabel('电位');

legend('精确数值解','数值解', '解析解');

N=5,10,20,30时的电位分布分别如下图：

（4）比较Ritz和Garlerkin法FEM的收敛性

两种方法均收敛。

收敛速度（使用Matlab的tic，toc计时积分、提取系数矩阵和解方程的时间）：

N=10时，Ritz用时0.266，Garlerkin用时0.297

N=20时，Ritz用时0.953，Garlerkin用时1.078

N=30时，Ritz用时2.047，Garlerkin用时2.219

所以，Ritz方法收敛速度要快于Garlerkin方法。