基于超单元结合虚拟变形法的桥梁静力有限元模型修正方法

技术领域

本发明涉及一种静力信息的实际运营桥梁结构的有限元模型修正方法。

背景技术

桥梁结构是连接公路的重要枢纽，其结构的安全性对保障公路交通的安全与畅通举足 轻重。有效、准确地评定桥梁结构的承载能力是保障桥梁运营安全的重要手段，而桥梁承 载能力的评定往往依靠一个能反映桥梁实际运营状况的精确基准有限元模型，基于静力测 试信息的桥梁有限元模型修正则是建立该基准模型的有效方法。

基于静力信息的桥梁结构有限元模型修正方法具有如下优点：模型修正结果不受结构 阻尼及质量信息的影响、测试结果信噪比较高、容易构造模型修正的超定优化问题等。但 是，对于实际桥梁结构而言，其结构体型巨大、自由度数量趋于无限，而且结构有限元模 型修正的优化问题求解往往需要大量的迭代运算，因此，如何解决大型复杂桥梁有限元模 型修正的计算效率问题仍然是一个具有挑战性的课题。

发明内容

本发明目的是为了解决现有桥梁结构有限元模型修正效率低，无法满足大型复杂桥梁 有限元模型修正的问题，提供了一种基于超单元结合虚拟变形法的桥梁静力有限元模型修 正方法。

本发明所述基于超单元结合虚拟变形法的桥梁静力有限元模型修正方法，该方法的具 体过程为：

步骤一：利用超单元技术对桥梁结构初始有限元模型进行缩聚，获取缩聚后的桥梁结 构有限元模型；

步骤二：根据步骤一获取的桥梁结构有限元缩聚模型，结合虚拟变形法建立桥梁结构 有限元模型的代理模型；

步骤三：根据步骤二获取的桥梁结构有限元模型的代理模型，进行基于静力信息的桥 梁结构有限元模型修正。

本发明的优点：

第一，采用超单元技术能够缩减桥梁结构初始有限元模型的规模；

第二，利用虚拟变形法能够建立大型复杂桥梁结构有限元模型的代理模型；

第三，该模型能够求解基于静力信息的桥梁有限元模型修正优化问题的迭代运算。

本发明所述基于超单元结合虚拟变形法的桥梁静力有限元模型修正方法能够大幅提 高基于静力信息的大型复杂桥梁结构有限元模型修正的计算效率，效率可以提高10～20 倍，同时还能够保证有限元模型修正的计算精度，代理模型与计算模型相对计算误差小于 0.1％，本发明适用于解决基于静力信息的实际运营桥梁结构有限元模型修正问题。

附图说明

图1是本发明所述基于超单元结合虚拟变形法的桥梁静力有限元模型修正方法的流 程框图；

图2是目标函数的迭代收敛曲线，曲线1代表工况1，曲线2代表工况2，曲线3代 表工况3，曲线4代表工况4；

图3是修正参数的迭代收敛曲线；

图4是三跨连续桥梁的结构示意图；

图5是主梁横断面图；

图6是静载试验过程中加载车辆轮距的平面布置图；

图7是94^#-95^#桥墩对称加载布置图；

图8是具体实施方式三中步骤二一中所述的空间梁单元的变形形式，a表示轴向变形， b表示弯曲和竖向变形，c表示弯曲和竖向变形，d表示弯曲变形，e表示扭曲变形，f表 示弯曲变形。

具体实施方式

具体实施方式一：下面结合图1说明本实施方式，本实施方式所述基于超单元结合虚 拟变形法的桥梁静力有限元模型修正方法，该方法的具体过程为：

步骤一：利用超单元技术对桥梁结构初始有限元模型进行缩聚，获取缩聚后的桥梁结 构有限元模型；

步骤二：根据步骤一获取的桥梁结构有限元缩聚模型，结合虚拟变形法建立桥梁结构 有限元模型的代理模型；

