本文分析和探讨使用 OpenSees 时可能会遇到的 5 个问题,分别是

  1. 如何在时程分析中改变步长
  2. twoNodeLink zeroLength 单元的正负向
  3. 如何施加预拉力
  4. twoNodeLink 单元的特性
  5. recorder 输出的几个问题

第 1 问 如何在时程分析中改变步长

减小时程分析步长可以改善收敛性

首先,建立一个简单的弹性时程分析模型。模型采用单自由度,两个节点,位置重合,用 twoNodeLink 单元相连。 twoNodeLink 单元中赋予材料属性,用单轴材料模型 Elastic 来实现。具体代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
model BasicBuilder -ndm 1 -ndf 1
node 0 0.
node 1 0.
fix 0 1
mass 1 1.
uniaxialMaterial Elastic 1 1.0e5
element twoNodeLink 1 0 1 -mat 1 -dir 1
timeSeries Path 1 -dt 0.01 -filePath accel.txt
pattern UniformExcitation 1 1 -accel 1 -fact 9800
rayleigh 0. 0. 0. [expr 2*0.02/pow([eigen -fullGenLapack 1],0.5)]
recorder Node -file node.txt -time -node 1 -dof 1 disp

constraints Transformation
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator Newmark 0.5 0.25
analysis Transient
analyze 3000 0.01

由于本例中结构为弹性,所以分析很快结束。但是当结构非线性较强时,可能会发生不能收敛的情况。例如,将上面代码中

1
uniaxialMaterial Elastic 1 1.0e5

改为

1
uniaxialMaterial Steel01 1 200. 1.0e5 0.1

即增加材料非线性,使 twoNodeLink 单元的内力达到 200 后,刚度褪化为初始刚度的 0.1 倍。这时再运行,得到以下错误信息

1
2
3
4
5
WARNING: CTestNormDispIncr::test() - failed to converge
after: 10 iterations current Norm: 0.0072 (max: 1e-08, Norm deltaR: 360)
NewtnRaphson::solveCurrentStep() -the ConvergenceTest object failed in test()
DirectIntegrationAnalysis::analyze() - the Algorithm failed at time 1.82
OpenSees > analyze failed, returned: -3 error flag

即迭代不收敛。这时很容易想到,由于增加了几何非线性,模型收敛更加困难。一个很容易的解决方法是,把原代码中的

1
analyze 3000 0.01

改成

1
analyze 6000 0.005

即缩小每一分析步的步长。这时发现,模型收敛。

计算时间与收敛性的平衡

尽管减小分析步步长可以改善收敛性,但是会增长计算时间,这就形成了一对矛盾。在地震作用较小的时段,结构非线性发展很小,可以用较大的步长快速算过。但是在地震作用较大的时段,结构可能产生较强非线性,所以要用足够小的步长来保证收敛。

发了平衡模型分析所用的总时间和模型的收敛性,就要在运算的过程中动态地改变步长,以达到平衡。

下面介绍一种动态改变分析步长的方法。代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
proc timeHistoryAnalyze {npts dt tol} {
# 输入参数:
# npts: 地震文件中数据点的个数
# dt: 地震文件中数据点之间的时间间隔
# tol: 最小允许的时间步长

# 当前分析时间
set current 0.0
# 是否收敛 收敛时ok=0
set ok 0
# 增加的步数,用于屏幕输出
set incr 0
# 在指定时间范围内循环,如果不收敛则退出
while {$ok == 0 && $current < [expr $npts * $dt] } {
# 设置初始时间间隔与地震记录时间间隔相同
set interval $dt
# 以$interval时间间隔分析一步
set ok [analyze 1 $interval]
# 计数器加1
set incr [expr $incr + 1]
# 如果刚才的分析步不收敛,减小步长,但不小于定义的最小步长
while {$ok != 0 && $interval>1.0e-7} {
# 步长减半
set interval [expr $interval/2]
puts "Fail to convenge at time $current, set a new step length = $interval"
# 重新分析
set ok [analyze 1 $interval]
# 计数器加1
set incr [expr $incr + 1]
}
# 当前时间增加
set current [expr $current + $interval]
# 每10步输入一次信息
if {[expr $incr % 10] == 0} {
puts "analyzing at time $current"
}
}
puts "analyze done with code $ok!"
}

