藏式古建筑木结构损伤识别的数值模拟

木结构

发布日期:

引 言

藏式古建筑木结构以其浓厚的地域特点和宗教色彩闻名于世。由于建造年代久远,在环境荷载、疲劳效应、腐蚀效应和材料老化等因素的影响下,大多数古建筑结构都存在不同程度的损伤。此外,很多藏式古建筑都作为旅游资源向公众开放,人群荷载对结构的影响也威胁着结构的安全。从遗产保护的角度开展针对藏式古建筑木结构的损伤识别和状态评估理论的研究是十分必要且急迫的。

 

近年来,已有学者对古建筑结构的健康监测和损伤识别进行了实践研究。例如,在上海市“沉香阁”部分梁、柱的临时支撑上布设了光纤光栅传感器,实时监测木结构节点的变形,利用监测结果检验临时支撑的有效性[1]。上海大世界的结构监测采用光纤光栅应变传感器、无线加速度传感器和无线应变传感器分别监测结构的局部裂缝、六角形塔楼的振型变化和部分主要构件的应变等[2]。宁波保国寺大殿也建立了一套长期监测系统,对大殿内各角落的温度、湿度、结构的变形与倾斜、沉降等进行实时监测[3]。上述研究是根据监测数据是否发生突变来预测结构是否发生损伤,具有信息滞后性且受传感器布设位置的影响,仅能对结构局部损伤进行预测与目前广泛应用于大型复杂土木工程结构领域的基于振动测试的结构损伤识别方法[4-9]相比,现有的应用于古建筑的损伤识别方法较为落后。文献[10]提出了基于振动响应灵敏度的损伤识别方法,利用模型修正识别结构损伤的位置和程度。该方法可充分利用动力试验测得的大量响应信息,对测点的位置没有特殊要求,具有一定的工程应用价值。徐伟华等[11-12]利用振动响应灵敏度方法实现了对拉索结构和杆结构的损伤识别。吕中荣等[13]利用响应灵敏度方法实现了对Gascogine天桥的有限元模型修正。笔者在对大量藏式古建筑木结构进行实地勘查的基础上,对梁、柱等构件的残损类型和残损机理进行了汇总和分析。以某典型藏式古建筑多层梁柱排架结构为损伤识别的研究对象,以振动响应灵敏度损伤识别方法为基本理论,采用修正的正则化方法求解识别方程,对该结构梁、柱构件的3种损伤状态进行了数值模拟识别。模拟分析的结果为此类结构实际的损伤识别提供了理论依据。

 

1 藏式古建筑木结构的结构特点

藏式古建筑木结构通常用梁柱组成数列纵向排架,梁上铺密椽,椽上敷设楼板,如图1所示。其中,梁、柱和椽等构件均由木材加工而来,墙体和楼板的主要制作材料为土、石。藏式古建筑木结构的柱式由柱础、柱身、柱头栌斗及上面的托木、弓木和梁等组成,如图2所示。柱脚直接搁置在柱础上,柱础提供竖向支撑力和水平摩擦力。柱础用石材制作而成。早期的建筑柱础是露出地面的,大约在元代以后,柱础大都埋入地面。柱、栌斗、垫木、弓木和梁等构件之间均采用暗销连接;梁与梁之间采用榫卯连接。栌斗、垫木和弓木的存在增大了梁柱之间的接触面,在不增加梁、柱高度的情况下增加了殿堂净空高度,同时弓木和垫木也提高了梁架的稳定性。

>

藏式古建筑木结构建造年代久远,内部构件存在一定的残损现象。通过对大量藏式古建筑木结构进行现场调研,对梁、柱构件的残损类型和残损机理进行了汇总和分析[14-16]。

 

1.1 柱的残损

柱是藏式建筑中最主要的承重构件之一。藏式古建筑中柱子通常采用“收分”截面,即柱截面的尺寸由柱底到柱顶逐渐减小。这样的构造增强了柱子及木构架的稳定性。藏式古建筑木结构柱的残损类型可以分为柱脚腐朽、柱身裂缝和柱身扭转等。