步骤三：根据步骤二获取的桥梁结构有限元模型的代理模型，进行基于静力信息的桥 梁结构有限元模型修正。

本实施方式中，所述虚拟变形法是一种快速结构重分析方法，只需要利用结构的影响 矩阵和初始有限元模型的计算结果就能够快速地计算出参数改变后的结构响应，利用这一 思想可有效建立桥梁结构有限元静力分析的代理模型，从而提高提高桥梁有限元模型修正 的计算效率。但是，虚拟变形法中的影响矩阵的维数是制约有限模型修正计算效率的关键 问题，而该维数与桥梁结构初始有限元模型的有限单元数量相关，因此，采用超单元技术 对桥梁初始有限元模型进行缩聚，从而解决了虚拟变形法中影响矩阵的降维问题。因此， 提出基于超单元结合虚拟变形法的桥梁静力有限元模型修正方法具有重要意义。

具体实施方式二：下面结合图1说明本实施方式，本实施方式对实施方式一作进一 步说明，步骤一所述获取缩聚后的桥梁结构有限元模型的具体过程为：

步骤一一、设桥梁结构含有n个自由度，其特征方程表示为：

式中，λ_i和分别表示桥梁结构的特征值和特征向量，其中i＝1,2,3…n；K和M分别表示 刚度矩阵与质量矩阵；

步骤一二、设x同时表示位移和特征向量，将位移向量分为主自由度方向x_m和从自由 度方向x_s，其中m和s分别代表主自由度和从自由度，将刚度矩阵按照自由度的不同划分 成四个子矩阵K_mm、K_ms、K_sm和K_ss，将质量矩阵按照自由度的不同划分成四个子矩阵M_mm、 M_ms、M_sm和M_ss，则式(1)转化为：

$\left(\begin{array}{cc}{K}_{\mathrm{mm}}& K_{\mathrm{ms}}\\ K_{\mathrm{sm}}& K_{\mathrm{ss}}\end{array}\right)\left(\begin{array}{c}{x}_m\\ x_s\end{array}\right)=\lambda \left(\begin{array}{cc}{M}_{\mathrm{mm}}& M_{\mathrm{ms}}\\ M_{\mathrm{sm}}& M_{\mathrm{ss}}\end{array}\right)\left(\begin{array}{c}{x}_m\\ x_s\end{array}\right)=\left(\begin{array}{c}{f}_m\\ f_s\end{array}\right)---\left(2\right)$

式中，f代表力向量，f_m表示主自由度方向上的力向量；f_s表示从自由度方向上的力向量， λ表示特征值；

步骤一三、忽略从自由度方向上的惯性力，即令f_s＝0，将式(2)展开，由第二行的 等式获取坐标变化式：

x_s＝-K_ss^-1K_smx_m＝C_Gx_m   (3)

式中，C_G＝-K_ss^-1K_sm；

步骤一四、建立坐标转换公式：

$x=\left(\begin{array}{c}{x}_m\\ x_s\end{array}\right)=\left(\begin{array}{c}{I}_m\\ C_G\end{array}\right)x_m=T_G{x}_m---\left(4\right)$

式中，I_m为单位阵，T_G为转换矩阵；

步骤一五、根据步骤一四的坐标转换公式(4)计算超单元的刚度矩阵K_super和质量矩 阵M_super：

$K_{\mathrm{super}}=T_G^T{\mathrm{KT}}_G---\left(5\right)$

$M_{\mathrm{super}}=T_G^T{\mathrm{MT}}_G---\left(6\right).$

具体实施方式三：下面结合图1说明本实施方式，本实施方式对实施方式一作进 一步说明，步骤二所述建立桥梁结构有限元模型的代理模型的具体过程为：

步骤二一、对空间梁单元的刚度矩阵进行特征值分析，获取空间梁单元的六种变形形 式，在局部坐标系下，单元的应变转换矩阵、单元节点位移与单元广义应变之间的关系为：