将上述代码放入文件中,就可以在分析中调用了。保存文件名为 opq1-dynamic-interval.tcl ,再将原代码中的

1
analyze 3000 0.01

改为

1
2
source opq1-dynamic-interval.tcl
timeHistoryAnalyze 3000 0.01 1.0e-4

运行,从控制台中可以看出,在某些分析步中第一次使用 0.01 步长没有收敛,程序自动将步长调整为 0.005 ,实现了收敛。

源码下载

本问的源代码可以点击下载

注意:除了减少时间步长之外,还有其它的改善收敛性方法,请参考:如何使 OpenSees 模型收敛?SmartAnalyze 帮助您

正负向问题何时产生

有的时候,一些简单的接触问题中,可能会采用简化的模拟方法,用只压不拉的模型来近似模拟接触。即在接触面所定义的两个点相互挤压时,单元产生一个较大的刚度,而两个点相互脱开时,刚度几乎没有。在 OpenSees 中,有一个简单的单轴材料模型,即 ENTElastic No Tension Material),可以较好地模拟这一特性。

通常,可以把这一材料模型作为力-位移关系模型来使用,赋于 zeroLength 单元或 twoNodeLink 单元中。这两个单元都支持 零长度 ,而通常接触面会被自然地视为是一个零长度的单元。那么此时就会产生问题,对于 零长度 的单元,如果采用 ENT 模型,什么方向是拉,什么方向是压呢?

算例

下面代码代表的模型,是平面模型,用两个 twoNodeLink 将两个位置相同的节点连接起来,其中一个赋予 ENT 材料,另一个赋予 Elastic 材料,以避免结构刚度矩阵奇异。两个节点中一个固定,另一个被施加正向的单位位移。观察结构的力-位移关系。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 0. 0.
fix 0 1 1 1
fix 1 0 1 1
uniaxialMaterial Elastic 1 10
uniaxialMaterial ENT 2 10
element twoNodeLink 1 0 1 -mat 1 -dir 1
element twoNodeLink 2 0 1 -mat 2 -dir 1
pattern Plain 1 Linear {
load 1 1. 0. 0.
}
recorder Node -file node.txt -time -node 1 -dof 1 disp
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator DisplacementControl 1 1 0.2
analysis Static
analyze 5

输出的文件 node.txt 的内容为

1
2
3
4
5
2 0.2
4 0.4
6 0.6
8 0.8
10 1

可以看出,此时结构刚度为 10,与 Elastic 材料的刚度相同。

现在把位移控制语句

1
integrator DisplacementControl 1 1 0.2

改为

1
integrator DisplacementControl 1 1 -0.2

即赋予一个负向的位移,这时的 node.txt 文件的内容为

1
2
3
4
5
-4 -0.2
-8 -0.4
-12 -0.6
-16 -0.8
-20 -1

可见此时结构的刚度发生了变化。为 ENTElastic 材料并联后的刚度。

那么同样是将两个重合的点拉开,为什么正向拉开被 ENT 材料认为是受拉,而负向拉开被认为是受压呢?很容易想到,这与定义单元时,两个节点 iNodejNode 的定义顺序有关。

为了验证,把原代码中

1
element twoNodeLink 2 0 1 -mat 2 -dir 1

改为

1
element twoNodeLink 2 1 0 -mat 2 -dir 1

仍采用负向加载,得到的 node.txt 结果为

1
2
3
4
5
-2 -0.2
-4 -0.4
-6 -0.6
-8 -0.8
-10 -1

可以发现,刚度确实发生了变化。

采用 zeroLength 单元替换 twoNodeLink 单元,得到的结果完全相同。

避免节点顺序的影响

虽然发现了这一特性,可以通过仔细定义两个节点的顺序来解决这一问题,但还是不尽人意。所以在实际使用中,建议将两个节点隔开一点距离,避免完全重合。这样,两个节点距离减小的方向就是压缩的方向,两个节点距离增大的方向就是拉伸的方向,与节点的定义顺序无关。

将上面代码中

1
node 1 0. 0.

改为

1
node 1 0.01 0.

