CN102306220A - 一种基于松耦合技术的索膜结构抗风设计方法 - Google Patents

Abstract

本发明提出一种基于松耦合技术的索膜结构抗风设计方法，包括下列步骤：步骤一：依据实际材料的物理性质和结构的初始预应力，通过找形方法建立索膜结构有限元模型；步骤二：根据步骤一中的找形结果和建筑周围环境，建立索膜结构流体计算模型；步骤三：利用风洞试验或实测数据或理论模拟的风速时程数据，作为流体计算的入口边界条件；步骤四：利用流固耦合协同平台进行膜结构抗风设计。本发明解决了索膜结构流固耦合传统计算中的收敛性问题，提高了计算速度。不但计算准确性相比以往简化方法和单边耦合方法有了显著的提高，而且计算速度非常快，使这项技术从科学研究走向实际应用。

Description

一种基于松耦合技术的索膜结构抗风设计方法

技术领域

本发明属于建筑行业索膜结构抗风分析与设计领域，且特别涉及一种基于松耦合技术的索膜结构抗风设计方法。

背景技术

流固耦合现象在斜拉桥、超高层建筑、塔桅结构、输电线、索膜结构等土木工程领域扮演着相当重要的角色(M.Gluck，M.Breuer，F.Durst，A.Halfmann，E.Rank.Computation of wind-induced vibrations of flexible shells and membranousstructures[J].Journal of Fluids and Structures 17(2003)739-765.)。过去由于缺乏对流固耦合现象的认识，曾发生过1879年苏格兰Tay大桥的垮塌，1948年西雅图Tacoma大桥风振倒塌，以及1965英格兰三座冷却塔倒塌等重大灾难(M.Gluck，M.Breuer，F.Durst，A.Halfmann，E.Rank.Computation of fluid-structureinteraction on lightweight structures[J].Journal of Wind Engineering and IndustrialAerodynamics 89(2001)1351-1368.)。如今，一方面诸如弛振、颤振等风致结构气动现象已为人们所熟知。另一方面由于大量轻盈结构的使用(例如亚特兰大七万人体育场屋顶采用膜结构，同时结构跨度达到240米)，建筑对风的响应越来越敏感。为了准确获得具有外形复杂或是受外界风场影响显著的建筑结构的风荷载，工程界已经进行了大量费用高昂的风洞实验或是采用半经验方法对流固耦合问题进行分析计算。

此外，工程界已经可以通过使用相对成熟的有限元软件对结构在风荷载下的应力、变形等效应进行精确计算。同时，也能够通过使用计算流体动力学软件对结构所受的风荷载进行准确估计。然而，对于进行流固耦合计算的方法和工程应用却少之又少。其主要原因有：(1)采用简单的松耦合方法会在膜结构上产生过大的加速度而导致求解不稳定(野村卓志·神明正人，基于ALE有限要素法的膜结构流固耦合分析，应用力学论文集.Vol.1，pp：241-251，1998)。(2)复杂的插值方法计算量巨大，缺乏工程可操作性(孙晓颖，薄膜结构风振响应中的流固耦合效应研究[D]，哈尔滨工业大学博士学位论文，2007年1月)。(3)不同软件间缺乏有效的消息传递与统一的数据传递接口。

本方法针对以上问题，提出了荷载分步修正技术，基于三角形重心坐标系的网格变形插值方法，以及利用第三方文件实现信息传递等技术，使基于松耦合的索膜结构抗风设计从理论变成现实。在此基础上，通过对世博轴索膜结构的实例分析，验证了上述方法的有效性和准确性。

发明内容

如背景技术所述，目前索膜结构工程无法使用先进的流固耦合技术，究其原因主要有：程序内部性能、数据外部处理与通信、计算效率等因素。因此，本发明提出一种基于松耦合技术的索膜结构抗风设计方法，其目的就是要通过解决以上问题，为索膜结构的抗风设计提供一种准确、高效的实用方法。

为了达到上述目的，本发明提出一种基于松耦合技术的索膜结构抗风设计方法，包括下列步骤：

步骤一：依据实际材料的物理性质和结构的初始预应力，通过找形方法建立索膜结构有限元模型；