$ϵ=\mathrm{Gu}=\left(\begin{array}{cccccccccccc}-\frac{1}{L}& 0& 0& 0& 0& 0& \frac{1}{L}& 0& 0& 0& 0& 0\\ 0& 0& \frac{6}{L^2}& 0& -\frac{3}{L}& 0& 0& 0& -\frac{6}{L^2}& 0& -\frac{3}{L}& 0\\ 0& -\frac{6}{L^2}& 0& 0& 0& -\frac{3}{L}& 0& \frac{6}{L^2}& 0& 0& 0& -\frac{3}{L}\\ 0& 0& 0& 0& -\frac{1}{L}& 0& 0& 0& 0& 0& \frac{1}{L}& 0\\ 0& 0& 0& -\frac{1}{L}& 0& 0& 0& 0& 0& \frac{1}{L}& 0& 0\\ 0& 0& 0& 0& 0& -\frac{1}{L}& 0& 0& 0& 0& 0& \frac{1}{L}\end{array}\right)\left(\begin{array}{c}{u}_1\\ v_1\\ {\omega }_1\\ {\theta }_{x1}\\ {\theta }_{y1}\\ {\theta }_{z1}\\ u_2\\ v_2\\ {\omega }_2\\ {\theta }_{x2}\\ {\theta }_{y2}\\ {\theta }_{z2}\end{array}\right)---\left(7\right)$

式中：ε表示单元的广义应变向量，形式为ε＝[ξ,χ^xz,χ^xy,κ^xz,γ,κ^xy]^T；

ξ表示单元沿x方向的轴向变形形式下的广义应变；

χ^xz表示单元在XOZ平面内的弯曲和竖向变形耦合的变形形式下的广义应变；

χ^xy表示单元在XOY平面内的弯曲和竖向变形耦合的变形形式下的广义应变；

κ^xz表示单元在XOZ平面内的弯曲变形形式下的广义应变；

γ表示单元绕X轴的扭转变形形式下的广义应变；

κ^xy表示单元在XOY平面内的弯曲变形形式下的广义应变；

u表示单元两个节点的六个自由度方向上的位移；

G表示单元广义应变ε与单元节点位移u之间的转换矩阵；

u代表单元节点位移，其中u₁、v₁、ω₁、θ_x1、θ_y1、θ_z1、u₂、v₂、ω₂、θ_x2、θ_y2、θ_z2分别代表单元节点位移沿各自由度方向上的位移分量；

步骤二二、对桥梁结构的一个单元施加单位广义应变，其等价形式是在单元上施加一 对大小相等、方向相反的自平衡力，由公式(7)计算出产生单位应变后节点位移的大小， 从而利用单元刚度矩阵获取各自由度上施加的力的大小，在这一对平衡力作用下计算出桥 梁结构各单元的广义应变和节点位移，分别作为应变影响矩阵和位移影响矩阵的一列；对 于由n个单元、m个节点组成的桥梁结构体系，依次对结构进行3n次有限元分析，得到其 广义应变影响矩阵D的维数为3n×3n，位移影响矩阵B的维数为3m×3n；

步骤二三、虚拟变形法计算过程中涉及到两种结构分别是修正结构和变形机构，虚拟 变形法的前提是假设修正结构与受外荷载作用的变形结构的变形和内力大小完全一致，若 在单元i上施加虚拟变形则前提假设表示为：

${ϵ}_i={ϵ}_i^L+{ϵ}_i^R---\left(8\right)$

$\overline{A_i}{\sigma }_i=A_i({\sigma }_i^L+{\sigma }_i^R)---\left(9\right)$

式中：ε_i、σ_i分别表示修正结构受外荷载作用下的应变和应力；

分别表示变形结构受外荷载作用下的应变和应力；

分别表示变形结构受虚拟变形作用下的应变和应力；

A_i分别表示修正结构和变形结构的横截面面积；

步骤二四、单元i在虚拟变形作用下的应变和应力为：

${ϵ}_i^R=D_{\mathrm{ii}}{ϵ}_i^0---\left(10\right)$

