工程数值分析---燕山大学

对线性问题有
通常上述元件方程是在各自不同的局部坐标中建立起来的， 通常上述元件方程是在各自不同的局部坐标中建立起来的，在进行下一步骤 之前要把元件方程转换到总体坐标系中表示。 之前要把元件方程转换到总体坐标系中表示。 3) 把所有元件的方程组装起来得到系统的总体方程
[ K ]{a} = {F }
应该注意，系统总体方程组装的物理意义是节点参数 的连续性和节点的平衡 的连续性和节点的平衡。 应该注意，系统总体方程组装的物理意义是节点参数{a}的连续性和节点的平衡。 一个系统中每个元件的节点数可以不同，但是每个节点的自由度数必须相同。 一个系统中每个元件的节点数可以不同，但是每个节点的自由度数必须相同。
N
i+1 = A E (
u i−1
Li
i
u (x ) ui
X 图 2-4
ui+1 −ui ) Li+1
2-1 有限单元法的概念
有限单元法求解直杆拉伸： 有限单元法求解直杆拉伸：
Ni
i
q (Li + Li+1 ) 2
4、以i结点为对象，列力的平衡方程 结点为对象，
∑Fx = 0
q (Li + Li+1) N i − N i+1 = 2 Li λi = Li+1
每个单元的节点号如下所示： 每个单元的节点号如下所示：
2-2 杆系结构有限单元法
杆的长度变化为（以受拉为正，受压为负） 杆的长度变化为（以受拉为正，受压为负）
在节点 j 处的端点轴向力为
其中
该力在x, 方向的分量表达式为 方向的分量表达式为： 该力在 y方向的分量表达式为：
2-2 杆系结构有限单元法
2-1 有限单元法的概念
有限单元法求解直杆拉伸： 有限单元法求解直杆拉伸： 假设线单元数为3个的情况， 假设线单元数为3个的情况， 平衡方程有3 平衡方程有3个： q 2 i=1时 i=1时， 2 u 1 −u 2 = a
0
L1 = a
1
u0
u1
L2 = a
2
L3 = a
u2
u3
3
图 2-6
EA q 2 − u 1 + 2 u 2 −u 3 = a i=2时 i=2时, EA q 2 − u 2 +u3 = a i=3时 i=3时， 2 EA
2-1 有限单元法的概念
有限单元法求解直杆拉伸： 有限单元法求解直杆拉伸：
1、离散化
2、外载荷集中到结点上，即把投影部 外载荷集中到结点上， 分的重量作用在结点i 分的重量作用在结点i上
1 2
L 1
L2 Li Li+1
i-1
i-1 i i+1
Li
i
Li+1
n-1 n
q (Li + Li+1 ) 2
2-2 杆系结构有限单元法
写成分块形式为
3. 组装各单元的方程得到结构各节点的平衡方程。 组装各单元的方程得到结构各节点的平衡方程。
4. 施加边界条件：{d1}={d2}={0}。施加这种边界条件最简单的方法是从 施加边界条件： 。 总体方程中划去对应行、 得到： 总体方程中划去对应行、列，得到：
2-2 杆系结构有限单元法
节点位移分量记为ui 节点位移分量记为 , vi , uj ,和vj ，对应的端点力为 Xie, Yie, Xje, 和Yje 。 和
2-2 杆系结构有限单元法
单元号
E
节点号
i, j
节点坐标
( xi , yi )；(xj , yj)
弹性模量
Ee
Байду номын сангаас
横截面积
Ae
杆件的长度可由下式计算
上式对任意的 {U }成立，得到 成立， 成立
这就是局部坐标系和总体坐标系中单元刚度矩阵之间的坐标转换关系。 这就是局部坐标系和总体坐标系中单元刚度矩阵之间的坐标转换关系。
2. 边界条件
总体刚度矩阵是奇异，就是结构存在刚体位移 刚体运动 刚体运动)。 总体刚度矩阵是奇异，就是结构存在刚体位移(刚体运动 。 必须施加适当的边界条件，消除可能的刚体运动，才能得到唯一解。 必须施加适当的边界条件，消除可能的刚体运动，才能得到唯一解。
第二章 有限元法介绍
2-1 有限单元法的概念 2-2 杆系结构有限单元法 2-3 需要注意的三个问题
2-1 有限单元法的概念
通过材料力学求解和有限元求解进行比较 例：等截面直杆在自重作用下的拉伸 图(a) 单位杆长重量为q 杆长为L 截面面积为A 弹性模量为E 单位杆长重量为q，杆长为L，截面面积为A，弹性模量为E
施加边界条件的另一种方法如下。例如：节点 的位移取指定的值 的位移取指定的值， 施加边界条件的另一种方法如下。例如：节点1的位移取指定的值，即
取一个很大的数á， 子阵上， 取一个很大的数 ，把 á[ I ] 加到总体刚度矩阵的 [K11 ]子阵上， 子阵上 对应的方程右端改写为á{ 对应的方程右端改写为 d1}，则 ，
i+1 图 2-3
图 2-2
2-1 有限单元法的概念
有限单元法求解直杆拉伸： 有限单元法求解直杆拉伸：
3、假设线单元上的位移为线性函数
u x
x i−1
i-1
u i −u i−1 u = u (x) = u i−1 + ( X − X i−1) Li
du ui −ui−1 εx = = dX Li ui −ui−1 σi = Eεi = E( ) Li ui −ui−1 N i = A σi = A E ( ) Li
5.已知外载荷 已知外载荷
求解方程组得到{d3},{d4} 。 求解方程组得到 6. 由求解得到的节点位移，可求解各杆件的内力，以及各固定点的支反力 由求解得到的节点位移，可求解各杆件的内力，
平面梁单元在局部坐标系中的刚度方程： 平面梁单元在局部坐标系中的刚度方程：
2.3需要注意的三个问题 需要注意的三个问题
(2-9)
在坐标系x′y′中单元每个节点只有一个自由度，即沿轴向的位移，所以 中单元每个节点只有一个自由度，即沿轴向的位移， 在坐标系 中单元每个节点只有一个自由度 每个单元共有2个自由度 每个单元共有 个自由度 ui′ , u′j ，每个节点对应的端点力也只有一个 分量， 分量，即单元的刚度方程变得十分简单
2以节点位移表示单元节点力
fi T {δ i } = = [ fi θi ] θi
(2-1)
T
{δ }
e
= fi θi f j θ j
(2-2)
{ p}
qi { pi } = mi
(2-3)
T
e
= qi mi q j m j
4） 施加适当的边界条件。 ） 施加适当的边界条件。 5）求解方程组，得到离散参数的值。 ）求解方程组，得到离散参数的值。 6）如果需要，进行附加计算，如计算元件的变形、内力、约束点的支反力等。 ）如果需要，进行附加计算，如计算元件的变形、内力、约束点的支反力等。
2.4例题 例题1 例题
所示直梁，已知E， ， ， 如图 所示直梁，已知 ，I，Z，M,AB=BC=CD=I, IAC=2I,ICD=I
如果á的数值比所有刚度系数大得多， 如果 的数值比所有刚度系数大得多，上式可近似写成 的数值比所有刚度系数大得多
3. 求解离散问题
一般步骤可归纳如下： 一般步骤可归纳如下： 1) 确定一组离散的参数 ，它们能够描述整个系统及其中各个单元的行为特点， 确定一组离散的参数{a}，它们能够描述整个系统及其中各个单元的行为特点， 叫做系统的参数(变量 变量)。 叫做系统的参数 变量 。 2) 对单元进行分析，由单元的物理性质，连续性和平衡条件，用系统参数 e} 对单元进行分析，由单元的物理性质，连续性和平衡条件，用系统参数{a 表示另一组量{q ， 表示另一组量 e}，它们之间的一般关系为
联立解得
5 qa2 u1 = 2 EA
8 qa2 u2 = 2 EA
9 qa2 u3 = 2 EA
与材料力学的精确解答在结点处完全相同
2-2 杆系结构有限单元法
考虑一个平面桁架结构，如图所示： 考虑一个平面桁架结构，如图所示：
每个节点处只有两个位移分量，整个问题的自由度数是节点数的二倍。 每个节点处只有两个位移分量，整个问题的自由度数是节点数的二倍。
1.划分单元 划分单元
两个节点之间的杆件构成一个单元，杆件结构的节点可以按以下原则选取 两个节点之间的杆件构成一个单元，杆件结构的节点可以按以下原则选取:
• • • • • •
(1)杆件的交点一定要取为节点。 杆件的交点一定要取为节点。 杆件的交点一定要取为节点 (2)阶梯形杆截面变化处一定要取为节点。 阶梯形杆截面变化处一定要取为节点。 阶梯形杆截面变化处一定要取为节点 (3)支撑点和自由端要取为节点。 支撑点和自由端要取为节点。 支撑点和自由端要取为节点 (4)集中载荷作用处要取为节点。 集中载荷作用处要取为节点。 集中载荷作用处要取为节点 (5)欲求位移的点要取为节点。 欲求位移的点要取为节点。 欲求位移的点要取为节点 (6)单元长度不要相差太多。 单元长度不要相差太多。 单元长度不要相差太多
令
N i+1
图 2-5
将位移和内力的关系代入得
(2-1)
q 1 2 −u i-1 + (1+ λi ) u i − λiu i+1 = (1+ ) Li λi 2EA
用结点位移表示的平衡方程，其中i=1， 用结点位移表示的平衡方程，其中i=1，2，… n有n个方程 i=1 未知数也有n 解方程组，得出结点位移， 未知数也有n个，解方程组，得出结点位移，进而计算应力
1.坐标转换 坐标转换 2.边界条件 2.边界条件 3.求解离散问题 求解离散问题
1. 坐标转换
在上面分析桁架结构时，每个单元有 个自由度 单元的刚度矩阵是4× 的 个自由度， 在上面分析桁架结构时，每个单元有4个自由度，单元的刚度矩阵是 ×4的 如果沿单元的轴向建立一个坐标系x′， ， 如果沿单元的轴向建立一个坐标系 ，y′， 如图2.10所示，x′轴与 轴的夹角为 。 所示， 轴与 轴的夹角为θ。 轴与x轴的夹角为 如图 所示
(2-7)
Z1 M 1 Z2 T M 2 {Q} = = [ Z1 M 1 Z 2 M 2 Z3 M 3 Z 4 M 4 ] Z3 M 3 Z4 M 4
(2-8)
qi = a11 fi + a12θi + a13 f j + a14θ j mi = a21 fi + a22θi + a23 f j + a24θ j q j = a31 fi + a32θi + a33 f j + a34θ j m = a f + a θ + a f + a θ 41 i 42 i 43 j 44 j j
∆(dx) =
EA
=
EA
x截面上的位移： 截面上的位移： 截面上的位移
u=∫
x 0 x q(L − x)dx N(x)dx q x2 =∫ = (Lx − ) 0 EA EA EA 2
根据几何方程求应变，物理方程求应力。 根据几何方程求应变，物理方程求应力。这里 du q εx = = (L − X) dX EA q σx = Eεx = (L − X) EA
写成矩阵形式有
每个单元的轴线方向不同，所以各自的轴向位移和端点力方向不同， 每个单元的轴线方向不同，所以各自的轴向位移和端点力方向不同， 各自的轴向位移和轴向端力的方向如下图所示： 各自的轴向位移和轴向端力的方向如下图所示：
局部坐标系和总体坐标系中节点位移之间的关系有
同样有
在两个不同坐标系中单元端点力在相应节点位移上所做的功相等
0
u
x L
L-x
N
N dx N
L 3 L 3 L 3
5 qa2 2 EA
L a= 3
8 qa2 2 EA
9 qa2 2 EA
x
X (a)
(b) 图 2-1
(c)
2-1 有限单元法的概念
材料力学方法求解直杆拉伸： 材料力学方法求解直杆拉伸： 图(b)---位移法 位移法 考虑微段dx， 考虑微段 ，内力 N=q (L-x) dx的伸长为 的伸长为 N(x)dx q(L − x)dx
(2-4)
节点载荷
Zi T {Qi } = = [ Zi M i ] M i
(2-5)
{Q}
e
= Zi M i Z j M j
T
(2-6)
f1 θ 1 f2 T θ 2 {δ } = = [ f1 θ1 f 2 θ 2 f3 θ3 f 4 θ 4 ] f3 θ3 f4 θ 4