步骤二：根据步骤一中的找形结果和建筑周围环境，建立索膜结构流体计算模型；

步骤三：利用风洞试验或实测数据或理论模拟的风速时程数据，作为流体计算的入口边界条件；

步骤四：利用流固耦合协同平台进行膜结构抗风设计。

进一步的，所述步骤一中找形方法为采用分步计数、逐步刚度逼近的方法建立索膜结构有限元模型。

进一步的，所述步骤三中数值风洞模型的建立仅对有限元模型所对应的膜面设置动网格属性，其余部分保持不动。

进一步的，所述步骤三中流体计算的非稳态计算采用LES模拟、隐式计算。

进一步的，所述有限元模型采用稀疏网格，所述流体计算模型采用稠密网格。

进一步的，所述步骤四中利用流固耦合协同平台进行膜结构抗风设计为考虑流固耦合作用进行紊流计算膜面的风振响应。

本发明通过：(1)在1个时间步上进行流体及结构间反复交替、不断修正求解，解决以往索膜结构流固耦合计算中的收敛性问题。(2)采用三角形重心坐标系进行节点定位，大幅提高数值计算的速度。(3)通过数据通讯管理和实时监控程序，实现索膜结构抗风自动化设计。本方法可以自动实现索膜结构抗设计的动力时程计算，可以考虑索膜结构风振响应中的气动阻尼、附加质量、气承刚度等机理问题，使计算结果更符合实际情况。目前，传统的流固耦合计算仅停留在小规模的科学研究，本方法提供的自动化技术手段和高效的算法，使这项技术在实际工程中的应用成为可能，并能改变过去只能采用简化手段进行索膜结构抗风设计的现象。

本发明解决了索膜结构流固耦合传统计算中的收敛性问题，提高了计算速度。不但计算准确性相比以往简化方法和单边耦合方法有了显著的提高，而且计算速度非常快，使这项技术从科学研究走向实际应用。

附图说明

图1所示为本发明较佳实施例的基于松耦合技术的索膜结构抗风设计方法流程图。

图2a所示为传统松耦合算法。

图2b所示为本发明荷载步修正算法。

图3所示为重心坐标计算原理图。

图4所示为变形跟踪系统实际使用效果。

图5所示为流固耦合技术方案流程。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

请参考图1，图1所示为本发明较佳实施例的基于松耦合技术的索膜结构抗风设计方法流程图。本发明提出一种基于松耦合技术的索膜结构抗风设计方法，包括下列步骤：

步骤一S100：依据实际材料的物理性质和结构的初始预应力，通过找形方法建立索膜结构有限元模型；

步骤二S200：根据步骤一中的找形结果和建筑周围环境，建立索膜结构流体计算模型；

步骤三S300：利用风洞试验或实测数据或理论模拟的风速时程数据，作为流体计算的入口边界条件；

步骤四S400：利用流固耦合协同平台进行膜结构抗风设计。

根据本发明较佳实施例，所述步骤一S100中找形方法为采用分步计数、逐步刚度逼近的方法建立索膜结构有限元模型，所述步骤三S300中数值风洞模型的建立仅对有限元模型所对应的膜面设置动网格属性，其余部分保持不动，所述步骤三S300中流体计算的非稳态计算采用LES模拟、隐式计算，所述有限元模型采用稀疏网格，所述流体计算模型采用稠密网格，所述步骤四S400中利用流固耦合协同平台进行膜结构抗风设计为考虑流固耦合作用进行紊流计算膜面的风振响应。

步骤四S400的协同平台中包含若干程序，这些程序因使用软件的不同而异，但其核心思路是相同的。

(1)收敛性问题的解决

传统的松耦合在一个时间步内仅需进行一个循环的流体计算及结构动力响应计算，这种方法对于像空气与建筑物这样相对质量比大的情况有效。但是，当由建筑物变形产生的附加空气质量不能忽略时，以上方法会由于在膜上生成过大加速度从而导致求解不稳定。图2a所示为传统松耦合算法。图中F_i为i时刻结构所受荷载，F_i+1为下一时刻i+1时刻需要加载的量。按照传统算法，及时荷载步被细分，如图中所示分为n个荷载步，每次增加dF/n的量，最终和一次加载的效果是一样的。图2b为改进方法，即荷载步被细分，但F_i+1也随时进行更新，从第二个荷载子步开始，每次增量dF/n变成dF’/n、dF”/n、dF”’/n......以此类推。因此，最终的加载效果如图2b中曲线所示。事实上，这种算法更符合流固耦合的实际发生机理，同时也使数值计算的运行性能得以改善，使计算过程更加容易收敛。

