有限元| 刚架临界荷载的新算法

fem178

发布于 2024-05-31 20:39:39

1140

发布于 2024-05-31 20:39:39

文章被收录于专栏：数值分析与有限元编程数值分析与有限元编程

梁单元新的弹性刚度矩阵和几何刚度矩阵

▲图1

图1所示为杆件单元上的任一微段

是其失稳之前的位置,

A'B'

为轴压力达到临界值时可能出现的分支平衡位置。设微段由平衡位置

A'B'

发生无限小的横向虚位移至

A''B''

。弹性稳定分析是二阶分析，虚位移应该在变形后。

ds1

和

ds2

分别为虚位移发生前、后微段的长度。微段的轴向虚应变可表示为

\delta \epsilon = \frac{ds2-ds1}{ds1} \quad (1)

由弧长公式

\begin{split} ds1 &= \sqrt{1+(\frac{dv}{dx})^2}dx\approx dx+\frac{1}{2}(\frac{dv}{dx})^2dx\\ ds2 &= \sqrt{1+[\frac{d(v+\delta v)}{dx}]^2}dx\approx dx+\frac{1}{2}(\frac{dv}{dx}+\frac{d\delta v}{dx})^2dx\\ \end{split} \quad (2)

(2)代入(1)，并略去高阶微量,得由横向虚位移产生的轴向虚应变为

\delta \epsilon = \frac{dv}{dx} \frac{d\delta v}{dx} \quad (3)

在梁单元的4个节点自由度中，

\omega

和

\theta

的单位不同，现在用

\hat{\theta} = l\theta \quad (4)

将单位统一起来。其中

是单元长度。梁的挠度

\begin{split} v &= N_1v_1 + N_2\hat{\theta_1} + N_3v_2 + N_4\hat{\theta_2} \\ &= \mathbf N \hat{\mathbf d}\\ \end{split} \quad (5)

其中

\begin{cases} N_1 = 1-3\xi^2 + 2\xi^3 \\ N_2 = -\xi + 2\xi^2 -\xi^3 \\ N_3 = 3\xi^2 - 2\xi^3 \\ N_4 = \xi^2 - \xi^3 \\ \end{cases}

\hat{\mathbf d} = \begin{Bmatrix} \omega_1 \\ \hat{\theta}_1 \\ \omega_2 \\ \hat{\theta}_2\\ \end{Bmatrix} \quad (6)

于是

\begin{split} \frac{dv}{dx} &= \frac{d \mathbf N }{dx}\hat{\mathbf d} \\ \frac{d\delta v}{dx} &= \frac{d \mathbf N }{dx}\delta \hat{\mathbf d}\\ \end{split} \quad (7)

两种坐标系的映射

\xi = \frac {x}{l} \quad (8)

由(9)(10)得

\begin{split} \frac{dv}{dx} & = \frac{d \mathbf N }{d\xi} \frac{d\xi }{dx}\hat{\mathbf d} \\ & = \frac{1}{l}[-6\xi+6\xi^2 \quad -1+4\xi-3\xi^2 \quad 6\xi-6\xi^2 \quad 2\xi-3\xi^2] \hat{\mathbf d}\\ \end{split} \quad (9)

\begin{split} \frac{d\delta v}{dx} & = \frac{d \mathbf N }{d\xi} \frac{d\xi }{dx}\delta \hat{\mathbf d} \\ & = \frac{1}{l}[-6\xi+6\xi^2 \quad -1+4\xi-3\xi^2 \quad 6\xi-6\xi^2 \quad 2\xi-3\xi^2] \delta \hat{\mathbf d}\\ \end{split} \quad (10)

梁横截面虚应变

\begin{split} \delta \hat{\epsilon} & = \frac {\partial ^2 \delta v}{\partial x^2 }\\ &= \frac {1}{l^2} \frac {\partial ^2 \delta v}{\partial \xi^2 }\\ &= \frac {1}{l^2} \hat{\mathbf B} \delta \hat{\mathbf d}\\ \end{split} \quad (11)

其中

\hat{\mathbf B} = [-6+12\xi \quad 4-6\xi \quad 6-12\xi \quad 2-6\xi]

▲图2

如图2所示，单元轴向荷载所作的虚功为

\delta W_e = F_P \Delta = F_P\int_{x=0}^l\delta \epsilon dx \quad (12)

式中单元的轴向荷载

F_P

以受拉为正。(3)代入(4)可得

\delta W_e = F_P\int_{x=0}^l \frac{dv}{dx} \frac{d\delta v}{dx} dx \quad (13)

内力虚功

