平面应力问题7
平面应力问题
平面应力问题的平面应力应变关系
\[\begin{Bmatrix}\sigma_{xx}\\\\\sigma_{yy}\\\\\tau_{xy}\end{Bmatrix}=\frac{E}{1-\gamma^2}\begin{bmatrix}1&\gamma&0\\\\\gamma&1&0\\\\0&0&\frac{(1-\gamma)}{2}\end{bmatrix}\begin{Bmatrix}\varepsilon_{xx}\\\\\varepsilon_{yy}\\\\\gamma_{xz}\end{Bmatrix}
\]
平面应变问题的应力-应变关系
\[\begin{Bmatrix}\sigma_{xx}\\\\\sigma_{yy}\\\\\tau_{xy}\end{Bmatrix}=\frac{E}{(1+v)(1-2v)}\begin{bmatrix}1-v&-v&0\\\\-v&1-v&0\\\\0&0&\frac{(1-2v)}{2}\end{bmatrix}\begin{Bmatrix}\epsilon_{xx}\\\\\epsilon_{yy}\\\\\gamma_{xy}\end{Bmatrix}
\]
小变形下的应力-位移关系
\[\begin{aligned}&\varepsilon_{xx}=\frac{\partial u}{\partial x}\\&\varepsilon_{yy}=\frac{\partial\nu}{\partial y}\\&\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial\nu}{\partial x}\end{aligned}
\]
写成矩阵形式
\[\begin{Bmatrix}\varepsilon_{xx}\\\\\varepsilon_{yy}\\\\\gamma_{xy}\end{Bmatrix}=\begin{bmatrix}\frac{\partial}{\partial x}&0\\\\0&\frac{\partial}{\partial y}\\\\\frac{\partial}{\partial y}&\frac{\partial}{\partial x}\end{bmatrix}\begin{Bmatrix}u\\\\\nu\end{Bmatrix}
\]
简写
\[\{\epsilon\}=[L]U
\]
对于n节点的单元，未知位移的节点近似插值函数为
\[u=N_{1}u_{1}+N_{2}u_{2}+\cdots+N_{n}u_{n}\\\nu=N_{1}\nu_{1}+N_{2}\nu_{2}+\cdots+N_{n}\nu_{n}
\]
矩阵形式
\[\begin{Bmatrix}u\\v\end{Bmatrix}=\begin{bmatrix}N_1&0&|&N_2&0&|&\ldots&|&N_n&0\\0&N_1&|&0&N_2&|&\ldots&|&0&N_n\end{bmatrix}\begin{Bmatrix}u_1\\\nu_1\\\nu_2\\u_2\\\nu_2\\\vdots\\u_n\\\nu_n\end{Bmatrix}
\]
简写
\[\{U\}=[N]\{a\}
\]
用上式中U代替 \(\{\epsilon\}=[L]U\)中的U
\[\{\epsilon\}=[L]U\\
\qquad=[L][N]\{a\}\\
\qquad=[B]\{a\}\\
\]
其中
\[[B]=\begin{bmatrix}\frac{\partial N_1}{\partial x}&0&|&\frac{\partial N_2}{\partial x}&0&|&\ldots&|&\frac{\partial N_n}{\partial x}&0\\\\0&\frac{\partial N_1}{\partial y}&|&0&\frac{\partial N_2}{\partial y}&|&\ldots&|&0&\frac{\partial N_n}{\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}&|&\ldots&|&\frac{\partial N_n}{\partial y}&\frac{\partial N_n}{\partial x}\end{bmatrix}
\]
矩阵[B] 被称为应变矩阵，可以通过插值函数\(N_i(x,y)\)求导得到
为了得到单元单元上的载荷与节点位移之间的关系，我们用虚功原理
\[\int\limits_{V_{e}}\delta\{\epsilon\}^{\mathrm{T}}\{\sigma\} \mathrm{d}V=\int\limits_{V_{e}}\delta\{U\}^{\mathrm{T}}\{b\} \mathrm{d}V+\int\limits_{\Gamma_{e}}\delta\{U\}^{\mathrm{T}}\{t\} \mathrm{d}\Gamma+\sum\limits_{i}\delta\{U\}_{(\{x\}=\{\overline{x}\})}^{\mathrm{T}}\{P\}_{i}
\]
式中：
\(\{\epsilon\}\)为应变矢量；\(\{\sigma\}\)为应力矢量；\(\{U\}\)为位移矢量
\(\{b\}\)为体力矢量；\(\{t\}\)为表面力矢量
\(\{P\}_i\)为作用在 \(\{x\}=\{\overline{x}\}\) 处的集中力矢量
d\(V\)为单元体积；d\(\Gamma\) 为表面力\(\{t\}\)所作用单元的单元边界
应变的变分 δ{ε}和位移的变分 δ{U} 可以分别表示为
\[\delta\{\epsilon\}=\delta([B]\{a\})=[B]\delta\{a\}\\
\delta\{U\}=\delta([N]\{a\})=[N]\delta\{a\}
\]
应力应变关系
\[\{\sigma\}=[D]\{\epsilon\}=[D][B]\{a\}
\]
代入虚功原理方程
\[\int\limits_{V_{e}}\delta\{a\}^{\mathrm{T}}[B]^{\mathrm{T}}[D][B]\{a\} \mathrm{d}V=\int\limits_{V_{e}}\delta\{a\}^{\mathrm{T}}[N]^{\mathrm{T}}\{b\} \mathrm{d}V+\int\limits_{\Gamma_{e}}\delta\{a\}^{\mathrm{T}}[N]^{\mathrm{T}}\{t\} \mathrm{d}\Gamma+\sum\limits_{i}\delta\{a\}^{\mathrm{T}}[N_{((x)=(\bar{x}))}]^{\mathrm{T}}\{P\}_{i}
\\
原方程\\
\int\limits_{V_{e}}\delta\{\epsilon\}^{\mathrm{T}}\{\sigma\} \mathrm{d}V=\int\limits_{V_{e}}\delta\{U\}^{\mathrm{T}}\{b\} \mathrm{d}V+\int\limits_{\Gamma_{e}}\delta\{U\}^{\mathrm{T}}\{t\} \mathrm{d}\Gamma+\sum\limits_{i}\delta\{U\}_{(\{x\}=\{\overline{x}\})}^{\mathrm{T}}\{P\}_{i}
\]
注意，对于平面单元，单元体积 d\(V\)和单元边界 dГ 可以分别写成 d\(\nu=t\)d\(A\)和 d\(\Gamma=t\)d\(l\) ,式中\(t\)
为单元的厚度，d\(A\)为单元的无限小面积，d\(l\)为单元的无限小边界长度。
由于\(\delta\{a\}\) 是节点位移的变分，因此关于坐标独立，可以把它从积分符号中拿出来并消除，方程就变为
\[\left[\int\limits_{\mathcal{A}_{c}}[B]^{\mathrm{T}}[D][B]t\mathrm{d}A\right]\{a\}=\int\limits_{\mathcal{A}_{c}}[N]^{\mathrm{T}}\{b\}t\mathrm{d}A+\int\limits_{\mathcal{L}_{c}}[N]^{\mathrm{T}}\{t\}t\mathrm{d}l+\sum\limits_{i}[N_{(\{x\}=\{\overline{x}\})}]^{\mathrm{T}}\{P\}_{i}
\]
写成矩阵形式
\[[K_{e}]\{a\}=f_{e} \\
\]
\[[K_e]=\left[\int\limits_{A_e}[B]^\mathrm{T}[D][B]t\mathrm{d}A\right]
\]
上式为刚度矩阵，其中
\[\{f_{e}\}=\int_{A_{e}}[N]^{\mathrm{T}}\{b\}t\mathrm{d}A+\int_{L_{e}}[N]^{\mathrm{T}}\{t\}t\mathrm{d}l+\sum_{i}[N_{((x)=\{\overline{x}\})}]^{\mathrm{T}}\{P\}_{i}
\]
上式为单位力矢量
空间离散
常应变三角形(CST)单元
单元的插值函数
\[\begin{aligned}
N_{1}(x,y)& =m_{11}+m_{12}x+m_{13}y \\
N_{2}(x,y)& =m_{21}+m_{22}x+m_{23}y \\
N_{3}(x,y)& =m_{31}+m_{32}x+m_{33}y
\end{aligned}
\]
其中
\[\begin{gathered}
m_{11} =\frac{x_{2}y_{3}-x_{3}y_{2}}{2A} m_{12} =\frac{y_{2}-y_{3}}{2A} m_{13} =\frac{x_{3}-x_{2}}{2A} \\
m_{21} =\frac{x_3y_1-x_1y_3}{2A} \text{m22} =\frac{y_3-y_1}{2A} \text{m23} =\frac{x_1-x_3}{2A} \\
m_{31} =\frac{x_{1}y_{2}-x_{2}y_{1}}{2A} m_{32} =\frac{y_{1}-y_{2}}{2A} m_{33} =\frac{x_2-x_1}{2A}
\end{gathered}
\]
\[A=\frac{1}{2}\text{det}\begin{bmatrix}1&x_1&y_1\\1&x_2&y_2\\1&x_3&y_3\end{bmatrix}
\]
位移场
三角形单元的位移场可以表示为
\[u=N_{1}u_{1}+N_{2}u_{2}+N_{3}u_{3}\\\nu=N_{1}\nu_{1}+N_{2}\nu_{2}+N_{3}\nu_{3}
\]
矩阵表示
\[\begin{Bmatrix}u\\\nu\end{Bmatrix}=\begin{bmatrix}N_1&0&|&N_2&0&|&N_3&0\\0&N_1&|&0&N_2&|&0&N_3\end{bmatrix}\begin{Bmatrix}u_1\\\nu_1\\u_2\\\nu_2\\u_3\\\nu_3\end{Bmatrix}
\]
简写为
\[\{U\}=[N]\{a\}
\]
应变矩阵
\[\{\epsilon\}=[B]\{a\}
\]
其中
\[\begin{gathered}[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}\end{gathered}
\]
代入得
\[[B]=\begin{bmatrix}m_{12}&0&|&m_{22}&0&|&m_{32}&0\\0&m_{13}&|&0&m_{23}&|&0&m_{33}\\m_{13}&m_{12}&|&m_{23}&m_{22}&|&m_{33}&m_{32}\end{bmatrix}
\]
备注：矩阵[B]与笛卡儿坐标系的\(x\)和 y 坐标无关，它只是节点坐标的函数，并且在整个单元内为常量，因此在整个单元内应变矢量也是常量。这就是为什么这种单元被称为“常应变三角形单元”的原因，本书中将常应变三角形单元简称为 CST(Constant Strain Triangle)单元。
由于矩阵[B]和[D]都是常数矩阵,因此刚度矩阵可表示为
\[[K_e]=[B]^\mathrm{T}[D][B]tA_e
\]
其中 \(A_e\)为单元的面积
单位力矢量
体力
考虑体力{b}是重力引起的，因此
\[\int\limits_{A_{\epsilon}}[N]^{\mathrm{T}}\{b\}t \mathrm{d}A=t\int\limits_{A_{\epsilon}}\begin{bmatrix}N_{1}&0\\0&N_{1}\\N_{2}&0\\0&N_{2}\\N_{3}&0\\0&N_{3}\end{bmatrix}\begin{Bmatrix}0\\0\\-\rho g\end{Bmatrix}\mathrm{d}A=t\begin{bmatrix}0\\-\int\limits_{A_{\epsilon}}N_{1}\rho g \mathrm{d}A\\0\\-\int\limits_{A_{\epsilon}}N_{2}\rho g \mathrm{d}A\\0\\-\int\limits_{A_{\epsilon}}N_{3}\rho g \mathrm{d}A\end{bmatrix}
\]
这里采用式(8.18)和式(8.19)给出的三角形积分公式，计算在整个单元上的插值函数的积分。运用这些公式，可将积分结果整理如下：
\[\int_{A_{e}}N_{1}\rho g \mathrm{d}A=\rho g\int_{A_{e}}N_{1}^{1}N_{2}^{0}N_{3}^{0} \mathrm{d}A=\rho g\frac{1!0!0!}{(1+0+0+2)!}2A_{e}=\rho g\frac{A_{e}}{3}\\\int_{A_{e}}N_{2}\rho g \mathrm{d}A=\rho g\int_{A_{e}}N_{1}^{0}N_{2}^{1}N_{3}^{0} \mathrm{d}A=\rho g\frac{0!1!0!}{(0+1+0+2)!}2A_{e}=\rho g\frac{A_{e}}{3}\\\int_{A_{e}}N_{3}\rho g \mathrm{d}A=\rho g\int_{A_{e}}N_{1}^{0}N_{2}^{0}N_{3}^{1} \mathrm{d}A=\rho g\frac{0!0!1!}{(0+0+1+2)!}2A_{e}=\rho g\frac{A_{e}}{3}
\]
代入得
\[\int_{A_e}[N]^\mathrm{T}\{b\}t \mathrm{d}A=-\frac t3\begin{Bmatrix}0\\\rho gA_e\\0\\\rho gA_e\\0\\\rho gA_e\end{Bmatrix}
\]
可以发现，单元自重在各节点是平均分配的。
表面力
如图 所示的单元作用有大小为\(q\)的均布载荷，该均布载荷作用在单元的边 2-3 上，并
且与 \(x\) 轴的夹角为\(\theta\) 。因此，表面力矢量可以表示为\(\{t\}=\{-q\cos\theta,-q\sin\theta\}^\mathrm{T}\circ\)
则表面力可以表示为
\[\int\limits_{L_{n}}[N]^{\mathrm{T}}\{t\}t \mathrm{d}l=\int\limits_{L_{23}}\begin{bmatrix}0&0\\0&0\\N_2&0\\0&N_2\\N_3&0\\0&N_3\end{bmatrix}\begin{Bmatrix}-q\cos\theta\\-q\sin\theta\\-q\sin\theta\end{Bmatrix}t \mathrm{d}l=t\int\limits_{L_{23}}\begin{Bmatrix}0\\0\\-N_2q\cos\theta\\\\-N_2q\sin\theta\\\\-N_3q\cos\theta\\\\-N_3q\sin\theta\end{Bmatrix}\mathrm{d}l
\]
注意,在边 2-3 上 \(N_1=0\)
采用式(8.18)给出的沿三角形边长的积分公式来计算沿长度方向的积分。运用这个公式上式可以表示为
\[\int\limits_{L_{e}}[N]^{\mathrm{T}}\{t\}t \mathrm{d}l=t\begin{Bmatrix}0\\0\\-q\cos\theta L_{2-3}/2\\-q\sin\theta L_{2-3}/2\\-q\cos\theta L_{2-3}/2\\-q\sin\theta L_{2-3}/2\end{Bmatrix}
\]
从上式可以看出,节点2和节点3平均分配了作用在它们之间的均布载荷\(qL_{2-3}\)。
集中力
% 文件：CST_COARSE_MESH_DATA.m
%
% 下列变量被声明为全局变量，以便被构成程序的所有函数（M文件）使用
%
global nnd nel nne nodof eldof n % 全局变量声明
global geom connec dee nf Nodal_loads % 包括节点数、元素数、节点每元素数、每节点自由度数等
format short e % 设置MATLAB的输出格式为科学计数法
% 定义模型参数
nnd = 21; % 节点数
nel = 24; % 元素数
nne = 3; % 每个元素的节点数
nodof = 2; % 每个节点的自由度数
eldof = nne * nodof; % 每个元素的自由度数
% 节点坐标x和y
geom = zeros(nnd, 3); % 初始化几何矩阵
% ... 此处省略了具体的节点坐标赋值 ...
% 元素连接信息
connec = zeros(nel, 3); % 初始化连接矩阵
# ... 此处省略了具体的元素连接信息 ...
% 材料属性
E = 200000.; % 弹性模量，单位MPa
vu = 0.3; % 泊松比
thick = 5.; % 梁厚度，单位mm
% 形成平面应力的弹性矩阵
dee = formdsig(E, vu); % 调用formdsig函数，该函数未在代码中给出
% 边界条件
nf = ones(nnd, nodof); % 初始化自由度矩阵
% ... 此处省略了具体的边界条件赋值 ...
% 计算自由度的总数
n = 0;
for i = 1:nnd
for j = 1:nodof
if nf(i, j) ~= 0
n = n + 1;
nf(i, j) = n;
end
end
end
% 载荷
Nodal_loads = zeros(nnd, 2); % 初始化节点载荷矩阵
% ... 此处省略了具体的载荷赋值 ...
%%%%%%%%%%%% 输入结束 %%%%%%%%%%%%%
浙公网安备 33010602011771号