CN104809356A - 一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法 - Google Patents

Abstract

本发明公开了一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，本发明采用原始不可压缩粘性流体N-S方程描述水体的运动，采用Level Set方法精确追踪大幅晃动的界面，在槽身与晃动水体之间建立合适的耦合条件，通过交错迭代法求解系统的流固耦合方程，从而建立了一种考虑三维水体大幅晃动的大型渡槽流固耦合地震反应分析方法。本发明抛弃无旋、无粘理想流体假定，而采用原始的不可压缩粘性N-S方程描述流体的运动，采用精度高、适应性强的Level Set方法捕捉自由表面，真实模拟大型渡槽遭受强烈地震时的非线性大幅度晃动。

Description

一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法

技术领域

本发明涉及渡槽结构抗震分析技术领域，特别是一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法。

背景技术

作为常见的输水建筑物，大型渡槽在长距离调水、水资源优化配置等方面发挥着巨大的作用，其重要性不言而喻。渡槽结构是一种特殊的薄壁输水建筑物，其上部水体的质量往往大于槽体本身的质量，这对于渡槽的抗震是极为不利的。渡槽遭受强烈地震时，结构的振动引发流体的运动，反过来流体的运动又改变着结构的振动，属于典型的流固耦合现象，流体的作用是不可回避且必须加以考虑的问题，是渡槽结构抗震研究的技术关键之一。

目前，渡槽结构的流固耦合分析方法主要有：Westergaard附加质量法、HOUSNER模型、势流体有限元流固耦合模型、边界元法、ALE有限元法等几类。附加质量法显然是极其简化的做法，既不能反映流固耦合作用，也不能反应槽内水体的晃动；HOUSNER模型简单易行，采用弹簧质量系统来近似计算复杂的运动液体对固体产生的脉冲压力和对流压力，但未能真正反映水体的晃动以及耦合机理；势流体有限元流固耦合模型考虑了流固耦合效果，但由于采用了无旋、无粘的理想流体假定，实际上用线性的方法解决非线性问题，忽略了水体大幅晃动这一重要影响因素；边界元法能够较好的模拟水体的晃动现象，但较难推广到三维情况；ALE有限元法采用原始的非线性流体控制方程，对于流体大幅晃动及其流固耦合问题取得了较好的效果，只是自由面的存在使得对网格更新提出了更高的要求，在自由面大幅度剧烈晃动的情况容易产生网格畸变而导致求解失败，对水波发生重叠或破碎的情况更是无能为力。

综合目前的研究，大型渡槽流固耦合地震反应分析实际上有三个方面的问题有待于进一步完善：1、应采用原始非线性流体控制方程描述水体的运动。渡槽内水体与坝前库水大范围水体相比，严格来说已不再适用无旋、无粘的理想流体的假定；2、大幅晃动问题的准确描述。已有研究结果表明，渡槽结构中水体的晃动具有明显的非线性特征，用线性化的方法描述强烈非线性的大幅晃动显然是不合适的；3、流固耦合系统的求解。单纯的自由面晃动问题已有不少丰富的研究成果，但考虑大幅度晃动下的流固耦合问题的求解则需进一步丰富。

发明内容

本发明的所要解决的技术问题是克服现有技术的不足而提供一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，采用原始的不可压缩粘性流体N-S方程描述水体的运动，采用Level Set方法精确追踪大幅晃动的界面，在槽身与晃动水体之间建立合适的耦合条件，通过交错迭代法求解系统方程。

本发明为解决上述技术问题采用以下技术方案：

根据本发明提出的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，包括如下步骤：

步骤1、建立固体结构动力学方程；

步骤2、采用原始的不可压缩粘性流体，建立流体N-S方程；

步骤3、采用Level Set方法追踪大幅晃动的界面，建立Level Set函数控制方程；

步骤4、建立流固耦合面Γ_fs上控制方程；

步骤5、采用交错迭代法求解步骤1、步骤2、步骤3和步骤4中所建立的方程。

作为本发明所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法进一步优化方案，所述步骤1中固体结构动力学方程为：

$M_s\stackrel{\·\·}{\mathrm{\δ}}+C\stackrel{\·}{\mathrm{\δ}}+\mathrm{K\δ}=F_e+F_f;$

