开源CAE实战系列（十一）：Code_Aster应用实例之混凝土大坝的结构抗震分析

1结构抗震理论

结构抗震理论的发展经历了由简单到复杂、由静力分析到动力分析的过程。早期的结构抗震设计主要采用静力分析法，其核心思想是将地震作用等效为静力荷载进行计算，因此形成了基于承载力的抗震设计方法。随着对地震作用和结构响应认识的深入，反应谱分析法逐渐发展起来，它能够考虑结构的动力特性，并在一定程度上反映地震作用下结构的最大响应，因此推动了基于承载力和延性的抗震设计方法的发展。进一步地，动力分析法成为更精细的抗震分析手段，它可以直接研究结构在地震波作用下随时间变化的响应过程，并结合时程分析法和推覆分析法等具体方法，对结构的弹塑性变形、损伤演化和抗震性能进行更全面的评价。由此，抗震设计理念也从最初单纯基于承载力的设计，逐步发展到兼顾承载力与延性的设计，并最终向基于性能的抗震设计方向演进。

1.1 静力分析法

静力分析法，也称为“底部剪力法”，假设结构物与地震动具有相同的振动，把结构物在地面运动加速度a作用下产生的惯性力视作静力作用于结构物上的抗震计算。惯性力计算公式为

a为地震动最大水平加速度；g为重力加速度，G=mg，为建筑物的重量，m为建筑物的质量；k为地面运动加速度峰值与重力加速度的比值，称为地震系数，其数值与结构动力特性无关。

静力分析法忽略了结构的动力学特性这一重要因素，把地震动加速度看作是结构地震破坏的单一因素，因此有很大的局限性，常导致对结构抗震能力的错误判断，仅当结构物在振动时几乎不产生变形而被当做刚体的时候，静力分析法才比较适用，因此现在已较少使用。

1.2 反应谱分析法

反应谱分析法，通过反应谱来计算结构动力特性（自振周期、振型和阻尼）所对应的共振效应，但其计算公式仍保留了早期静力法理论的形式。地震时结构所受到的最大水平基底剪力，即总水平地震作用为

其中，k为地震系数，β(T)则是加速度反应谱与地震动最大加速度的比值，它表示地震时结构振动加速度的放大倍数。反应谱理论尽管考虑了结构的动力特性，然而在结构设计中，它仍然把地震惯性力作为静力来对待，所以它只能称为准动力理论。

1.3 时程分析法

时程分析法是利用实际地震地面运动加速度记录进行抗震分析与设计，能考虑结构的弹性和弹塑性性态，能反映地震动的三大要素（频谱、振幅和持时），因而在复杂工程结构抗震分析和设计中，得到了广泛的应用。

它具有如下特点：

（1）输入地震动参数要求采用符合场地情况，具有概率意义的加速度时程，对于复杂结构要求考虑地震动三个分量的时间过程及其空间相关性；

（2）要求采用的结构和构件的动力模型应接近实际情况，应考虑结构和构件的非线性恢复力特性；

（3）该方法能计算出结构反应的全过程，包括变形和能量损耗的积累。

时程分析法是先进的，但是在应用上还有很多困难和局限，尤其是工程设计时的应用，例如输入地震波的不确定性、结构性能的近似假定和模拟等，使分析结果的可信度受到限制；该方法需要专门的程序和应用知识，输入输出数据量大，计算技术复杂，一般需要专门的技术人员进行分析。因此，我国规范只要求少数重要超高或有薄弱位置的结构，采用时程分析方法进行多发地震下的补充计算或罕见地震作用下薄弱层的弹塑性变形验算。

2 模态分析理论

模态分析理论有三个基本假设：

• 时不变性假设：系统的动态特性不会随着时间的变化而相应的变化
• 线性假设：任意的激励引起的结构的响应等于各个激励分别作用于结构所引起响应的线性叠加
• 可观测性假设：全部的数据信息都可以观测得到，即确定系统的动态特性所需要的所有数据为已知

结构动力学基本方程

通常，结构自身的动力学特性与外载荷激励无关，令{F(t)}=0得到

而阻尼对于结构自身的动力学特性影响不大，工程上往往忽略阻尼项，由此得到

在有限元分析中，常用的阻尼模型主要包括粘性阻尼、材料阻尼和摩擦阻尼三类。

