计算科学 - 在 matlab pde 工具箱中提取 FEM 矩阵 - 吾爱随笔录

我正在尝试遵循 Matlab 中的动态线性弹性，链接在这里。例如，我的目标是使用 matlab 中的函数 assembleFEMatrices 提取有限元矩阵，并通过 Backward Euler 求解二阶 ODE 的结果系统。

我没有成功将solvepde的结果与我的代码的结果相协调。

让 ''model'' 是根据 Matlab 示例指定的 PDE 模型。

为了提取有限元矩阵，我这样做：

FEM = assembleFEMatrices(model);
M = FEM.M;
K = FEM.K;
F = FEM.F;
G = FEM.G;

这是遵循 Matlab关于 FEM、编写 PDE和提取 FE 矩阵的文档。

让 $x,y$ 成为位移 $x,y$ 方向，分别。得到的 ODE 系统应该是

$M \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + K \begin{bmatrix} x \\ y \end{bmatrix} = G+F$ .

这里， $M$ 是质量矩阵， $K$ 是刚度矩阵， $G$ 来自 Neumann BC 而 $F$ 是强制项。

然后我把上面的系统写成

$\begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + M^{-1}K \begin{bmatrix} x \\ y \end{bmatrix} = M^{-1}(G+F)$ .

然后为了将其简化为一阶 ODE，我引入 $z_1 = \begin{bmatrix} x \\ y \end{bmatrix}$ 和 $z_2 = \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix}$

获得

$\frac{d}{dt} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} 0 & I \\ -M^{-1}K & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} + \begin{bmatrix} 0 \\ M^{-1}(G+F) \end{bmatrix}$ .

然后我使用反向欧拉解决上述系统并进行比较 $z_1$ 结果我从 Matlab 的 solvepde 得到，但它们不匹配。我不认为问题出在我的解决方法上，因为我尝试了 Matlab 的 lsim，得到的结果几乎与后向欧拉相同。我是否错误地提取和误解了 FE 矩阵？

建议表示赞赏。

编辑：下面的代码。我稍微更改了 Matlab 示例中的一些常量，以使代码运行得更快。我还考虑了全光束而不是半光束。在下面的代码中，mySoln 是我从 assembleFEMatrices 获得的解决方案，并将二阶 ODE 转换为一阶 ODE，如上所述。results.NodalSolution 是从 Matlab 的 solvepde 得到的解。我无法让这与 mySoln 协调一致。

% Create PDE model

N = 2; % Two PDE in plane stress elasticity
model = createpde(N);
blength = 1.5; % Beam length, in.
height = .01; % Thickness of the beam, in.

l2 = blength/2;
h2 = height/2;
rect = [3 4 0 l2 l2 0 -h2 -h2  h2 h2]';
g = decsg(rect,'R1',('R1')');
pg = geometryFromEdges(model,g);

figure
pdegplot(g,'EdgeLabels','on');
title('Geometry With Edge Labels Displayed');
axis([-.1 1.1*l2 -5*h2 5*h2]);

E = 3.0e7; % Young's modulus of the material, lbs/in^2
gnu = .3; % Poisson's ratio of the material
rho = .3/386; % Density of the material
G = E/(2.*(1+gnu));
mu = 2*G*gnu/(1-gnu);
c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
f = [0 0]'; % No body forces
specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0, 'f', f);

% BCs
symBC = applyBoundaryCondition(model,'dirichlet','Edge',4,'u',[0 0]);
clampedBC = applyBoundaryCondition(model,'dirichlet','Edge',2,'u',[0 0]);
sigma = 5e1;
presBC = applyBoundaryCondition(model,'neumann','Edge',3,'g',[0 sigma]);
setInitialConditions(model,0,0);
generateMesh(model,'Hmax',height/2,'MesherVersion','R2013a');

% Solve pde using matlab pde toolbox
tfinal = 3e-3;
tlist = linspace(0,tfinal,100);
result = solvepde(model,tlist);

%%
% plot appearance of beam

figure;

X_nodes = result.Mesh.Nodes(1,:);
X_nodes = X_nodes(:);
Y_nodes = result.Mesh.Nodes(2,:);
Y_nodes = Y_nodes(:);

for k = 1:size(result.SolutionTimes,2)

    % Extract displacement at time k
    disp_X = result.NodalSolution(:,1,k);
    disp_X = disp_X(:);

    disp_Y = result.NodalSolution(:,2,k);
    disp_Y = disp_Y(:);

    plot(X_nodes + disp_X, Y_nodes + disp_Y,'*')
    pause(0.1)

end

%%

% try to recover solution by assembling the FEM matrix and time-stepping it
FEM = assembleFEMatrices(model);
Mmat = FEM.M;
Kmat = FEM.K;
Fmat = FEM.F;
Gmat = FEM.G;

nNodes = size(model.Mesh.Nodes,2);

Amat = [zeros(2*nNodes,2*nNodes) eye(2*nNodes); -Mmat\Kmat zeros(2*nNodes,2*nNodes)];
Bmat = [zeros(2*nNodes,1); Mmat\(Fmat + Gmat)];

% Do backward euler here

deltaT = tlist(2) - tlist(1);
IminushA = sparse(eye(4*nNodes)) - deltaT*Amat;

mySoln = zeros(4*nNodes,length(tlist));

for k = 2:length(tlist)
   fprintf('%d...\n',k)
   mySoln(:,k) = IminushA\(mySoln(:,k-1) + deltaT*Bmat);
end