(2)运行性能的改善

虽然，目前的固体软件和流体软件均能在各自的领域做出精确的分析，但当考虑两者的相互耦合计算时却碰到了许多困难。其中之一便是：网格数量多的一方(流体网格)需要跟踪网格数量少的一方(固体网格)进行变形。以往的方法，如样条插值法、守恒法等，虽然保证了数据传递的准确性，但由于计算方法过于繁琐，导致计算速度成倍下降，使得这种先进的技术无法应用于实际工程，只能停留在小规模的科学研究或由于“计算量巨大，缺乏工程可操作性”，只好采用“简化的数值模拟方法”。本方法中对于网格插值问题的考虑巧妙地采用三角形重心坐标原理，使得原先位于一个稀疏网格(固体网格)三角形单元内部的稠密网格(流体网格)节点会随着稀疏网格(固体网格)节点的变形依然贴附于原三角形内。三角形内的点会随着三角形在平面内、外发生大变形，但这些点对所从属的整个三角形的拓扑关系不发生变化，如图2所示。三角形重心原理在本领域内的使用尚属首次，且真是由于该原理算法朴素、计算简单，很容易使本领域的普通技术人员能够掌握和实施。其原理如下：

图3中三角形ABC代表稀疏网格(固体网格)中的某一三角形单元，A、B、C分别为从属于该单元的节点，点P代表落在三角形ABC内部的稠密网格(流体网格)的一个节点。假设以点A作为起点，那么点B相当于在AB方向移动一段距离得到，而点C相当于在AC方向移动一段距离得到。则对于平面内任意一点，可由如下方程来表示：

P＝A+u×(C-A)+v×(B-A)                             (1)

整理方程(1)得到

P-A＝u×(C-A)+v×(B-A)                             (2)

令v₀＝C-A，V₁＝B-A，v₂＝P-A则，

v₂＝u×v₀+v×v₁                                    (3)

等式(3)两边分别点乘v₀和v₁得到两个等式：

(v₂)□v₀＝(u×v₀+v×v₁)□v₀                        (4)

(v₂)□v₁＝(u×v₀+v×v₁)□v₁                        (5)

求解上述方程(4)、(5)得到：

如果系数u或v为负值，那么相当于朝相反方向移动，即BA或CA。因此，如果：u≥0，v≥0，且u+v≤1则可以判断点P位于三角形ABC内部。

(3)用于外部数据处理与通信的措施

松耦合计算需要调用至少两个独立的程序，通常这两个程序往往由不同的公司开发，他们之间的信息互通以及数据传递需要编制专门的接口程序。其中最关键的要素是信息的互通，即当程序A运算完毕后需要暂停运行，同时通知程序B进行相应的计算，而当程序B完成计算后同样需要暂停运行等待程序A的运行。本发明中采用第三方通知文件的方法完成这一操作，即：当程序A结束一个时间步的运算后，先暂停计算，然后发消息给文件C，而程序B则始终保持与文件C的联系，一旦文件C得到程序A结束运算的信号，则马上告知程序B进行相应操作，反之依然，由此完成程序A与B之间的信息互通。图4所示为变形跟踪系统实际使用效果。图5所示为流固耦合技术方案流程。

作为2010年上海世博会立体交通组织重要载体的世博轴，由两层地下空间和一层高架步道组成。高架步道上覆连续张拉式索膜结构屋顶，总长约840m，最大跨度约97m，总面积近64000m²。包括支点系统和膜面系统两部分。

支点系统通过水平索将31组外桅杆及背索、19个下拉点，18个与阳光谷的拉结点连接成稳定体系。膜结构轻盈飘逸，膜面采用国家A级PTFE膜材，主要由脊索、边索、谷索和膜形成连续的三角倒锥形膜单元。