再进行尝试,发现不论顺序如何,拉和压正常了。

源码下载

本问的源代码可以点击下载

第 3 问 如何施加预拉力

用材料模型施加预拉力

在一些有限元软件中,施加预拉力可以通过降温法方便地完成。然而,由于 OpenSees 的热学模块还不完善,不能用这一方法来施加预拉力,但是, OpenSees 中的单轴材料模型 Steel02 可以帮助用户施加预拉力。

Steel02 代表Menegotto-Pinto模型,有很多参数可以选择。用于施加预拉力时,最重要的参数是 $sigInit ,也就是初始应力。定义了初始应力,也就相应地定义了施加的预拉力。通常 Steel02 可以与 trusscorotTrusstwoNodeLink 等单元配合使用。

预拉力定义方法

下面通过一个算例来介绍有预拉力模型的建立方法。在平面两自由度体系中,设置两个节点,用一个 truss 单元来模拟受压结构,再用另一个 truss 单元来模拟一个后张法预应力钢绞线,张拉于简支梁的两端。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
model BasicBuilder -ndm 2 -ndf 2
node 0 0. 0.
node 1 2000. 0.
fix 0 1 1
fix 1 0 1
uniaxialMaterial Elastic 1 2.0e5
# Steel02 Arguments #tag #Fy #E #b #R0 #cR1 #cR2 #a1 #a2 #a3 #a4 #sigInit
uniaxialMaterial Steel02 2 2.0e4 2.0e6 0.1 15 0.925 0.15 1.0 1.0 1.0 1.0 2.0e3
# truss #tag iNode #jNode #Area #Mat
element truss 1 0 1 1.0e5 1
element truss 2 0 1 500. 2
recorder Node -file node.txt -time -node 1 -dof 1 disp
recorder Element -file ele1.txt -time -ele 1 force
recorder Element -file ele2.txt -time -ele 2 force
constraints Plain
numberer Plain
system BandGeneral
test EnergyIncr 1.0e-6 200
algorithm Newton
integrator LoadControl 1
analysis Static
analyze 1

对于 Steel02 的参数,重要的几个是 FyEsigInit 。由于不希望预应力钢绞线出现屈服,所以一些与材料非线性有关的参数在这里暂不细致考虑。用 truss 单元来模拟钢绞线,通过分析 1 步,可以实现预拉力的施加。

在本模型中,两个 truss 的刚度分别为

\[
K_1 = \frac{E_1A_1}{L} = \frac{2.0\times 10^5\times 1\times 10^5}{2000} =1\times 10^7
\]
\[
K_2 = \frac{E_2A_2}{L} = \frac{2\times 10^6\times 500}{2000}=5\times 10^5
\]

定义了钢绞线的初应力为

\[
\sigma_p = 2\times 10^3
\]

在施加预拉力时,钢铰线产生了初始变形

\[
\Delta_p = \frac{F_p}{K_2}=\frac{\sigma_p A_2}{K_2}
\]

张拉后,预拉力由 1 单元承担。根据力的平衡,设\(F_1\)为 1 单元中的压力,\(F_2\)为 2 单元中的拉力,则有

\[
F_1 = F_2
\]

在此过程中,1 单元和 2 单元同时产生了压缩变形\(\Delta\),则此时钢铰线的实际变形为

\[
\Delta_2 = \Delta_p -\Delta
\]

根据物理关系,有

\[
F_1 = K_1 \Delta
\]
\[
F_2 = K_2 \left(\Delta_p-\Delta\right)
\]

代入上式,得

\[
K_1 \Delta = K_2 \left(\frac{\sigma_p A_2}{K_2}-\Delta\right)
\]

解得

\[
\Delta = \frac{\sigma_p A_2}{K_1+K_2}=0.0952
\]
\[
F_1 = F_2 = K_1 \Delta=952000
\]

查看输出的 node.txt 文件

1
1 -0.0952381

ele1.txt 文件

1
1 952381 0 -952381 -0

与计算的结果吻合。

施加预拉力后的加载

在实际结构分析中,在施加预拉力后,通常要进行进一步加载。在本问中以轴向加载为例,探讨进一步加载的方法。

在上一小节模型的基础上,添加如下代码