粘性阻尼是有限元法中最常使用的阻尼模型，特点是阻尼力的大小与结构运动速度成正比。在工程分析中，常采用线性阻尼模型表示粘性阻尼，典型形式为瑞利阻尼，即：

材料阻尼主要由材料内部能量耗散产生，取决于材料的本构关系。在这种模型中，阻尼应力通常与应变率成正比。例如，Voigt模型可以表示为：

材料阻尼中粘弹性阻尼模型目前在工程中被广泛使用。

摩擦阻尼则来源于结构或构件接触面之间的相对滑移和摩擦耗能。最简单的是库仑摩擦阻尼模型，该模型通常属于非线性阻尼模型，能够反映摩擦力与相对运动方向之间的关系。

3 大坝未蓄水时的模态分析

整体上可以分为以下步骤：

（1）读取网格，定义单元

• LIRE_MAILIAGE：读取网格，格式med
• AFFE_MODELE：定义有限元模型的单元类型，定义拱坝为3D模型

（2）定义材料，赋予材料

• DEFI_MATERIAU：定义拱坝相关的材料参数，包括瑞利阻尼系数，为后续抗震分析做准备
• AFFE_MATERIAU：将定义好的材料属性赋予对应的有限元模型

（3）边界条件

• AFFE_CHAR_MECA：定义拱坝与地基交界面处的边界条件，LIAISON_UNIF可以连接单个节点和一个节点组，从而通过单个节点控制节点组的自由度

load = AFFE_CHAR_MECA(DDL_IMPO=_F(DX=0.0, DY=0.0, DZ=0.0, GROUP_NO=('arch_bottom_1p', )), LIAISON_UNIF=_F(DDL=('DX', 'DY', 'DZ'), GROUP_NO=('arch_bottom', )), MODELE=model)

（4）计算质量、刚度矩阵

• ASSEMBLAGE：计算拱坝结构的质量矩阵、刚度矩阵和阻尼矩阵，并定义一个数据结构NUME_DDL，这里的阻尼矩阵定义是为了后续的地震分析做准备

ASSEMBLAGE(CHAM_MATER=fieldmat, CHARGE=(load, ), MATR_ASSE=(_F(MATRICE=CO('M'), OPTION='MASS_MECA'), _F(MATRICE=CO('K'), OPTION='RIGI_MECA'), _F(MATRICE=CO('C'), OPTION='AMOR_MECA')), MODELE=model, NUME_DDL=CO('NDDL'))

（5）模态计算

• CALC_MODES：计算20阶模态

（6）结果后处理及保存

• IMPR_RESU：输出计算结果，查看振型

4 大坝蓄水后的模态分析

蓄水后的分析步骤包括：

（1）矩阵投影

• PROJ_MATR_BASE / PROJ_BASE：将质量矩阵、刚度矩阵和阻尼矩阵投影到前面计算的未蓄水时的模态空间中，两个命令中，前者一次只能完成一个投影操作，后者可以同时进行多个矩阵投影操作

PROJ_BASE(BASE=modes, MATR_ASSE_GENE=(_F(MATRICE=CO('M1'), MATR_ASSE=M), _F(MATRICE=CO('K1'), MATR_ASSE=K), _F(MATRICE=CO('C1'), MATR_ASSE=C)), NUME_DDL_GENE=CO('NDDL1'))

NUME_DDL_GENE：基于已计算的模态空间生成新的方程编号，即新的模态特征值对应的方程从1~20进行编号，不再与节点和物理自由度有关，只与前面计算的特征模态有关。

（2）库水附加质量

• MACRO_MATR_AJOU：计算未蓄水大坝的库水附加质量，并给库水赋予材料属性，注意需要指定自由水面的压力值0

MACRO_MATR_AJOU(DDL_IMPO=_F(GROUP_NO=('water_surface', ), PRES_FLUIDE=0.0), FLUIDE=_F(RHO=1000.0, TOUT='OUI'), GROUP_MA_FLUIDE=('water', ), GROUP_MA_INTERF=('arch_face', ), MAILLAGE=mesh, MATR_MASS_AJOU=CO('WATER'), MODELISATION='3D', MODE_MECA=modes, NUME_DDL_GENE=NDDL1)

（3）质量矩阵合并

• COMB_MATR_ASSE：合并大坝结构质量矩阵和库水的附加质量矩阵