世博轴地处台风多发的上海地区(50年重现期基本风压0.55kN/m²)，项目现场位于黄浦江畔，周边地形开阔无高大建筑阻挡，属C类地貌。膜结构最大标高38m，平均标高超出10m平台约15m，属于低矮房屋。对于此类建筑，《荷载规范》既无体系系数可供参考，更无风振系数的计算方法。

形状确定是膜结构设计的第一步，也是决定一个膜结构设计是否成功的关键。相比于普通小品膜结构而言，世博轴索膜结构具有跨度更大、结构更柔的特点，采用有限元方法分析存在收敛困难等问题。此处采用分步计数、逐步刚度逼近的方法成功解决这一问题。

为考虑周边建筑对膜面风荷载的影响，数值风洞模型包含了阳光谷、10m平台等其他建筑。为提高计算效率，仅对有限元模型所对应的膜面设置动网格属性，其余部分保持“不动”。非稳态计算采用LES模拟、隐式计算。

为了研究不同网格体系的计算精度和效率，设计了以下三套网格工况：工况一：有限元模型和流体力学模型均采用具有3924个单元的一致网格体系；工况二：有限元模型和流体力学模型均采用具有11935个单元的一致网格体系；工况三：有限元模型采用具有3924个单元的稀疏网格，流体力学模型采用具有11935个单元的稠密网格，两者之间采用上文所述的变形协调系统传递数据。

CFD计算在来流为10m/s的均匀流条件下A点的位移。计算开始的初始1min结构并不参与计算，即网格不动，直至流体计算稳定后才释放结构使网格发生运动。因此，结构在被释放后将有一个脉冲运动。

计算表明，当单元尺寸小于0.1m时计算结果与网格尺度无关，A点处的风压为0.03Pa。本算例中稀疏网格的单元尺度基本在5m左右，最大单元尺度接近8m；而稠密网格的单元尺寸在0.5m左右。非稳态计算稳定时刻CFD模型的稀疏网格和稠密网格在A点处的风压发现，前者在A点处的风压为0.04Pa，后者为0.032Pa，相对误差分别为33％和7％。由此可见，CFD网格密度对风压计算结果影响较大。

对于有限元模型而言，单元数量对索膜结构有限元计算速度有着非常重要的影响，随着有限元单元数量增多，非线性计算的收敛性变差，需要更多的迭代步和更小的步长，从而导致了计算速度的下降。另一方面单元数量对于计算精度提高则有限。下表为具有不同单元数量的有限元模型模态计算结果和耗费的机时。

表1不同单元数量的有限元模型计算结果比较

因此，为了既保证数值计算的准确性，同时也为了提高计算效率，实际计算中可以采取有限元模型采用稀疏网格，CFD计算采用稠密网格的策略。这也说明了采用“非一致网格”的必要性。本文采用稀疏一致网格、稠密一致网格和非一致网格对世博轴索膜结构算例进行时长为2min的动力时程计算。

列出四种工况：工况一代表FEM单元数量和CFD单元数量均为3924的一致网格；工况二代表FEM单元数量为3924和CFD单元数量为11935的非一致网格；工况三代表FEM单元数量和CFD单元数量均为11935的一致网格；工况四代表FEM单元数量为3924和CFD单元数量为11935的非一致网格。工况二和工况四的区别在于，工况二采用在计算初始化时就已经设定的FEM单元与CFD单元之间的对应关系，在之后的每一次网格更新计算中无需重新建立这种关系；而工况四却没有采取这一初始化操作，它需要在每一次网格更新计算中重新建立一次非一致网格间的对应关系。

从工况一与工况三的比较发现，单元数量对索膜结构有限元计算速度有着非常重要的影响。这是由于随着有限元单元数量增多，非线性计算的收敛性变差，需要更多的迭代步和更小的步长，从而导致了计算速度的下降。从工况二与工况四的比较发现，非一致网格间对应关系的建立是一个非常耗费机时的运算。以本文为例，要在3924个FEM单元和11935个CFD单元之间建立一次对应关系就需要进行一次3924×11935次运算。因此，在初始化中即建立好非一致网格间的对应关系对于提高流固耦合的计算效率非常重要。