1
2
3
4
5
pattern Plain 1 Linear {
load 1 1. 0.
}
integrator DisplacementControl 1 1 0.01
analyze 5

即将结构的 1 节点向正向按位移加载,此时得到 node.txt 文件为

1
2
3
4
5
6
1 -0.0952381
78750.5 -0.0852381
170625 -0.0752381
269063 -0.0652381
370781 -0.0552381
474141 -0.0452381

我们知道,两个单元并联,则并联后的刚度为

\[
K=K_1+K_2=1.05\times 10^7
\]

在第一次移动 0.01 位移的时候,施加的外力为

\[
F=Kd_1=1.05\times 10^7\times 0.01 = 1.05\times 10^5
\]

但是从 node.txt 中得到的外力为 78750.5 ,与计算的结果不一致。其中的错误可能是因为,在加载过程中,在每一分析步运行时,预应力都被通过某种机制重复地施加了一次。详细原因这里不作细究。解决这一问题的方法很简单,只要在分析之间加上一句话

1
loadConst -time 0.0

这时再运行,得到 node.txt

1
2
3
4
5
6
1 -0.0952381
105000 -0.0852381
210000 -0.0752381
315000 -0.0652381
420000 -0.0552381
525000 -0.0452381

与计算的结果一致。

以上分析表明,含预拉力模型实际加载前要将 time 归零

源码下载

本问的源代码可以点击下载

OpenSees 的单元库中,目前的 连接单元 目录下只有一个单元 twoNodeLink 。这个单元可以把单轴应力-应变材料关系模型直接作为结构的力-位移关系模型,有广泛的应用情景。但是这一单元有很多难以把握的特性,如果不能充分理解,很容易出错。因此对 twoNodeLink 单元的特性要正确地理解。

低维度空间中指定局部坐标轴方向

在有限单元法中,每一个单元都有自己的局部坐标系。如果定义为零长单元,则单元的局部坐标方向与整体坐标方向相同。如果为非零长单元,则单元的局部 1 轴为 iNode 指向 jNode 方向。但是局部 2 轴的方向没有指定。这就需要用户来指定。语句为

1
-orient <x1 x2 x3> y1 y2 y3

如果用户不指定 y 方向,程序会输出警告信息

1
2
WARNING TwoNodeLink::setUp() - element: 1
no local y vector specified

在指定方向的语句中, y1y2y3 分别代表 2 轴方向向量的坐标(在全局坐标系下)。用户也可以用同样的方式指定局部 1 轴的坐标,即赋予语句中 x1x2x3 的值。

那么这时产生了一个问题。对于三维空间中的模型,定义一个方向向量需要三个坐标,但是如果在二维空间中,如何自定义 twoNodeLink 单元的方向呢?

下面我们建立一个算例,来探讨这个问题。

在二维三自由度空间中,建立一个两节点模型,用 twoNodeLink 单元相连。代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 1. 0.
fix 0 1 1 1
fix 1 0 1 1
uniaxialMaterial Elastic 1 1.
element twoNodeLink 1 0 1 -mat 1 -dir 1
recorder Node -file node.txt -time -node 1 -dof 1 2 3 disp
pattern Plain 1 Linear {
load 1 1. 0. 0.
}
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator LoadControl 1
analysis Static
analyze 5

查看 node.txt 文件,得到

1
2
3
4
5
1 1 0 0
2 2 0 0
3 3 0 0
4 4 0 0
5 5 0 0

这个模型没有太多实际意义,只是告诉我们,当前没有自定方方向时,局部 1 轴是两个节点的连线方向。那么如果想指定局部 1 轴的方向怎样处理?我们建立一个新的模型。

在这个模型中,希望把 twoNodeLink 单元的局部 1 轴坐标赋于全局坐标的 2 方向上。那就要用两个方向向量来分别定义局部 1 轴和局部 2 轴。在二维空间中,每个方向向量有两个坐标。把代码改为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 1. 0.
fix 0 1 1 1
fix 1 1 0 1
uniaxialMaterial Elastic 1 1.
element twoNodeLink 1 0 1 -mat 1 -dir 1 -orient 0. 1. 1. 0.
recorder Node -file node.txt -time -node 1 -dof 1 disp
pattern Plain 1 Linear {
load 1 0. 1. 0.
}
constraints Plain
numberer Plain
system BandGeneral
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator LoadControl 1
analysis Static
analyze 5