其中，M_s是质量矩阵，C是阻尼矩阵，K为总体刚度矩阵，C＝αM_s+βK，α是第一瑞利比例阻尼系数，β为第二瑞利比例阻尼系数，δ为位移列阵，为速度列阵，为加速度列阵，F_e为外荷载向量，F_f为流体作用于固体的荷载。

作为本发明所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法进一步优化方案，所述步骤2中流体N-S方程为：

$\▿\·u=0;$

$\frac{\∂u}{\∂t}+\▿\·\left(\mathrm{uu}\right)=\▿\·\frac{1}{\mathrm{\ρ}}\mathrm{\τ}-\frac{1}{\mathrm{\ρ}}\▿p+f_b;$

其中，ρ为密度，τ为粘性应力张量，μ为动力粘性系数，T为转置符号，f_b为体力，g为重力加速度，为外加激励加速度向量，t为时间，u为流场速度，p为流场压力。

作为本发明所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法进一步优化方案，所述步骤3中Level Set函数控制方程为：

其中，为距离函数且用来指示不同的流体区域，定义为某空间点到交界面的符号距离，在交界面以上为正，以下为负，而零值处代表了交界面；

其中，S( )是符号函数，是方程在tⁿ时刻的解，初始条件为n表示第n时间步。

作为本发明所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法进一步优化方案，所述步骤4中流固耦合面上的控制方程为：

$\left.\begin{array}{c}\stackrel{\‾}{\stackrel{\·}{\mathrm{\δ}}}=\stackrel{\‾}{u}\\ \stackrel{\‾}{t_s}=\stackrel{\‾}{t_f}\end{array}\right\}$ 在Γ_fs上

其中，Γ_fs为流固耦合面，为流固耦合面上固体的速度，为流固耦合面上流体的速度，为流固耦合面上固体的表面应力，为流固耦合面上流体的表面应力，根据得到固体动力学方程中F_f的表达式：其中，ds为流固耦合面上微元的面积。

作为本发明所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法进一步优化方案，所述步骤5具体如下：

(1)考虑流体作用于固体的荷载F_f，求解步骤1中已建立的固体结构动力学方程，得到位移δ、速度响应量；

(2)以流固耦合面Γ_fs上的固体运动速度为边界条件，计算步骤2中已建立的N-S方程，得到流场速度u，流场压力p；

(3)计算步骤3中已建立的Level Set函数控制方程，捕捉自由表面；

(4)根据所述的初始条件对步骤3中已建立的Level Set函数控制方程进行重新初始化；

(5)根据所述的计算流体作用于固体的荷载F_f，返回(1)进行迭代计算，直到收敛后进入下一时间步的计算。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：抛弃无旋、无粘理想流体假定，而采用原始的不可压缩粘性N-S方程描述流体的运动，采用精度高、适应性强的LEVEL SET方法捕捉自由表面，真实模拟大型渡槽遭受强烈地震时的非线性大幅度晃动。

附图说明

图1是渡槽槽身截面几何尺寸图。

图2是渡槽槽身有限元模型。

图3a是跨中截面时间为t＝0s、t＝1s、t＝3s时渡槽内水体自由面位置。

图3b是跨中截面时间为t＝5s、t＝8s、t＝9s时渡槽内水体自由面位置。

图3c是跨中截面时间为t＝10s、t＝11s、t＝12s时渡槽内水体自由面位置。

图4是左右两侧槽壁波高时程曲线。

图5是不同计算模型得到的槽顶横向水平位移时程曲线。

图中的附图标记解释为：1—支撑结构，2—槽体，3—水体，4—空气。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

本发明考虑三维水体大幅晃动的大型渡槽流固耦合地震反应分析方法，包括如下步骤：

步骤1、建立固体结构动力学方程，其方程为：

$M_s\stackrel{\·\·}{\mathrm{\δ}}+C\stackrel{\·}{\mathrm{\δ}}+\mathrm{K\δ}=F_e+F_f;$

其中，M_s是质量矩阵，C是阻尼矩阵，K为总体刚度矩阵，C＝αM_s+βK，α是第一瑞利比例阻尼系数，β为第二瑞利比例阻尼系数，为两个待定常数，δ为位移列阵，为速度列阵，为加速度列阵，F_e为外荷载向量，F_f为流体作用于固体的荷载。

