MATLAB 等腰直角三角形单元计算：μ=1/4，杨氏模量 E，厚度 t

等腰直角三角形单元的节点坐标为：

$$\begin{aligned}\x_1 &= 0, &y_1 &= 0 \x_2 &= 0, &y_2 &= h \x_3 &= h, &y_3 &= 0\end{aligned}$$

其中$h$为等腰直角三角形的直角边长。

杨氏模量为$E$，厚度为$t$，则等腰直角三角形单元的球形函数矩阵为：

$$[N] = \frac{1}{2S}\begin{bmatrix}y_2 - y_3 & 0 & y_3 - y_1 & 0 & y_1 - y_2 & 0 \0 & x_3 - x_2 & 0 & x_1 - x_3 & 0 & x_2 - x_1 \x_3 - x_2 & y_2 - y_3 & x_1 - x_3 & y_3 - y_1 & x_2 - x_1 & y_1 - y_2\end{bmatrix}$$

其中$S$为等腰直角三角形的面积，$S = \frac{1}{2} h^2$。

应变矩阵为：

$$[B] = \begin{bmatrix}\frac{\partial N_1}{\partial x} & 0 & \frac{\partial N_2}{\partial x} & 0 & \frac{\partial N_3}{\partial x} & 0 \0 & \frac{\partial N_1}{\partial y} & 0 & \frac{\partial N_2}{\partial y} & 0 & \frac{\partial N_3}{\partial y} \ \frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial y} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial y} & \frac{\partial N_3}{\partial x}\end{bmatrix}$$

单元刚度矩阵为：

$$[K] = t [B]^T [D] [B] [N] [N]^T$$

其中$[D]$为弹性矩阵，对于各向同性材料，$[D]$为：

$$[D] = \frac{E}{1-\mu^2}\begin{bmatrix}1 & \mu & 0 \ \mu & 1 & 0 \ 0 & 0 & \frac{1-\mu}{2}\end{bmatrix}$$

MATLAB代码如下：

% 等腰直角三角形单元的节点坐标和材料参数
h = 1;
E = 1e7;
mu = 1/4;
t = 0.1;

% 计算等腰直角三角形的面积和节点坐标
S = 0.5 * h^2;
x1 = 0; y1 = 0;
x2 = 0; y2 = h;
x3 = h; y3 = 0;

% 计算球形函数矩阵
N = 1/(2*S) * [
y2-y3, 0, y3-y1, 0, y1-y2, 0;
0, x3-x2, 0, x1-x3, 0, x2-x1;
x3-x2, y2-y3, x1-x3, y3-y1, x2-x1, y1-y2
];

% 计算应变矩阵
B = [
diff(N(1,:), 'x'), 0, diff(N(2,:), 'x'), 0, diff(N(3,:), 'x'), 0;
0, diff(N(1,:), 'y'), 0, diff(N(2,:), 'y'), 0, diff(N(3,:), 'y');
diff(N(1,:), 'y'), diff(N(1,:), 'x'), diff(N(2,:), 'y'), diff(N(2,:), 'x'), diff(N(3,:), 'y'), diff(N(3,:), 'x')
];

% 计算弹性矩阵
D = E/(1-mu^2) * [
1, mu, 0;
mu, 1, 0;
0, 0, (1-mu)/2
];

% 计算单元刚度矩阵
K = t * B.' * D * B * N * N.';