这时程序会报错

1
2
WARNING insufficient arguments after -orient flag
twoNodeLink element: 1

既然提示输入的变量不够,在二维空间中,默认所有 z 坐标都为 0,那就把 z 坐标为 0 加入进去,把单元定义语句改为

1
element twoNodeLink 1 0 1 -mat 1 -dir 1 -orient 0. 1. 0. 1. 0. 0.

此时再运行,控制台输出的信息表明,已经指定了方向。

1
2
WARNING TwoNodeLink::setUp() - element: 1
ignoring nodes and using specified local x vector to determine orientation

查看 node.txt 文件,发现结果正确

1
2
3
4
5
1 0 1 0
2 0 2 0
3 0 3 0
4 0 4 0
5 0 5 0

进一步探究发现,如果 z 坐标不为 0,得到的结果会出现错误。因此,在对 twoNodeLink 单元赋方向时,如果坐标数小于 3,应对不存在的坐标值赋 0

用户定义方向坐标不垂直

如果用户定义的两个方向坐标不垂直,程序如何处理呢?

在上述二维模型的基础上做一些小的修改,输出结构的初始刚度矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 1. 0.
fix 0 1 1 1
uniaxialMaterial Elastic 1 1.
uniaxialMaterial Elastic 2 2.
uniaxialMaterial Elastic 3 3.
element twoNodeLink 1 0 1 -mat 1 2 3 -dir 1 2 3 -orient 1. 0. 0. 0. 1. 0.
constraints Plain
numberer Plain
system FullGeneral
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator LoadControl 1
analysis Static
analyze 1
printA

得到

1
2
3
1 0 0
0 2 -1
0 -1 3.5

调整 y1y2 两个坐标的值,使局部两轴不正交,发现得到的刚度矩阵不变。调整 x1x2 两个坐标的值,发现刚度矩阵改变。

上面的小例子说明,如果定义的两个方向不正交,程序是以局部 1 轴的定义为准的。实际上,程序是采用定义的两个向量叉积作为局部 3 轴,再用局部 3 轴与局部 1 轴的叉积作为局部 2 轴的。与 geomTransf 定义坐标的方法类似。

观察上一小节中得到的刚度矩阵,可以看出, twoNodeLink 单元的局部 2、3 自由度是耦合的,这在使用的过程中要特别注意。如果不希望发生自由度耦合的话,应该定义多个在 1 方向上的 twoNodeLink 单元,而不是使用一个单元的多个方向。

对于程序给出的刚度矩阵

\[
A_{ij}=\left[\begin{matrix}1&0&0\\0&2&-1\\0&-1&3.5\end{matrix}\right]
\]

我们对矩阵中的元素一一进行考查。
对于\(A{11}\)和\(A{22}\),很明显是我们定义的两个刚度。但是对于\(A{33}\),这个刚度与我们定义的刚度是不同的,因为发生了自由度耦合。这一点要特别注意。\(A{23}\)和\(A_{32}\)元素不为 0,也是由于自由度耦合造成的。那么这几个刚度值是如何计算出来的呢?

在这里,有一个 shear distance 的概念,直接影响了结构的刚度矩阵。 twoNodeLink 的原理如下图所示, shear distance 主要影响的是代表局部 2 轴刚度,也就是图中\(u_{b,2}\)弹簧的位置。在默认状态下, shear distance 定义为单元的中部,即如图所示位置。

twoNodeLink

此时,如果 iNode 固定,而 jNode 产生单位转角时,可以看到,对于\(u{b,2}\)弹簧所连接的两个刚臂,与 iNode 一端相连的刚臂是固定的,而与 jNode 一端相连的刚臂产生了一个转角。由于 jNode 的平动自由度被约束,在\(u{b,2}\)弹簧位置就会产生一个相对位移,大小为转角与转动刚臂长度的乘积(小变形情况下)。在本例中,刚臂长度为 0.5,\(u_{b,2}\)弹簧的刚度为 2,因此刚度矩阵中