本案例比较了考虑与不考虑流固耦合作用，对世博轴索膜结构风振响应的影响。振动受到结构非线性阻尼的影响。在振动的前60s，由于振幅较大，非线性阻尼也较大，计算阻尼比约为0.15％。随着运动幅值的减小，60s之后的阻尼比也变得非常小。振动频率约为0.6Hz，与膜结构的基频接近。考虑流固耦合的膜结构在自由衰减过程中的总阻尼比为4％，远大于不考虑流固耦合时的阻尼比。振动频率约0.14Hz。由此可见，流固耦合作用使膜结构的动力特性发生了较大的变化。

采用来流为10m/s的紊流计算膜面的风振响应。索膜结构的风致响应可以分解为三个部分：静态成分、准静态成分和瞬态成分。其中静态耦合对应于风荷载的平均响应，准静态耦合对应风荷载的背景响应，瞬态耦合对应风荷载的共振响应。表2列出了考虑与不考虑流固耦合的膜结构风振响应中各成分的分量。

表2振动分量的分解

工况	静态	准静态	瞬态
				不考虑耦合	0.4	0.025	0.05
考虑耦合	0.45	0.025	0

从表中数据可以看出：(1)静态耦合反映了几何形状的变化对结构分析结果的影响，本算例中结构的几何变形使得A点位移趋于增大。(2)准静态耦合反映了脉动风中大尺度涡对结构响应的影响。从本例中可以发现，流固耦合基本不影响大尺度涡对结构的作用。(3)瞬态耦合反映了脉动风的高频部分(小尺度涡)对结构响应的影响。从计算结果来看，风荷载的高频部分对膜结构的风振响应影响很弱以至于可以忽略不计。

实际工程中，设计人员往往喜欢用风振系数的概念来描述结构的风致动力效应。本文如果采用位移风振响应的定义来比较A点的风振系数，则可得到：

${\mathrm{\β}}_0=\frac{D_{\max }}{D_{\mathrm{stead}}}=\frac{0.49}{0.4}=1.225$

${\mathrm{\β}}_c=\frac{D_{\max }}{D_{\mathrm{stead}}}=\frac{0.5}{0.4}=1.25$

式中，β₀和β_c分别代表A点不考虑和考虑流固耦合作用的风振系数；D_max表示位移最大值；D_stead表示平均风荷载下的位移。由此可知，对于世博轴索膜结构而言，如果不考虑流固耦合效应，而采用普通的动力时程计算，结果将偏于不安全。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。

Claims (6)

1.一种基于松耦合技术的索膜结构抗风设计方法，其特征在于，包括下列步骤：

步骤一：依据实际材料的物理性质和结构的初始预应力，通过找形方法建立索膜结构有限元模型；

步骤二：根据步骤一中的找形结果和建筑周围环境，建立索膜结构流体计算模型；

步骤三：利用风洞试验或实测数据或理论模拟的风速时程数据，作为流体计算的入口边界条件；

步骤四：利用流固耦合协同平台进行膜结构抗风设计。

2.根据权利要求1所述的基于松耦合技术的索膜结构抗风设计方法，其特征在于，所述步骤一中找形方法为采用分步计数、逐步刚度逼近的方法建立索膜结构有限元模型。

3.根据权利要求1所述的基于松耦合技术的索膜结构抗风设计方法，其特征在于，所述步骤三中数值风洞模型的建立仅对有限元模型所对应的膜面设置动网格属性，其余部分保持不动。

4.根据权利要求1所述的基于松耦合技术的索膜结构抗风设计方法，其特征在于，所述步骤三中流体计算的非稳态计算采用LES模拟、隐式计算。

5.根据权利要求1所述的基于松耦合技术的索膜结构抗风设计方法，其特征在于，所述有限元模型采用稀疏网格，所述流体计算模型采用稠密网格。

6.根据权利要求1所述的基于松耦合技术的索膜结构抗风设计方法，其特征在于，所述步骤四中利用流固耦合协同平台进行膜结构抗风设计为考虑流固耦合作用进行紊流计算膜面的风振响应。

Images