柱脚腐朽如图3(a)所示。由于藏式古建筑的柱础一般都埋在楼地面内,而柱脚则直接浮搁在础石上,因此柱脚位置处很容易受潮腐朽。

柱身裂缝如图3(b)所示。很多柱子都存在柱身径向裂缝,其宽度一般较小,但有些可能会形成径向贯通裂缝。对于宽度较大的裂缝,为了防止柱身被进一步侵蚀,可采用塞木条灌胶或设置铁箍的方法来限制裂缝的开展。柱身裂缝产生的主要原因可归纳为以下几个方面:a.外界作用超过了柱子的承载能力;b.相邻构件的损坏和变形造成柱子偏心受压,当偏心距加大时,柱子的承载能力降低,从而导致柱身开裂;c.环境因素等的影响使木材的材料性能降低,导致构件在长期荷载作用下产生裂缝。柱的扭转如图3(c)所示。由于藏式建筑的结构特点,一些上部荷载分布不均匀的角柱往往会出现一定程度的扭转。由于柱础对柱脚的磨擦作用限制了这种扭转变形,从而导致柱身出现扭转斜裂缝,一般从下至上呈螺旋分布。

>

 

1.2 梁的残损

在藏式古建筑木结构中,梁与柱共同组成纵向排架。梁截面较大,一般情况下都有很高的承载冗余度。梁的残损类型主要有梁体裂缝和梁体过大挠曲。

梁体裂缝如图4(a)所示。梁身的裂缝一般为水平方向,发生于梁的侧面或底面,裂缝宽度一般不大。梁体产生裂缝的主要原因为:a.由于藏式古建筑的构造特点,梁与梁之间采用榫卯连接,因此梁头截面一般要被削弱,容易导致局部应力集中现象;b.木材表面原有的干缩裂缝在受弯状态下会进一步扩展。

梁体过大挠曲:弓木和垫木的存在使梁的净跨减小,因此一般情况下梁的挠曲较小。导致梁产生过大挠度的原因主要有:a.上部荷载过大,如放置经书架、佛像等的位置以及回廊等游人频繁活动的位置处;b.在梁构件中,因原有木材有疵病(如木节、斜纹和扭纹等),随天长日久疵病影响逐渐加大,削弱了截面的抵抗能力,造成挠曲过大。对挠曲过大的梁构件一般采取“支顶”的方法来减小其挠度,如图4(b)所示。

>

 

2 基于响应灵敏度的损伤识别方法

基于振动响应灵敏度的损伤识别方法是利用模型修正来识别结构损伤的位置和程度。假定系统质量参数可以获得,且在模型修正过程中保持不变,把刚度参数作为待修正的物理参数。将结构系统的振动响应展开成刚度参数的泰勒级数,并对振动响应灵敏度进行分析,利用振动响应的试验值与理论值的误差进行迭代求解结构的刚度参数,不断更新结构有限元模型,实现对结构损伤位置和程度的识别。对所识别的结构进行动力实验,可以测得结构的加速度响应¨xm、速度响应xm和位移响应xm。根据所需识别的结构的预估参考模型和实验荷载,采用Newmark-β等直接积分法可以计算出结构的振动响应¨x,x 和x。

M¨x(t)+Cx(t)+Kx(t)=LF(t) (1)

其中:M,C,K 分别为参考模型的质量、阻尼和刚度矩阵;L 为映射向量,将实验荷载F 作用的自由度映射到对应的结构自由度上。

采用瑞利阻尼模型定义阻尼矩阵,即C=

a1M+a2K,其中:a1,a2

为可由模态阻尼比和模态频率求得的常数。

由于参考模型中的物理参数与实际结构相比会存在差别,所以实测的振动响应值与计算的响应值也会不一致。因此,需对结构参考模型进行修正,即不断地更新参考模型的物理参数,使响应实测值与计算值的残差最小,从而迭代求解出损伤结构的物理参数。

