矩形薄板弯曲偏微分方程的重心插值配点法求解

前面的笔记克罗内克积及二维重心插值微分矩阵已经记录了二维情况下微分矩阵的推导,同时根据笔记Euler梁弯曲问题的重心无网格配点法求解中的步骤,易得矩形薄板弯曲问题的求解程序如下,ocatve及matlab均可运行:

%% 本代码实现边界为CCCC的矩形薄板弯曲问题的求解
clc;clear;

%% 主程序
Nx = 20
Ny = 20

[x,wx] = chebyshevpoints(Nx,1,1); % 0~2
[y,wy] = chebyshevpoints(Ny,1,2); %0~1
Dsx = barylag(x,4,wx); % 求出关于x的四阶微分矩阵
Dx1 = Dsx(:,:,1);
Dx2 = Dsx(:,:,2);
Dx3 = Dsx(:,:,3);
Dx4 = Dsx(:,:,4);

Dsy = barylag(y,4,wy); % 求出关于y的四阶微分矩阵
Dy1 = Dsy(:,:,1);
Dy2 = Dsy(:,:,2);
Dy3 = Dsy(:,:,3);
Dy4 = Dsy(:,:,4);

[Y,X] = meshgrid(y,x) % 这里是先固定x方向 再对y方向离散 与python的np.meshgrid相反 为了与python对应 交换X,Y的位置

X = X(:) % flatten
Y = Y(:)

%% 构造偏微分方程组
%% 这里采用kron运算
Inx = eye(Nx+1);
Iny = eye(Ny+1);
Lx4 = kron(Iny,Dx4);
% Lx4 = kron(Dx4,Iny);
Ly4 = kron(Dy4,Inx);

Lx2y2 = kron(Dy2,Dx2)
Lx1 = kron(Iny,Dx1);
Ly1 = kron(Dy1,Inx);
L = Lx4 + Ly4 + 2*Lx2y2;
NN = size(X,1);
q = -1;
bendingstiff = 1;
F = ones(NN,1)*q/bendingstiff
%% B.C
II = eye(NN);
bcx0 = find(X==0);
bcx1 = find(X==2);
bcy0 = find(Y==0);
bcy1 = find(Y==1);
%% 施加第一类边界条件
L(bcx0,:) = II(bcx0,:);F(bcx0,:) = 0.0;
L(bcx1,:) = II(bcx1,:);F(bcx1,:) = 0.0;
L(bcy0,:) = II(bcy0,:);F(bcy0,:) = 0.0;
L(bcy1,:) = II(bcy1,:);F(bcy1,:) = 0.0;
%%
bcs_x = [bcx0;bcx1];bcs_y = [bcy0;bcy1];
L0 = Lx1(bcs_x,:);
F0 = zeros(size(bcs_x,1),1);
L1 = Ly1(bcs_y,:);F1 = zeros(size(bcs_y,1),1);
L = [L;L0;L1];F=[F;F0;F1];

%% solve
U = L \ F
minu = min(U)
SHAPE = [Nx+1,Ny+1]
XX = reshape(X,SHAPE)
YY = reshape(Y,SHAPE)
UU = reshape(U,SHAPE)
figure;scatter3(X(:),Y(:),U(:));
figure;surf(XX,YY,UU)


效果如下(Octave编辑器中运行):

阅读原文查看效果演示视频

简介:开源CAE技术爱好者,编程爱好者,计算力学,欢迎志同道合的朋友一起交流。欢迎关注微信公众号:挨踢的土木佬
(0)
打赏 喜欢就点个赞支持下吧 喜欢就点个赞支持下吧

声明:本文来自“挨踢的土木佬”,分享链接:https://www.zyxiao.com/p/302194    侵权投诉

网站客服
网站客服
内容投稿 侵权处理
分享本页
返回顶部