以下是一个六节点三角形单元有限元matlab求解含4个单元的矩型薄板变形的程序示例:
```matlab
clear;clc;close all;
% 定义常数
E = 2.1e11; % 弹性模量
nu = 0.3; % 泊松比
h = 0.01; % 板厚
q = 5000; % 载荷
% 定义节点坐标
x(1) = 0; y(1) = 0;
x(2) = 0.5; y(2) = 0;
x(3) = 1; y(3) = 0;
x(4) = 0.25; y(4) = 0.5;
x(5) = 0.75; y(5) = 0.5;
x(6) = 0.5; y(6) = 1;
% 定义单元节点编号
IEN = [1 4 2 5 6 3; 4 2 5 1 3 6;
2 4 5 1 6 3; 4 5 3 2 6 1];
% 定义初始位移和载荷向量
U = zeros(2*length(x),1);
F = zeros(2*length(x),1);
% 定义刚度矩阵和载荷向量
K = zeros(2*length(x));
for e = 1:4
Ke = zeros(12);
xe = x(IEN(e,:));
ye = y(IEN(e,:));
for i = 1:6
B = [ye(mod(i,6)+1)-ye(mod(i-2,6)+1), 0;
0, xe(mod(i-2,6)+1)-xe(mod(i,6)+1);
xe(mod(i-2,6)+1)-xe(mod(i,6)+1), ye(mod(i,6)+1)-ye(mod(i-2,6)+1)];
detJ = det(B)/2;
invB = inv(B);
D = E/(1-nu^2)*[1, nu, 0;
nu, 1, 0;
0, 0, (1-nu)/2];
Be = [invB(1,1), 0, invB(2,1), 0, invB(3,1), 0;
0, invB(2,2), 0, invB(3,2), 0, invB(1,2);
invB(2,2), invB(1,1), invB(3,2), invB(2,1), invB(1,2), invB(3,1)];
Ke(2*i-1:2*i,2*i-1:2*i) = Ke(2*i-1:2*i,2*i-1:2*i) + Be'*D*Be*detJ*h;
end
K([2*IEN(e,:)-1; 2*IEN(e,:)], [2*IEN(e,:)-1; 2*IEN(e,:)]) = K([2*IEN(e,:)-1; 2*IEN(e,:)], [2*IEN(e,:)-1; 2*IEN(e,:)]) + Ke;
end
F(2*6-1) = -q*h*0.5;
% 处理边界条件
K([1:2:2*length(x)],:) = 0;
K([1:2:2*length(x)], [1:2:2*length(x)]) = eye(length(x));
F(1:2:2*length(x)) = 0;
% 解方程
U = K\F;
% 输出结果
disp('节点编号 x位移 y位移');
for i = 1:length(x)
disp([i, U(2*i-1), U(2*i)]);
end
```
该程序采用线性三角形单元进行离散化,节点坐标和单元节点编号需要根据具体问题进行修改。在4个单元的情况下,需要将单元节点编号矩阵IEN的行数修改为4。