假定结构质量可以获得,待识别的结构物理参数以结构刚度为表征。设真实的结构刚度K′与预估参考模型的结构刚度K 的差值为ΔK,其表达式为

K′-K =ΔK =Σ

ne

i=1

αiKi (2)

其中:ne表示单元个数;Ki

表示第i 个单元的刚度

矩阵;αi

为单元刚度差异参数,αi

的存在使得结构响应的计算值与实测值之间存在差异。

由于在结构动力实验中加速度响应易于得到,因此多采用加速度响应来识别结构刚度差异参数

αi。将Δ¨x 定义为加速度响应实测值¨xm

与计算值¨x

之间的差值,即

Δ¨x =¨xm -¨x (3)

Δ¨x 可以用结构刚度差异参数α 近似表示为

Δ¨x =Sα (4)

该式为加速度响应对结构刚度差异参数的一阶

灵敏度方程。由式(4)可知,若已知Δ¨x 和S,则可求解出α,下面介绍确定加速度灵敏度矩阵S的方法。

 

2.1 振动响应对物理参数的灵敏度

为了得到加速度响应对结构刚度差异参数的灵敏度,将式(1)对结构的刚度差异参数αi求导,有

M¨x(t)

αi

+Cx(t)

αi

+Kx(t)

αi

-K

αi

x(t)-a2K

αi

x(t) (5)

其中:,¨x(t)/αi,x(t)/αi

和x(t)/αi

分别为加速度、速度、位移响应对结构刚度差异参数αi的灵敏度。

由于x 和x 可根据式(1)获得,K

αi

Σ n

iαiKi

αi

Σ n

i=1Ki,对比式(1),式(5)的右边项可以看作输入的荷

载向量,因此等式左边的量¨x(t)/αi,x(t)/αi和

x(t)/αi

也可以由直接积分法求得。式(4)中只需要测试自由度的加速度响应对结构刚度差异参数的灵敏度S。由于灵敏度矩阵S 里面的元素是时变的,所以加速度响应对刚度差异参数的灵敏度矩阵也是时变的。对应某一时刻tl,灵敏度矩阵S(l)可以表示为

S(l)=

¨x1(tl)

α1

¨x1(tl)

α2

… ¨x1(tl)

αi

… ¨x1(tl)

αne

¨x2(tl)

α1

¨x2(tl)

α2

… ¨x2(tl)

αi

… ¨x2__________(tl)

αne

    … 

¨xj(tl)

α1

¨xj(tl)

α2

… ¨xj(tl)

αi

… ¨xj(tl)

αne

    … 

¨xNsensor(tl)

α1

¨xNsensor(tl)

α2

…¨xNsensor(tl)

αi

…¨xNsensor(tl)

αne其中:i为第i个单元;j为第j 个传感器;Nsensor

为传感器的个数。

 

2.2 物理参数识别的迭代算法

获得灵敏度矩阵后,根据式(4)得到刚度差异参数α,将其代入式(2),得到参考模型与实际结构间的刚度差异,进而可以修正结构参考模型的刚度。由于式(4)采用的是一阶Taylor近似,因此得到的刚度差异参数α是不准确的,需要利用修正后的结构参考模型重新计算Δ¨x和S,再次利用式(4)求解α,此迭代过程一直进行到满足一定的收敛条件为止。以前、后两个迭代步(第k步和第k+1步)计算出来的αk+1和αk 的相对误差作为迭代过程的收敛条件,即

αk+1-αk

αk+1 ≤Ttolerance (7)

其中:Ttolerance

为接近0的极小值。

最终识别的第i个单元刚度差异参数为各迭代步识别结果的总和,即

αi =α0i

+α1i

+α2i

+…+αn

i (8)

其中:α0i

=0;上标n为迭代步数。