${\sigma }_i^R=E_i(D_{\mathrm{ii}}-1){ϵ}_i^0---\left(11\right)$

式中，D_ii表示应变影响矩阵，E_i为第i个单元的弹性模量；表示虚拟变形向量ε⁰的第i个元素；

步骤二五、假设截面参数发生变化，而结构的弹性模量保持不变，然后将式(10)和 式(11)代入式(8)和式(9)中，得到：

$\overline{A_i}{ϵ}_i=A_i({ϵ}_i-{ϵ}_i^0)---\left(12\right)$

步骤二六、采用杆单元的刚度变化比值或者材料弹性模量的变化比值的效果相同，两 者都能够作为结构的截面特性参数，这一比值称为截面特性比例因子，表示为：

${\mu }_i=\frac{{ϵ}_i-{ϵ}_i^0}{{ϵ}_i}=\frac{\overline{A_i}}{A_i}---\left(13\right)$

步骤二七、利用结构参数改变前后应变相同的条件，将式(10)和式(13)代入式(8) 中，得到：

$[1-(1-{\mu }_i)D_{\mathrm{ii}}]{ϵ}_i^0=(1-{\mu }_i){ϵ}_i^L---\left(14\right)$

式中：δ_ij表示克罗内克符号；μ_i表示截面特性比例因子向量μ的第i个元素；表示虚拟 变形向量ε⁰的第j个元素；

步骤二八、将一个单元产生虚拟变形的情况推广到n个单元，即在结构上施加虚拟变 形的单元有n个，则得到下式以进行虚拟应变ε⁰的求解：

$\underset{j=1}{\overset{n}{\Sigma }}[{\delta }_{\mathrm{ij}}-(1-{\mu }_i)D_{\mathrm{ij}}]{ϵ}_j^0=(1-{\mu }_i){ϵ}_i^L---\left(15\right)$

步骤二九、利用变形结构的有限元模型，得到变形结构在外荷载作用下的节点位移u^L， 获取计算结构位移的代理模型为：

u＝u^L+Bε⁰   (16)

式中，B为位移影响矩阵。

本实施方式中，步骤二三中所述的修正结构是指结构截面特性发生改变的结构； 变形结构是指截面特性未发生任何改变的结构。

具体实施方式四：下面结合图1说明本实施方式，本实施方式对实施方式一作进一步 说明，步骤三所述进行基于静力信息的桥梁结构有限元模型修正的具体过程为：

步骤三一、根据步骤二获取的代理模型，通过改变修正参数的大小计算桥梁结构的位 移，由此构造桥梁结构位移残差向量得到模型修正的目标函数：

J＝(u^m-u^a)^TW(u^m-u^a)   (17)

式中：u^m表示结构的测试位移；u^a表示结构的理论计算位移，u^a＝u^L+Bε⁰；u^L表示变形 结构在外荷载作用下利用有限元模型计算出的位移；W表示模型修正过程中的权子矩阵；

步骤三二、利用基于复摄动理论的灵敏度计算方法，将位移对修正参数的一阶导数确 定在复数域，通过取得计算位移虚部的办法得到结构的位移灵敏度；

将位移函数u(x)在x＝μ+i△μ处进行泰勒级数的展开，得到，

$u(\mu +\mathrm{i\Delta \mu })=u\left(\mu \right)+u^{\prime }\left(\mu \right)\mathrm{\Delta \mu }·i-\frac{u^{\prime \prime }\left(\mu \right)}{2!}{\left(\mathrm{\Delta \mu }\right)}^2-\frac{u^{\prime \prime \prime }\left(\mu \right)}{3!}{\left(\mathrm{\Delta \mu }\right)}^3·i+...---\left(18\right)$

式中：μ表示结构有限元模型修正参数；△μ表示结构有限元模型修正参数的摄动量；u′(μ) 表示结构位移对修正参数的一阶导数，即是位移灵敏度；u″(μ)表示结构位移对修正参数的 二阶导数；u″′(μ)表示结构位移对修正参数的三阶导数；

