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)式,加权余量为:
令
得到含有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方法。