\begin{split} \delta W_i & = \int_V \delta \epsilon^T \sigma dV \\ & = \int_V \delta \hat{\mathbf d}^{T} \hat{\mathbf B}^T E \hat{\mathbf B}\hat{\mathbf d} dV \\ & = \delta \hat{\mathbf d}^{T} \int_V \hat{\mathbf B}^T E \hat{\mathbf B} dV \hat{\mathbf d} \end{split} \quad \cdots (14)

由虚功原理

\delta W_e = \delta W_i

可得

\int_V \hat{\mathbf B}^T E \hat{\mathbf B}dV \hat{\mathbf d} = F_P\int_{x=0}^l (\frac{d \mathbf N }{dx})^T \frac{d \mathbf N }{dx} dx \hat{\mathbf d} \quad (15)

记

\mathbf k^e = \int_V \hat{\mathbf B}^T E \hat{\mathbf B} dV

\mathbf k_{\sigma}^e = F_P\int_{x=0}^l (\frac{d \mathbf N }{dx})^T \frac{d \mathbf N }{dx} dx

则有

\mathbf k^e \hat{\mathbf d} = F_P \mathbf k_{\sigma}^e \hat{\mathbf d} \quad (16)

上式是一个求广义特征值问题。

经积分计算后可得

\mathbf k^e = \frac {2EI}{l^3}\begin{bmatrix} 6& 3& -6& -3 \\ -3& 2& 3& 1 \\ -6& 3& 6& 3 \\ -3& 1& 3& 2 \\ \end{bmatrix} \quad (17)

\mathbf k_{\sigma}^e = \frac{F_P}{10l} \begin{bmatrix} 12 & -1& -12 & -1 \\ -1 & \frac{4}{3} & 1 & -\frac{1}{3} \\ -12 & 1 & 12 & 1\\ -1 & -\frac{1}{3} & 1 & \frac{4}{3} \end{bmatrix}\quad (18)

算例

▲图3

已知简支梁的长度为

，划分2个单元，求其临界荷载。

\mathbf k^e = \frac {16EI}{L^3}\begin{bmatrix} 6& 3& -6& -3 \\ -3& 2& 3& 1 \\ -6& 3& 6& 3 \\ -3& 1& 3& 2 \\ \end{bmatrix}

\mathbf k_{\sigma}^e = \frac{1}{5L} \begin{bmatrix} 12 & -1& -12 & -1 \\ -1 & \frac{4}{3} & 1 & -\frac{1}{3} \\ -12 & 1 & 12 & 1\\ -1 & -\frac{1}{3} & 1 & \frac{4}{3} \end{bmatrix}

考虑边界条件之后，整体弹性刚度矩阵

\mathbf K = \frac {16EI}{L^3}\begin{bmatrix} 2& 3& 1& 0\\ 3& 12& 0& -3 \\ 1& 0& 4& 1 \\ 0& -3& 1& 2 \\ \end{bmatrix}

整体几何刚度矩阵

\mathbf K_{\sigma} = \frac{1}{5L} \begin{bmatrix} \frac{4}{3} & 1& -\frac{1}{3} & 0 \\ 1 & 24 & 0& -1 \\ -\frac{1}{3} & 0 & \frac{8}{3} & -\frac{1}{3} \\ 0& -1 & -\frac{1}{3} & \frac{4}{3} \end{bmatrix}

由

\mathbf K \hat{\mathbf d} = F_P \mathbf K_{\sigma} \hat{\mathbf d}

得

\begin{bmatrix} 2& 3& 1& 0\\ 3& 12& 0& -3 \\ 1& 0& 4& 1 \\ 0& -3& 1& 2 \\ \end{bmatrix} \begin{Bmatrix} q_2\\ q_3 \\ q_4 \\ q_6\\ \end{Bmatrix} = \lambda \begin{bmatrix} \frac{4}{3} & 1& -\frac{1}{3} & 0 \\ 1 & 24 & 0& -1 \\ -\frac{1}{3} & 0 & \frac{8}{3} & -\frac{1}{3} \\ 0& -1 & -\frac{1}{3} & \frac{4}{3} \end{bmatrix} \begin{Bmatrix} q_2\\ q_3 \\ q_4 \\ q_6\\ \end{Bmatrix}

其中

\lambda = \frac{80F_PL^2}{EI}

，在MATLAB中求特征值与特征向量的代码如下

A = [2,3,1,0;
     3,12,0,-3;
     1,0,4,1;
     0,-3,1,2];

B = [4/3, 1, -1/3 0;
     1, 24, 0, -1;
     -1/3, 0, 8/3, -1/3;
     0, -1, -1/3, 4/3];

[X,D] = eig(A,B)
% D是对角矩阵，对角线上的元素是特征值。