步骤2、采用原始的不可压缩粘性流体，建立流体N-S方程，其方程为：

$\▿\·u=0$

$\frac{\∂u}{\∂t}+\▿\·\left(\mathrm{uu}\right)=\▿\·\frac{1}{\mathrm{\ρ}}\mathrm{\τ}-\frac{1}{\mathrm{\ρ}}\▿p+f_b;$

其中，ρ为密度，τ为粘性应力张量，μ为动力粘性系数，T为转置符号，f_b为体力，g为重力加速度，为外加激励加速度向量，t为时间，u为流场速度，p为流场压力。

步骤3、采用Level Set方法追踪大幅晃动的界面，Level Set函数控制方程为：

其中，为距离函数且用来指示不同的流体区域，定义为某空间点到交界面的符号距离，在交界面以上为正，以下为负，而零值处代表了交界面；

为了保证为距离函数，要使满足方程：

其中，S( )是符号函数，是方程在tⁿ时刻的解，初始条件为n表示第n时间步。

步骤4、建立流固耦合面上控制方程，其控制方程为：

$\left.\begin{array}{c}\stackrel{\‾}{\stackrel{\·}{\mathrm{\δ}}}=\stackrel{\‾}{u}\\ \stackrel{\‾}{t_s}=\stackrel{\‾}{t_f}\end{array}\right\}$ 在Γ_fs上

其中，Γ_fs为流固耦合面，为流固耦合面上固体的速度，为流固耦合面上流体的速度，为流固耦合面上固体的表面应力，为流固耦合面上流体的表面应力，根据得到固体动力学方程中F_f的表达式：其中，ds为流固耦合面上微元的面积。

步骤5、采用交错迭代法求解系统的流固耦合方程，具体步骤为：

(1)考虑流体作用于固体的荷载F_f，求解步骤1中已建立的固体结构动力学方程，得到位移δ、速度响应量；

(2)以流固耦合面Γ_fs上的固体运动速度为边界条件，计算步骤2中已建立的N-S方程，得到流场速度u，流场压力p；

(3)计算步骤3中已建立的Level Set函数控制方程，捕捉自由表面；

(4)根据所述的初始条件对步骤3中已建立的Level Set函数控制方程进行重新初始化；

(5)根据所述的计算流体作用于固体的荷载F_f，返回(1)进行迭代计算，直到收敛后进入下一时间步的计算。

下面以一个具体实例来说明本发明：

某大型U型渡槽单跨长17.3m，槽体2两端搁置在支撑结构1上，不考虑支撑结构1的变形。渡槽结构截面几何尺寸见图1，其三维有限元模型见图2，共划分为固体单元16736个，水体3单元31159个和空气4单元23800个。槽体1本身为C40混凝土，弹性模量为32.5GPa,泊松比为0.167，容重取为25.0kN/m³，槽内水深按3.45m。地震峰值加速度0.3g，特征周期0.25s，考虑水平横向人工加速度地震波的作用。

图3a，图3b，图3c分别给出了用本发明计算得到的跨中截面时间为t＝0s、t＝1s、t＝3s、t＝5s、t＝8s、t＝9s、t＝10s、t＝11s、t＝12s时渡槽内水体自由面位置；图3a是跨中截面时间为t＝0s、t＝1s、t＝3s时渡槽内水体自由面位置，图3b是跨中截面时间为t＝5s、t＝8s、t＝9s时渡槽内水体自由面位置，图3c跨中截面时间为t＝10s、t＝11s、t＝12s时渡槽内水体自由面位置。

图4给出了左右两侧槽壁波高时程曲线，从中可以看出，在地震激励作用下，一开始自由面微幅晃动，波峰和波谷峰值基本相等。随着时间的增长而振幅逐渐增加，且波峰幅值大于波谷幅值，自由面最大晃动高度在0.58～0.64m范围内，随时间增长表现出强非线性特性。在地震结束之后，槽内水体呈自由晃动，晃动圆频率为2.09rad/s，与理论值2.05rad/s非常吻和，这也从一个方面说明了本发明方法的正确性。此外，对比在地震激励下水体在渡槽的动力作用下晃动频率可以发现，槽内水体的晃动频率基本不会随地震激励的频率发生变化。