2.3 求解灵敏度方程的修正正则化方法式(4)中的未知数小于方程数,即式(4)属于超定方程,仅具有唯一的最小二乘解。结构参数识别等动力学反问题一般都具有不适定性,即矩阵Sk 的条件数很大,Δ¨xk 的微小扰动将会引起方程解的极大跳动,且条件数对于超定方程解稳定性的影响更为明显。在实际计算中由于测量噪声的存在,这种扰动是不可避免的,因此使用常规最小二乘方法将得不到稳定解。工程中常采用Tikhonov正则化方法[17]求解此类不适定的动力学反问题。其基本思想就是利用一簇与原不适定问题相近的适定问题的解去逼近原有不适定问题的解,即在方程最小二乘解的基础上增加一个解的约束条件,使方程的病态状况得到改善。笔者采用修正的Tikhnov正则化方法来求解式(4)。修正正则化方法对离散性较大的物理参数的识别在收敛性和准确性上都优于传统Tikhonov方法。两种方法的不同点在于修正方法考虑了前k+1步解的范围,把对方程当前步解的约束转变成了对方程前k+1步解的约束,保证了解的物理意义。利用正则化参数λ将约束问题转化为无约束的优化问题,即J(Δαk+1,λ)= Skαk+1-Δ¨xk 22

+λ2 Σ

k+1

j=1

αj-αk,* 2

(9)

当k=0时,α0,*=0

当k≥1时

αk,* =

lup Σ

j=1

( α ) j

( i >lup)

Σ k

j=1

( α ) j

i    llow ≤ Σ

j=1

( α ) j

( i ≤lup)

llow Σ

j=1

( α ) i

( i <l )

low

(10)

其中:Σ

k+1

j=1

αj 为前k+1步刚度差异参数之和;lup

llow

分别为单元刚度差异参数的上、下限值,通常情

况下满足-1≤llow≤0,0≤lup≤1。

对Sk 进行奇异值分解

Sk =UΣVT (11)

其中:U 为nt×nt阶酉矩阵,满足UTU=Int;nt为采

样点总数;V 为ne×ne阶酉矩阵,满足VTV=Ine;Σ

为nt×ne阶矩阵,其对角线元素为降序排列的奇异值σi(i=1,2,...,ne),即σ1≥σ2≥...≥σne≥0,其他元素均为0。

将式(11)代入式(9),得到方程的正则解

Δαk+1 = (ST

kSk+λ2I)-1[ST

kΔ¨xk-λ2(αk-αk,*)]=

Σ n

i=1

(Vifixw -Vi(1-fi)xv) (12)

其中:xw =UT

iĨxk

σi

;xv =VT αk-αk, ( ) * ;fi =

σ2i

σ2i

( +λ ) 2 (i=1,2,…,ne)被定义为过滤因子。

不适定方程解的不稳定根源在于病态矩阵Sk的小奇异值趋于零的特性。正则化方法相当于利用过滤因子对病态矩阵的奇异值进行了修正,保证了解的稳定性。正则化参数λ 可采用L-curve的方法[18]求得。

识别误差为

Eerror= αid -αreal

αreal

×100% (13)

其中:αid

为识别的单元刚度差异参数;αreal

为实际单

元刚度差异参数。

αi确定之后,根据式(2)可确定待识别结构的刚度。

2.4 损伤识别方法实施步骤

1)对待识别结构进行几何量测和物理性能实验,获取必要的几何尺寸及基本的原始数据,建立结构的参考模型。

2)对待识别结构进行动力实验,布设一定数量的加速度传感器,测量结构的加速度响应¨xm。

3)根据式(1),利用步骤1中的参考模型和步骤2中的实验荷载值,采用Newmark-β直接积分法计算出结构的加速度响应¨x。根据式(5)计算出振动响应对结构刚度差异差数的灵敏度,形成加速度灵敏度矩阵S。

4)根据式(3)计算步骤2中测量的加速度响应与步骤3中计算的加速度响应的差值Δ¨x。