\[
A_{23}=A_{32}=-0.5\times2=-1
\]

同时,刚臂转动后\(u_{b,2}\)弹簧产生的力会产生一个附加弯矩,作用在 jNode 上,大小为

\[
M_s = 1\times 0.5=0.5
\]

与定义的转角弹簧刚度叠加,得到刚度矩阵中

\[
A_{33}=3+0.5=3.5
\]

剪切距离的定义

剪切距离 shear distance 还可以由用户自定义。用剪切中心点与 iNode 之间的距离与单元长度之比来确定。在二维空间中定义一个区间\([0, 1]\)中的数,在三维空间中定义两个区间\([0, 1]\)中的数。默认状态下这个数取值为 0.5。

将上面代码中单元定义的部分改成

1
element twoNodeLink 1 0 1 -mat 1 2 3 -dir 1 2 3 -orient 1. 0. 0. 0. 1. 0. -shearDist 0.2

得到输出的刚度矩阵为

1
2
3
1 0 0
0 2 -1.6
0 -1.6 4.28

发现与上一小节中得出的结论一致。

源码下载

本问的源代码可以点击下载

第 5 问 recorder 输出的几个问题

特征值与特征向量输出

想要获得结构的特征向量,需要使用 recorder 输出。以一个三维空间中的两自由度体系为例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 0. 2000.
node 2 0. 4000.
fix 0 1 1 1
mass 1 200. 200. 0.
mass 2 100. 100. 0.
geomTransf Linear 1
element elasticBeamColumn 1 0 1 1.0e4 2.0e6 1.0e8 1
element elasticBeamColumn 2 1 2 1.0e4 2.0e6 1.0e8 1
puts [eigen 2]
recorder Node -file eigen1.txt -node 1 2 -dof 1 2 3 eigen1
recorder Node -file eigen2.txt -node 1 2 -dof 1 2 3 eigen2
record

运行,控制台中输出了两个数字

1
7.782659e+01  2.065031e+03

这就是结构的两个特征值。在 OpenSees 中,特征值\(\lambda\)的定义为

\[
\left(K-\lambda M\right)\phi=0
\]

即我们通常所说的\(\omega^2\)。它与结构周期的关系为

\[
T=\frac{2\pi}{\omega}=\frac{2\pi}{\sqrt{\lambda}}
\]

因此,算例结构的两阶周期为

\[
T_1 = \frac{2\pi}{\sqrt{77.83}}=0.71s
\]

\[
T_2 = \frac{2\pi}{\sqrt{2065}}=0.14s
\]

对应两个周期的振型,可以从 recorder 生成的两个文件中找到。打开 eigen1.txt ,内容为

1
0.029707 0 -2.58115e-05 0.0907468 -1.35525e-20 -3.28741e-05

由于 recorder 中输出了两个节点,每个节点输出了三个自由度,所以共输出了 6 个数字,前三个为节点 1 在三个自由度上对应 1 阶周期的振型,后三个为节点 2 的振型。可以看出,在 1 阶振型中,2 方向和 3 方向的位移都远小于 1 方向,可以忽略。只关注在 1 方向上的位移即可。这个位移值与振型的值是一致的。即 1 阶振型为

\[
\phi_1 = \left\{ \begin{matrix}0.029707&0.0907468\end{matrix}\right\}^T
\]

可以根据用户的需要进行归一化。从 eigen2.txt 文件中可以找到 2 阶振型

\[
\phi_2 = \left\{ \begin{matrix}-0.0641677&0.0420121\end{matrix}\right\}^T
\]

值得一提的是,如果所求的特征值数与结构的振动自由度数相同,要将特征值输入语句改为

1
eigen -fullGenLapack 2

否则程序会报错。这在小模型算例中可能会用到。

节点加速度输出

在地震响应分析中,我们关心结构的位移和速度通常是相对于地面的相对位移相对速度,而关心结构的加速度通常是绝对加速度,因为只有绝对加速度才能与牛顿第二定律建立关联。

在用 recorder 输出时程分析结果时,应该特别注意输出的是相对值还是绝对值。对于位移和速度,输出的都是相对值,和我们关注的一样。而对于加速度,要特别注意。

如果采用语句