M1 = COMB_MATR_ASSE(reuse=M1, COMB_R=(_F(COEF_R=1.0, MATR_ASSE=WATER), _F(COEF_R=1.0, MATR_ASSE=M1)))

（4）模态计算

• CALC_MODES：运用原来的刚度矩阵和合并后的质量矩阵计算新的模态

（5）转换模态空间基

• REST_GENE_PHYS：用附加质量的方法计算的模态结果是基于合并后的模态空间坐标系，该命令将计算结果恢复为广义的物理空间坐标系，以便输出模态振型结果

modes1 = CALC_MODES(CALC_FREQ=_F(NMAX_FREQ=20), MATR_MASS=M1, MATR_RIGI=K1, OPTION='PLUS_PETITE')

restran = REST_GENE_PHYS(RESU_GENE=modes1)

（6）结果后处理及保存

• IMPR_RESU：输出模态振型结果

5 大坝抗震分析

（1）静模态计算

• MODE_STATIQUE：计算坝基位移约束相关的静态模态

modes0 = MODE_STATIQUE(MATR_MASS=M, MATR_RIGI=K, MODE_STAT=_F(AVEC_CMP=('DX', 'DY', 'DZ'), GROUP_NO=('arch_bottom_1p', )))

（2）计算地震谱

• CALC_CHAR_SEISME：基于上面的静模态计算模型，生成X/Y/Z三个方向地震载荷向量
• LIRE_FONCTION：以函数的方式读入地震载荷谱数据

X = CALC_CHAR_SEISME(DIRECTION=(1.0, 0.0, 0.0), GROUP_NO=('arch_bottom_1p', ), MATR_MASS=M, MODE_STAT=modes0)

ax = LIRE_FONCTION(INDIC_PARA=(1, 1), INDIC_RESU=(1, 2), INTERPOL=('LIN', ), NOM_PARA='INST', PROL_DROITE='CONSTANT', PROL_GAUCHE='CONSTANT', TYPE='FONCTION', UNITE=2)

（3）矩阵投影

• PROJ_BASE：将质量矩阵、刚度矩阵和阻尼矩阵投影到蓄水模态计算空间中，包括三个方向的地震载荷向量投影

PROJ_BASE(BASE=restran, MATR_ASSE_GENE=(_F(MATRICE=CO('M3'), MATR_ASSE=M), _F(MATRICE=CO('K3'), MATR_ASSE=K), _F(MATRICE=CO('C3'), MATR_ASSE=C)), NUME_DDL_GENE=CO('DDL3'), VECT_ASSE_GENE=(_F(TYPE_VECT='FORC', VECTEUR=CO('X1'), VECT_ASSE=X), _F(TYPE_VECT='FORC', VECTEUR=CO('Y1'), VECT_ASSE=Y), _F(TYPE_VECT='FORC', VECTEUR=CO('Z1'), VECT_ASSE=Z)))

（4）反应谱分析

• DYNA_VIBRA：计算基于蓄水大坝模态空间的瞬态结构响应：1.选择合适的时间步；2.对比积分方案（欧拉积分与龙格库塔积分）；3.施加地震载荷

resharm = DYNA_VIBRA(BASE_CALCUL='GENE', EXCIT=(_F(FONC_MULT=ax, VECT_ASSE_GENE=X1), _F(FONC_MULT=ay, VECT_ASSE_GENE=Y1), _F(FONC_MULT=az, VECT_ASSE_GENE=Z1)), INCREMENT=_F(INST_FIN=4.0, INST_INIT=0.0, PAS=0.002), MATR_AMOR=C3, MATR_MASS=M3, MATR_RIGI=K3, TYPE_CALCUL='TRAN')

（5）转换模态空间基

• REST_GENE_PHYS：用附加质量的方法计算的模态结果是基于合并后的模态空间坐标系，该命令将计算结果恢复为广义的物理空间坐标系，以便输出模态振型结果

（6）结果后处理及保存

• RECU_FONCTION：后处理，提取三个位移分量DX/DY/DZ，生成时间历程曲线
• IMPR_RESU：输出模态振型结果

DX = RECU_FONCTION(GROUP_NO=('arch_edge_1p', ), NOM_CHAM='DEPL', NOM_CMP='DX', RESULTAT=restran0)