5)根据式(12)求解出结构刚度差异参数α,再根据式(2)计算结构刚度变化ΔK,利用式K′=K+ΔK 得到修正后的刚度K′。

6)重复步骤3~5,直到前后两迭代步的相对误差满足式(7)为止。

7)根据式(8)计算出最终识别的单元刚度差异参数,得到结构的损伤位置和损伤程度。

3 藏式古建筑木结构的损伤识别

模拟图5为一典型藏式古建筑木结构多层梁柱排架。根据藏式古建筑结构特点,取其中一榀排架作为损伤识别的研究对象,其结构简化有限元模型如图6所示。其中:“①”表示节点编号;“1”表示单元编号。由于柱脚直接搁置在柱础内,柱础提供竖向支撑力和水平摩擦力,故将柱脚设为铰接;梁端搭在墙上,简化为竖向铰接;梁柱节点假设为刚接。柱截面为270 mm×270 mm,梁截面为200 mm×300mm。木材的密度为418kg/m3,弹性模量为7.08×103 MPa[16]。楼板和椽子简化为1 000kg的集中质量块施加在每个梁节点上。

>

>

为模拟古建筑结构构件物理参数离散性较大的特点,为每个结构单元的弹性模量添加一个服从正态分布的范围在[-5%E,+5%E]内的随机误差,其中,E 为实验测得的藏式古建筑所用木材的弹性模量。

结构的前5阶固有频率分别为2.77,3.54,4.43,9.60和10.10Hz。结构的阻尼假设为瑞利阻尼,前两阶振型的阻尼比取为0.05。如图7所示,采用幅值缩减10倍的东西向EL-Centro地震波沿水平方向对结构进行激振(时长为5s)。此时,结构的动力方程可表示为M¨x+Cx+Kx =-M¨xg (14)其中:¨xg为地面加速度。

采用Newmark时程分析方法计算结构的加速度响应,计算过程中积分步长取0.01s,即加速度时程的采样频率为100Hz。在节点8,10,14和18上布设水平方向和竖直方向的加速度传感器。

为考虑实际中测量噪声对识别结果的影响,在“测量”的加速度信号上加入人工模拟噪声,即

¨xm =¨x′

m +EPNnoiseσ(¨x′

m) (15)

其中:¨x′

m,¨xm

分别为测点添加噪声前、后的加速度信号;Ep为噪声水平;Nnoise

表示噪声服从标准正态分布,具有0均值和单位标准差;σ(¨x′m)为加速度响应时程的均方差。

笔者研究了测量信号无噪声干扰、低水平噪声干扰(含1%噪声)和高水平噪声干扰(含5%噪声)3种情况对识别结果的影响。共模拟了3种损伤状态,分别为:a.14号柱单元损伤(弹性模量降低15%E);b.18号梁单元损伤(弹性模量降低15%E);c.

15号柱单元和9号梁单元同时损伤(弹性模量均降低15%E)。考虑结构单元刚度差异参数的物理意义,将式(10)中的lup设为1,llow设为-1。

表1列出了3种损伤状态在不同噪声水平干扰条件下的结构整体刚度的识别误差。可以看出,随着噪声水平的提高,识别误差也随之增大。对于第1种损伤状态,在测试信号没有噪声干扰时,经过52次迭代,结果达到收敛。如图8所示,各个单元的识别误差几乎为0,结构离散的物理参数、损伤位置和损伤程度均可以被准确地识别出来。在测试信号含1%噪声干扰时经过7次迭代,结果达到收敛。识别结果如图9所示。在测试信号含5%噪声干扰时,经过7次迭代,结果达到收敛。此时识别误差达到1.31%。由图10可以看出,虽然测量噪声的存在对识别结果有一定的负面影响,但是识别精度仍满足工程要求。

>

>

>

损伤状态3是多个构件存在损伤的状态,在5%的噪声水平下,识别误差为1.19%,与状态1相比,识别误差并没有随着损伤构件的增多而增大。以上结果证明,该识别方法既能识别出单个构件的损伤,也能识别出多个构件的损伤。在测量噪声达到5%的情况下,识别结果仍然满足要求。

 