然后将式(18)的实部和虚部提取出来，得到：

R(u(μ+i△μ))＝u(μ)+O((△μ)²)   (19)

I(u(μ+i△μ))＝u′(μ)△μ+O((△μ)³)   (20)

式中：R()和I()代表实部与虚部，O()为高阶无穷小；

再将式(20)的两边同时除以摄动量，将一阶导数分离出来，即为结构的位移灵敏度，

$u^{\prime }\left(\mu \right)=\frac{I\left(u(\mu +\mathrm{i\Delta \mu })\right)}{\mathrm{\Delta \mu }}+O\left({\left(\mathrm{\Delta \mu }\right)}^2\right)---\left(21\right)$

在对摄动量取值足够小的情况下，认为该误差忽略不计，于是结构位移灵敏度的计算 公式为：

$u^{\prime }\left(\mu \right)=\frac{I\left(u(\mu +\mathrm{i\Delta \mu })\right)}{\mathrm{\Delta \mu }}---\left(22\right)$

步骤三三、对每一个参数增加一个复摄动量，计算其对应的灵敏度，得到位移灵敏度 矩阵，计算迭代步长，

△μ＝(S^TWS)^-1S^TW(u^m-u^a)    (23)

式中：G表示结构位移灵敏度矩阵；

步骤三四、确定迭代收敛的条件，即完成有限元模型修正的优化求解。

下面结合下述试验验证本发明的效果：

以图4所示的三跨连续梁桥结构为例，图4为一个三跨连续梁桥结构，其主体结构为 混凝土材料，横截面具体尺寸如图5所示，该桥梁结构的静力加载方案如图6及图7所示。 为了便于进行结构有限元模型修正，将该模型共划分为260个等尺寸的有限单元，且选择 9个测点，分别是桥梁各跨的四分之一截面、跨中截面和四分之三截面。选用的修正参数 为绕z轴的截面惯性矩比例因子μ₁、μ₂、μ₃、μ₄、μ₅和μ₆，分别表示1-40单元、41-80 单元、81-130单元、131-180单元、181-220单元、221-260单元的截面惯性矩比例因子， 即认为由于损伤程度不同的原因，每一跨的左右半边截面特性变化不同。

运用超单元方法对初始有限元模型进行缩聚，通过计算分别利用2个、4个、5个和 10个单元作为一个超单元的结构静力位移和频率计算精度，得到频率计算相对误差如表1 所示，静力位移的相对误差几乎可以忽略不计，于是建立8个或3个单元缩聚为一个单元 的两种超单元，将结构有限元模型的单元数降为40个。

表1 缩聚模型的频率误差

记有限元模型计算位移为u^FEM，由超单元结合虚拟变形法所构建代理模型的计算位移 为u^VDM，结构位移计算如表2所示，说明了代理模型用于结构静力位移的计算是可行的。

表2 三种方法的计算位移

确定修正参数的一个初值μ＝[0.9,0.85,0.82,0.86,0.91,0.95]，通过有限元模型计算出此 时的桥梁控制截面的位移大小。通过迭代优化过程，得到目标函数的收敛结果如图2所示， 修正参数的收敛结果如图3所示，修正参数的具体收敛数值如表3所示，为了更好地反映 有限元模型修正的效果，将修正前后结构位移对试验测试位移的相对误差列于表4，测试 位移与修正前的位移的最大相对误差为15.074％，经过有限元模型修正之后，测试位移与 修正后的位移的最大相对误差为0.116％，说明修正过程取得较好的效果，能够让修正后的 有限元模型更接近于实际的结构。

表3 工况1有限元模型修正前、后位移比较

表4 模型修正参数修正前、后比较

本发明提出的基于超单元结合虚拟变形法的桥梁静力有限元模型修正方法完成本实施 例的耗时为15.902秒，而运用常规有限元模型进行模型修正的方法则耗时362.751秒，从 而证明了本发明所提方法大大提高了计算效率。