1
recorder Node -file node.txt -node 1 -dof 1 accel

则输出的是相对加速度。如果希望输出绝对加速度,则必须将语句改为

1
recorder Node -file node.txt -timeSeries 1 -node 1 -dof 1 accel

这时输出的是绝对加速度。

拟时间输出

将第一节中的模型稍作修改,得到一个悬臂柱水平位移加载的模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 0. 2000.
node 2 0. 4000.
fix 0 1 1 1
geomTransf Linear 1
element elasticBeamColumn 1 0 1 1.0e4 2.0e6 1.0e8 1
element elasticBeamColumn 2 1 2 1.0e4 2.0e6 1.0e8 1
pattern Plain 1 Linear {
load 2 1000. 0. 0.
}
recorder Node -file node.txt -time -node 2 -dof 1 disp
recorder Element -file ele.txt -time -ele 1 force
constraints Plain
numberer Plain
system BandGeneral
test EnergyIncr 1.0e-6 200
algorithm Newton
integrator DisplacementControl 2 1 0.1
analysis Static
analyze 5

注意模型中对 pattern 的定义是必不可少的。在 OpenSees 中, DisplacementControl 的实现采用了一种很巧妙的机制。首先,要用 pattern 命令定义结构中的荷载,这些定义为参照荷载。用\(P{ref, i}\)来表示。在本例中,仅定义了一个水平向,大小为 1000 的参照荷载\(P{ref, 1}\)。每个参照荷载可以定义一个对应的荷载放大系数\(\lambda_i\)。在本例中没有指定放大系数,则\(\lambda_1=1\)。

对于结构所受的实际外荷载,可以用下式描述。

\[
P_f = P_c + \lambda\sum_{i=1}^{np}{\lambda_iP_{ref, i}}
\]

式中的\(\lambda\)实际上就是所谓的拟时间(pseudo-time),也就是 recorder-time 输出的对应值。对于\(P_c\),实际上是指不参与被拟时间放大的荷载,比如第三问中的钢铰线预拉力,再比如低周往复加载时施加的重力荷载。这可以通过以下语句实现。

1
loadConst -time 0.0

在该语句之前的荷载,不再被拟时间放大。同时拟时间也被归零。其它荷载都到乘上拟时间之后施加在结构上。

打开 node.txt ,可以看到输出结果

1
2
3
4
5
0.9375 0.1
1.875 0.2
2.8125 0.3
3.75 0.4
4.6875 0.5

第一列表示的是拟时间。也就是所有荷载的乘数。第二列是节点 2 的 1 方向位移。第一行数据的意思是,当节点 2 在 1 方向产生 0.1 的位移是,拟时间为 0.9375,也就是结构所受的水平力为 0.9375×1000=937.5。依次类推。

通过 node.txt 文件可以输出结构的荷载-位移曲线。

单元内力输出

打开上一小节中生成的 ele.txt 文件,内容如下。

1
2
3
4
5
0.9375 -937.5 0 3.75e+06 937.5 0 -1.875e+06
1.875 -1875 0 7.5e+06 1875 0 -3.75e+06
2.8125 -2812.5 0 1.125e+07 2812.5 0 -5.625e+06
3.75 -3750 0 1.5e+07 3750 0 -7.5e+06
4.6875 -4687.5 0 1.875e+07 4687.5 0 -9.375e+06

其中第 1 列为拟时间,第 2 到 7 列为单元受力。前三列为 iNode 节点受的等效外力,后三列为 jNode 节点受的等效外力。分别按 1、2、3 三个方向排序。这些力是定义在全局坐标系下的。单元力输出时,输了上节中使用的 force 一般还有另外两种选择: localForcebasicForce

下面我们把上一小节中 elasticBeamColumn 单元替换成 dispBeamColumn ,并赋予一个 250×250×14×14 的 H 形纤维截面。一种快速生成纤维截面的方法是使用本站[结构工具箱][9]中 OpenSees 前处理模块的[纤维截面生成工具箱][11]。

生成的纤维截面的代码为