4 结 论

1)藏式古建筑木结构梁、柱构件的主要残损类型包括:梁体裂缝、梁体过大挠曲、柱脚腐朽、柱身扭转和柱身裂缝等。导致损伤的主要原因可以概括为:a.木材的材质缺陷及材性退化的不利影响;b.结构本身构造的影响;c.环境及外荷载作用的影响等。

2)基于振动响应灵敏度理论,采用修正正则化方法的损伤识别方法适用于构件之间物理参数离散性较大的古建筑木结构的损伤识别。该方法不仅可以识别出单一构件的损伤,还可以识别多个构件的损伤。

3)损伤识别误差会随着测量信号中噪声含量的增加而增大,损伤构件的增多并不一定会导致识别误差的增大。在噪声水平为5%时,损伤位置和损伤程度仍能被较好地识别出来,识别结果可以满足工程要求。

4)本研究仅做了藏式古建筑木结构损伤识别的数值模拟工作,若将本研究提出的方法付诸实用,理论和技术问题需要进一步探讨。

 

 

参  考  文  献

[1] 陆铮.光纤光栅用于木结构健康监测技术[J].住宅科技,2005,11:25-28.

[2] 史学涛.结构健康监测系统的研究[D].上海:同济大学,2006.

[3] 王路.专家盛赞保国寺大殿科技保护监测系统[N].宁波日报,2008-06-23(A08).

[4] Doebling S W,Peterson L D,Alvin K F.Estimationof reciprocal residual flexibility from experimental modal data[J].Journal of Aeronautics and Astronautics,1996,34:1678-1685.

[5] 欧进萍,肖仪清,黄虎杰,等.海洋平台结构实时安全监测系统[J].海洋工程,2001,19(2):1-6.

[6] 李宏男,李东升.土木工程结构安全性评估、健康监测及诊断述评[J].地震工程与工程振动,2002,22(3):82-90.

[7] 李惠,周文松,欧进萍,等.大型桥梁结构智能健康监测系统集成技术研究[J].土木工程学报,2006,39

[8] 刘宗政,陈恳,郭隆德,等.基于环境激励的桥梁模态参数识别[J].振动、测试与诊断,2010,30(3):300-303.

[9] 马乾瑛,王社良,朱军强.大跨空间结构智能监测优化设计及信号处理[J].振动、测试与诊断,2011,31(3):286-290.

[10]Lu Z R,Law S S.Features of dynamic response sensitivity and its application_______ in damage detection[J].Journal of Sound and Vibration,2007,303(1-2):305-329.

[11]徐伟华,刘济科,吕中荣.基于振动响应的弦结构损伤检测[J].振动与冲击,2009,28(6):29-31.

[12]徐伟华,吕中荣,刘济科.基于振动响应的杆结构损伤检测[J].固体力学学报,2010,31(1):48-51.

[13]吕中荣,罗绍湘,刘济科.利用响应灵敏度修正Gascogine天桥的有限元模型[J].中山大学学报:自然科学版,2006,45(3):13-16.

[14]周乾,闫维明,纪金豹.明清古建筑木结构典型抗震构造问题研究[J].文物保护与考古科学,2011,23(2):36-48.

[15]李铁英.应县木塔现状结构残损要点及机理分析[D].太原:太原理工大学,2004.

[16]李鹏.藏式古建筑木构架梁柱节点力学机理研究[D].北京:北京交通大学,2009.

[17]Tikhonov A N.Solution of incorrectly formulated problems and the regularization method [J].Soviet Mathematics Doklady,1963,4:1035-1038.

[18]Hansen P C.Analysis of discrete ill-posed problems by means of the L-curve[J].Society for Industrial and Applied Mathematics,1992,34(4):561-580.

木结构 关注微信
手机点击拨打:13282000900

24小时客服电话

设计建造:13282000900