图5给出了不同计算模型得到的槽顶横向水平位移时程曲线。从中可以看出，槽内无水时当地震激励结束后，槽顶位移幅值迅速减小为零，整个过程中最大幅值为0.77cm。考虑流固耦合以及大幅晃动效应后，结构上部集中质量增加，振幅平均增加0.96cm，增幅为24.7％。在8s后，两条曲线随时间变化的规律明显不同，槽内无的模型槽体结构经过短暂的小幅振动趋于静止；而流固耦合模型在地震结束后槽顶位移呈周期性振动，振动幅度约为0.13cm。

综上所述，本发明抛弃无旋、无粘理想流体假定，而采用原始的不可压缩粘性N-S方程描述流体的运动，采用精度高、适应性强的LEVEL SET方法捕捉自由表面，能够真实模拟大型渡槽遭受强烈地震时的非线性大幅度晃动。

Claims (6)

1.一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，其特征在于，包括如下步骤：

步骤1、建立固体结构动力学方程；

步骤2、采用原始的不可压缩粘性流体，建立流体N-S方程；

步骤3、采用Level Set方法追踪大幅晃动的界面，建立Level Set函数控制方程；

步骤4、建立流固耦合面Γ_fs上控制方程；

步骤5、采用交错迭代法求解步骤1、步骤2、步骤3和步骤4中所建立的方程。

2.根据权利要求1所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，其特征在于，所述步骤1中固体结构动力学方程为：

$M_s\stackrel{..}{\mathrm{\δ}}+C\stackrel{.}{\mathrm{\δ}}+\mathrm{K\δ}=F_e+F_f;$

其中，M_s是质量矩阵，C是阻尼矩阵，K为总体刚度矩阵，C＝αM_s+βK，α是第一瑞利比例阻尼系数，β为第二瑞利比例阻尼系数，δ为位移列阵，为速度列阵，为加速度列阵，F_e为外荷载向量，F_f为流体作用于固体的荷载。

3.根据权利要求2所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，其特征在于，所述步骤2中流体N-S方程为：

$\▿\·u=0;$

$\frac{\∂u}{\∂t}+\▿\·\left(\mathrm{uu}\right)=\▿\·\frac{1}{\mathrm{\ρ}}\mathrm{\τ}-\frac{1}{\mathrm{\ρ}}\▿p+f_b;$

其中，ρ为密度，τ为粘性应力张量，μ为动力粘性系数，T为转置符号，f_b为体力，g为重力加速度，为外加激励加速度向量，t为时间，u为流场速度，p为流场压力。

4.根据权利要求3所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，其特征在于，所述步骤3中Level Set函数控制方程为：

其中，为距离函数且用来指示不同的流体区域，定义为某空间点到交界面的符号距离，在交界面以上为正，以下为负，而零值处代表了交界面；

其中，S()是符号函数，是方程在tⁿ时刻的解，初始条件为n表示第n时间步。

5.根据权利要求4所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，其特征在于，所述步骤4中流固耦合面上的控制方程为：

$\left.\begin{array}{c}\stackrel{\stackrel{\‾}{.}}{\mathrm{\δ}}=\stackrel{\‾}{u}\\ \stackrel{\‾}{t_s}=\stackrel{\‾}{t_f}\end{array}\right\}$ 在Γ_fs上

其中，Γ_fs为流固耦合面，为流固耦合面上固体的速度，为流固耦合面上流体的速度，为流固耦合面上固体的表面应力，为流固耦合面上流体的表面应力，根据得到固体动力学方程中F_f的表达式：其中，ds为流固耦合面上微元的面积。

6.根据权利要求5所述的一种考虑水体大幅晃动的渡槽流固耦合地震反应分析方法，其特征在于，所述步骤五具体如下：

(1)考虑流体作用于固体的荷载F_f，求解步骤1中已建立的固体结构动力学方程，得到位移δ、速度响应量；

(2)以流固耦合面Γ_fs上的固体运动速度为边界条件，计算步骤2中已建立的N-S方程，得到流场速度u，流场压力p；

(3)计算步骤3中已建立的Level Set函数控制方程，捕捉自由表面；

(4)根据所述的初始条件对步骤3中已建立的Level Set函数控制方程进行重新初始化；

(5)根据所述的计算流体作用于固体的荷载F_f，返回(1)进行迭代计算，直到收敛后进入下一时间步的计算。