1
2
3
4
5
6
7
8
9
10
11
12
13
section Fiber 1001 {
#200x200x14x14
fiber 97.6667 -95.0000 46.6667 1
fiber 93.0000 -95.0000 46.6667 1
fiber 88.3333 -95.0000 46.6667 1
...
fiber -93.0000 95.0000 46.6667 1
fiber -88.3333 95.0000 46.6667 1
}
uniaxialMaterial Elastic 2001 404391249000
uniaxialMaterial Elastic 4001 2153846153.8462
uniaxialMaterial Elastic 3001 3589743589.7436
section Aggregator 1 3001 Vy 4001 Vz 2001 T -section 1001

该截面中聚合了抗剪和抗扭刚度。用这一截面替换上一小节的弹性截面,得到模型代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
model BasicBuilder -ndm 2 -ndf 3
node 0 0. 0.
node 1 0. 2000.
node 2 0. 4000.
fix 0 1 1 1
geomTransf Linear 1
uniaxialMaterial Elastic 1 2.0e5
source opq05-section-200x200x14x14.tcl
element dispBeamColumn 1 0 1 5 1 1
element dispBeamColumn 2 1 2 5 1 1
pattern Plain 1 Linear {
load 2 1000. 0. 0.
}
recorder Node -file node.txt -time -node 2 -dof 1 disp
recorder Element -file ele.txt -ele 1 force
recorder Element -file ele-local.txt -ele 1 localForce
recorder Element -file ele-basic.txt -ele 1 basicForce
constraints Plain
numberer Plain
system BandGeneral
test EnergyIncr 1.0e-6 200
algorithm Newton
integrator DisplacementControl 2 1 0.1
analysis Static
analyze 5

模型中使用了三种单元 recorder ,分别输出 forcelocalForcebasicForce

ele.txt 文件中内容为

1
2
3
4
5
-45.4835 2.74536e-14 181934 45.4835 -2.74536e-14 -90967
-90.967 5.95345e-14 363868 90.967 -5.95345e-14 -181934
-136.451 1.29617e-13 545802 136.451 -1.29617e-13 -272901
-181.934 9.29857e-14 727736 181.934 -9.29857e-14 -363868
-227.418 -8.92754e-14 909670 227.418 8.92754e-14 -454835

ele-local.txt 文件中内容为

1
2
3
4
5
2.74536e-14 45.4835 181934 -2.74536e-14 -45.4835 -90967
5.95345e-14 90.967 363868 -5.95345e-14 -90.967 -181934
1.29617e-13 136.451 545802 -1.29617e-13 -136.451 -272901
9.29857e-14 181.934 727736 -9.29857e-14 -181.934 -363868
-8.92754e-14 227.418 909670 8.92754e-14 -227.418 -454835

ele-basic.txt 文件中内容为

1
2
3
4
5
-2.74536e-14 181934 -90967
-5.95345e-14 363868 -181934
-1.29617e-13 545802 -272901
-9.29857e-14 727736 -363868
8.92754e-14 909670 -454835

不难看出, force 语句得到的等效力是建立在全局坐标系下的, localForce 语句得到的等效力是建立在单元局部坐标系下的。而 basicForce 中输出的力为轴力和两个端弯矩,同样也是在局部坐标系下建立。

对于 forcelocalForce 的输出都很好理解,那么为什么 basicForce 只输出了两个三个量呢?为什么剪力没有了呢?这与 OpenSees 中采用的有限单元法有关。这一部分比较复杂,将在其它文章中进行介绍。

纤维截面应力应变输出

如果需要知道纤维截面的应力,可以通过 recorder 进行输出。在上例的代码中增加一句

1
recorder Element -file stressStrain.txt -ele 1 section 1 fiber 93.00 -95.00 stressStrain

即可实现纤维截面应力和应变的输出。在增加的语句中, section 1 表示单元的积分截面。由于在单元 1 的定义中取了 5 个积分截面,所以 section 后面的数字可以为 1~5,表示 5 个积分截面。 fiber 后面的两个数字表示所考察纤维的局部坐标。

stressStrain.txt 文件如下

1
2
3
4
5
0.34057 1.70285e-06
0.68114 3.4057e-06
1.02171 5.10855e-06
1.36228 6.8114e-06
1.70285 8.51425e-06

第一列表示应力,第二列表示应变。

源码下载

本问的源代码可以点击下载