CN113779831A - 一种基于区域分解的缩聚feti工程数值方法 - Google Patents

一种基于区域分解的缩聚feti工程数值方法 Download PDF

Info

Publication number
CN113779831A
CN113779831A CN202111025339.3A CN202111025339A CN113779831A CN 113779831 A CN113779831 A CN 113779831A CN 202111025339 A CN202111025339 A CN 202111025339A CN 113779831 A CN113779831 A CN 113779831A
Authority
CN
China
Prior art keywords
matrix
boundary
subdomain
rigid
freedom
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111025339.3A
Other languages
English (en)
Other versions
CN113779831B (zh
Inventor
戎宇飞
齐雪超
袁伟
张晓虎
孙秦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202111025339.3A priority Critical patent/CN113779831B/zh
Publication of CN113779831A publication Critical patent/CN113779831A/zh
Application granted granted Critical
Publication of CN113779831B publication Critical patent/CN113779831B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开一种基于区域分解的缩聚FETI工程数值方法，通过构造只与漂浮子域边界自由度相关的刚阵及其刚体模态阵，求出广义柔度矩阵及边界柔度阵。形成的新的“帽子”矩阵不仅能充分利用并保持漂浮子域刚阵的稀疏性，且其阶数等于漂浮子域边界自由度个数，远远小于原始“帽子”矩阵的规模，从而大大减少存储空间，提高广义刚度矩阵三角分解求解效率，使得局部柔度阵的计算变得经济可行。同时，新的“帽子”矩阵计算的中间结果恰好为框架方程组的PCPG迭代求解提供一种具有更优收敛特性的Dirichlet预处理器。

Description

一种基于区域分解的缩聚FETI工程数值方法
技术领域
本发明适用于飞机结构数值仿真领域，具体涉及一种基于区域分解的缩聚FETI工程数值方法。
背景技术
线性方程组的求解是科学与工程计算中数值模拟的关键技术之一。近年来，国家在航空航天领域发展迅速，这些高精尖武器装备的研制，对工程结构的设计提出了精细化、轻量化、可靠性等很高的技术要求。在结构有限元分析时，网格质量和模型精度均要进一步提高，使得线性方程组的规模急剧上升，这对求解的精确性和实时性提出了更大的挑战。伴随着计算机技术的革新及计算机集群系统的发展，高性能计算备受瞩目，逐步成为科技创新的重要手段及推动工程解决方案进步的重要因素，被广泛应用于科研领域之中。而区域分解作为一种适应现代计算机硬件体系的并行算法，其“分而治之”的思想符合高性能计算体系特点，逐渐成为求解大规模科学技术和工程问题的一种有效算法，在结构力学领域中受到广泛的关注。
区域分解是指将原始计算区域分解形成若干较小的子域，从而将原始域中偏微分方程的求解分割成若干子域偏微分方程的耦合求解，与传统的计算方法相比，区域分解法的优点主要表现在：1)将原来的大问题化为若干小问题，缩小计算规模；2)把不规则区域上的问题化成规则区域上的问题求解；3)可以在各子区域使用局部最优网格，而不需要用全局一致网格；4)允许在各子区域使用不同的数学模型，以便更真准确地描述物理现象；5)计算的主要步骤是在各子区域内独立进行的，算法可以高度并行。区域分解思想发展出很多不同的并行求解算法，其中基于拉格朗日乘子混合公式的有限元撕裂合并方法(简称FETI，Finite element tearing and interconnecting)是近年来快速发展的一种算法。文献“Alagrange multiplier based divide and conquer finite element algorithm”给出了FETI算法的思路，即先在分割界面上引入拉格朗日乘子连续性条件以便各子域内部独立解在分割界面上相匹配，再将所有的子问题均衡地映射到多个处理器上，最后，通过多个处理器间的协同运行来获得整个区域的解。在有限元撕裂合并方法中框架方程组的求解是子结构并行分析的关键，其中计算漂浮子域边界柔度阵时，需要利用漂浮子域刚体模态阵构造“帽子”矩阵来消除其刚度矩阵的奇异性，然而“帽子”矩阵的阶数与漂浮子域的自由度数一致且为近似满阵，存储并计算“帽子”矩阵以及对消除了矩阵奇异性的广义刚度矩阵进行三角分解等都会消耗巨大的存储空间和计算时间，是不经济甚至不可能的，这严重限制了FETI算法的发展。
发明内容
为克服FETI算法中，当漂浮结构规模较大时，“帽子”矩阵存储量大且难以计算的问题，本发明提出了一种基于区域分解的缩聚FETI工程数值方法，从减少计算量及存储量的角度出发，以满足工程设计人员现代高性能电脑为目的。
一种基于区域分解的缩聚FETI工程数值方法，包括如下步骤：
步骤1)：原始结构划分为若干子域结构以及分区框架；
步骤2)：根据步骤1)中若干子域结构生成各子域有限元方程；
步骤3)：对步骤2)中各子域有限元方程进行缩聚，得到缩聚边界刚阵、等价载荷；
步骤4)：根据步骤3)得到的缩聚边界刚阵、等价载荷生成各子域局部刚体模态阵；
步骤5)：利用缩聚边界刚阵、局部刚体模态阵计算各子域局部柔度阵；
步骤6)：根据步骤2)至5)获得各个矩阵构建各子域缩聚FETI算法的框架方程组；
步骤7)：利用Dirichlet预处理共轭投影梯度法求解框架方程组，得到各子域界面拉格朗日力和边界刚体模态幅值；
步骤8)：各子域界面拉格朗日力和边界刚体模态幅值回代，得到各子域位移解。
优选地，所述步骤5)，具体包括如下步骤：利用漂浮结构的局部刚体模态阵来消除其缩聚边界刚度矩阵的奇异性，从而得到漂浮结构的一种可逆的广义缩聚边界刚度矩阵：
Figure BDA0003243160580000021
式中，
Figure BDA0003243160580000022
为漂浮子结构可逆的广义缩聚边界刚度矩阵，
Figure BDA0003243160580000023
为子域局部刚体模态，
Figure BDA0003243160580000031
为新的“帽子”矩阵，
利用广义缩聚边界刚度矩阵的逆即可得到漂浮结构边界在外载荷作用下的总位移解
Figure BDA0003243160580000032
再利用其局部刚体模态阵构造投影矩阵：
Figure BDA0003243160580000033
式中，P(s)为投影矩阵，N为子域边界自由度个数，IN为N阶单位矩阵，
将漂浮结构的总位移解
Figure BDA0003243160580000034
左乘投影矩阵P(s)即可分离出漂浮结构边界的弹性变形位移：
Figure BDA0003243160580000035
式中，
Figure BDA0003243160580000036
漂浮结构边界的弹性变形位移，
由式(16)可知，在外载荷作用下漂浮结构产生弹性变形的过程与一般固定结构不同，其呈现出的是一种广义柔度矩阵，子域边界柔度阵与缩聚后边界刚阵的如下关系式：
Figure BDA0003243160580000037
式中，
Figure BDA0003243160580000038
为子域边界柔度阵，
该刚度矩阵缩聚处理技术同样可用于固定子域，由于固定子域的刚度矩阵一般为可逆矩阵，则可将(17)式写为：
Figure BDA0003243160580000039
式中，N为子域边界自由度个数，IN为N阶单位矩阵，
提出从解方程的角度来求解式(18)，其中
Figure BDA00032431605800000310
为系数矩阵，单位阵IN的每一列可看作一个右端项，局部柔度矩阵
Figure BDA00032431605800000311
的每一列可看作对应右端项的解向量，这样局部柔度阵的计算就可看成多右端项线性方程组的求解问题，如(19)式所示，
Figure BDA00032431605800000312
其中，
Figure BDA00032431605800000313
为局部柔度矩阵
Figure BDA00032431605800000314
的每一列向量，Ii为单位阵IN的每一列向量，
将系数矩阵
Figure BDA00032431605800000315
进行直接分解，从单位阵IN中提取出于该子域边界某一自由度对应的一列，可以看作在相应边界自由度上施加单位载荷；再利用三角回代求出解向量，此解向量为在单位载荷Ii作用下结构所有自由度的响应；进一步利用提取算子从解向量中提取出仅与该子域边界相关的元素，即构成
Figure BDA00032431605800000316
的某一列，重复上述过程，依次计算出该子域所有边界自由度下解向量，即完成
Figure BDA0003243160580000041
的求解。
优选地，所述步骤6)，具体包括如下步骤：子域分割边界自由度构造的局部刚体模态阵与由等价载荷和框架作用到子域分割边界上的拉格朗日力构成的合外力需要满足正交关系，即子域自平衡条件：
Figure BDA0003243160580000042
式中，
Figure BDA0003243160580000043
为子域局部刚体模态阵的转置，
Figure BDA0003243160580000044
为等价载荷，
Figure BDA0003243160580000045
表示框架作用于子域s边界上的拉格朗日力，
分区框架是一个虚拟的界面连接结构，起到了将分裂的各个子域缝合到一起的作用，即通过分区框架上的耦合节点使相邻子域的界面位移协调一致，对任一子域s，界面位移协调方程表示如下，
Figure BDA0003243160580000046
式中，
Figure BDA0003243160580000047
为界面位移协调算子，是一个Boolean矩阵，其元素值只有0和1，相当于一个装配阵，ugb为分区框架上所有节点位移，
分区框架的另一个作用是实现相邻子域的界面节点力之间的耦合，分区框架自平衡方程可写为：
Figure BDA0003243160580000048
式中，
Figure BDA0003243160580000049
利用缩聚有限元方程的位移通解式、子域自平衡条件式(20)以及框架位移协调条件式(21)和框架载荷平衡条件式(22)，可以形成缩聚FETI算法的框架方程组：
Figure BDA00032431605800000410
式中，
Figure BDA0003243160580000051
优选地，所述缩聚有限元方程的位移通解式为：
Figure BDA0003243160580000052
式中，
Figure BDA0003243160580000053
为与子域边界自由度相关的局部柔度阵，
Figure BDA0003243160580000054
为等价载荷，
Figure BDA0003243160580000055
表示框架作用于子域s边界上的拉格朗日力，
Figure BDA0003243160580000056
为与子域边界自由度相关的局部刚体模态阵，
Figure BDA0003243160580000057
为子域边界刚体模态幅值。
优选地，所述步骤7)中，Dirichlet预处理矩阵为：
Figure BDA0003243160580000058
式中，
Figure BDA0003243160580000059
为子域缩聚后边界刚阵，
Figure BDA00032431605800000510
为子域边界-边界自由度对应的刚阵，
Figure BDA00032431605800000511
为子域内部-内部自由度对应的刚阵，
Figure BDA00032431605800000512
为子域内部-边界自由度对应的刚阵，
Figure BDA00032431605800000513
为子域边界-内部自由度对应的刚阵，
利用Dirichlet预处理共轭投影梯度法对式(24)进行求解后得到各子域界面拉格朗日力和漂浮子域边界刚体模态幅值。
基于上述技术方案，本发明的有益效果是：本发明利用单位阵高效计算局部柔度阵的直接求解技术，并结合刚度矩阵缩聚求解技术及局部刚体模态阵的概念，提出一种计算漂浮结构局部柔度阵的新思路，即通过构造只与漂浮子域边界自由度相关的刚阵及其刚体模态阵，求出广义柔度矩阵及边界柔度阵。这样形成的新的“帽子”矩阵不仅能充分利用并保持漂浮子域刚阵的稀疏性，且其阶数等于漂浮子域边界自由度个数，远远小于原始“帽子”矩阵的规模，从而大大减少存储空间，提高广义刚度矩阵三角分解求解效率，使得局部柔度阵的计算变得经济可行。与此同时，新的“帽子”矩阵计算的中间结果恰好为框架方程组的PCPG迭代求解提供一种具有更优收敛特性的Dirichlet预处理器。
附图说明
图1是本发明的缩聚FETI算法流程示意图；
图2是平板有限元模型示意图；
图3是平板有限元模型2子域分区示意图；
图4是原FETI算法中板弯模型各子域位移变化示意图；
图5是缩聚FETI策略中板弯模型各子域位移变化示意图；
图6是机翼模型内部梁、肋结构的分布示意图；
图7是机翼平面形状主要几何尺寸；
图8是机翼模型几何分区及各子域有限元模型示意图；
图9是缩聚FETI-Dirichlet预处理并行策略下PCPG迭代收敛曲线。
具体实施方式
本发明从减少计算量及存储量的角度出发，以满足工程设计人员现代高性能电脑为目的，提出了利用单位阵高效计算局部柔度阵的直接求解技术，并结合刚度矩阵缩聚求解技术及局部刚体模态阵的概念，提出一种计算漂浮结构局部柔度阵的新思路，即通过构造只与漂浮子域边界自由度相关的刚阵及其刚体模态阵，求出广义柔度矩阵及边界柔度阵。这样形成的新的“帽子”矩阵不仅能充分利用并保持漂浮子域刚阵的稀疏性，且其阶数等于漂浮子域边界自由度个数，远远小于原始“帽子”矩阵的规模，从而大大减少存储空间，提高广义刚度矩阵三角分解求解效率，使得局部柔度阵的计算变得经济可行。与此同时，新的“帽子”矩阵计算的中间结果恰好为框架方程组的PCPG迭代求解提供一种具有更优收敛特性的Dirichlet预处理器。如图1所示，主要步骤包括：
步骤1：有限元结构模型子域划分。有限元撕裂合并算法将结构模型划分为由“分区框架”上的耦合节点间接地联系在一起的若干互不直接相连的Ns个独立子区域。“分区框架”是一种虚拟的、由各子域界面节点共同组成的连接器，它受到由相邻子域传递过来的界面节点力的作用，且每个耦合节点的位移都与相邻界面节点的位移一致。同时，各个子域的节点编号、单元编号等完全独立，因而各子域的网格划分可独立并行进行。
步骤2：子域有限元方程。完成步骤1子域划分后，每一个子域可视为一个独立的子结构，每一个子域的位移边界和受力载荷均直接继承自原结构的部分边界约束。对于每一子域，有下面的有限元方程：
K(s)u(s)＝p(s)(s＝1,2,…Ns) (1)
式中，Ns为子域总数，K(s)表示子域s的刚度矩阵，u(s)表示子域的节点位移向量，p(s)表示子域的节点总载荷。
子域s的总体载荷向量p(s)可表示为原结构中该子域节点载荷和框架作用于子域界面节点力之和：
Figure BDA0003243160580000071
式中，
Figure BDA0003243160580000072
为原始结构中该子域节点载荷，
Figure BDA0003243160580000073
表示框架作用于子域s边界上的拉格朗日力，B(s)为Boolean矩阵，反映了子域的界面节点位移向量在子域的总体节点位移向量中的装配关系，其元素值仅包含0和1。利用B(s)T可以实现从子域全体节点向量中提取出子域的界面节点向量：
Figure BDA0003243160580000074
式中，
Figure BDA0003243160580000075
为子域的界面节点位移。
步骤3：矩阵缩聚求解技术。将式(1)中子域位移u(s)按内部自由度i和边界自由度b组织成如下形式：
Figure BDA0003243160580000076
式中，u(s)为子域节点位移，
Figure BDA0003243160580000077
为子域内部节点位移，
Figure BDA0003243160580000078
为边界节点位移。则子域有限元方程可以写成：
Figure BDA0003243160580000079
式中，
Figure BDA00032431605800000710
为子域内部-内部自由度对应的刚阵，
Figure BDA00032431605800000711
为子域内部-边界自由度对应的刚阵，
Figure BDA0003243160580000081
为子域边界-内部自由度对应的刚阵，
Figure BDA0003243160580000082
为子域边界-边界自由度对应的刚阵，fi (s)为子域内部自由度节点力，
Figure BDA0003243160580000083
为子域边界自由度节点力。
利用式(5)的第一行解得子域内部节点位移如下：
Figure BDA0003243160580000084
将式(6)代入式(5)第二式并消除子域内部节点位移可得缩聚平衡方程：
Figure BDA0003243160580000085
其中，
Figure BDA0003243160580000086
为子域缩聚后边界刚阵，
Figure BDA0003243160580000087
为等价载荷，具体地：
Figure BDA0003243160580000088
Figure BDA0003243160580000089
进行区域分解后，将出现两种形式的子域结构：固定子域和漂浮子域。式(7)的子域边界节点位移通解可以写成：
Figure BDA00032431605800000810
式中，
Figure BDA00032431605800000811
为与子域边界自由度相关的局部柔度阵，
Figure BDA00032431605800000812
为与子域边界自由度相关的局部刚体模态阵，
Figure BDA00032431605800000813
为子域边界刚体模态幅值。对于固定子域可视其刚体模态阵和刚体位移都为0，最终求解时再将0项舍去。
步骤4：子域局部刚体模态阵。考虑三维空间中完全自由的子域，其有且仅有6个刚体运动自由度，节点i相对于求矩参考点的运动可以写为：
Figure BDA00032431605800000814
式中，I3为3×3的单位阵，(xi,yi,zi)和(x0,y0,z0)分别为节点i和求矩参考点的坐标。
因此，子域刚体模态可以写为：
Figure BDA00032431605800000815
式中，R(s)为子域刚体模态阵，N为子域节点数。
进一步可得子域局部刚体模态阵：
Figure BDA00032431605800000816
式中，
Figure BDA0003243160580000091
为子域局部刚体模态。
步骤5：子域局部柔度阵。对于约束完备的固定结构来说，柔度矩阵可由刚度矩阵直接求逆得到。然而，对于漂浮结构来说，因其含有刚体位移，其刚度矩阵不可逆，其刚度矩阵和柔度矩阵不是简单的互逆关系，无法通过其刚度矩阵求逆来得到其弹性变形位移解。需要利用漂浮结构的局部刚体模态阵来消除其缩聚边界刚度矩阵的奇异性，从而得到漂浮结构的一种可逆的广义缩聚边界刚度矩阵：
Figure BDA0003243160580000092
式中，
Figure BDA0003243160580000093
为漂浮子结构可逆的广义缩聚边界刚度矩阵，
Figure BDA0003243160580000094
为子域局部刚体模态，
Figure BDA0003243160580000095
为本发明新的“帽子”矩阵。
分析发现，漂浮结构受外力作用产生的位移包含结构弹性变形位移和结构刚体位移。弹性变形位移代表了结构的承载能力，是结构静力分析主要的关注对象，需要把弹性变形位移从漂浮结构总位移中分离出来。首先，利用广义缩聚边界刚度矩阵的逆即可得到漂浮结构边界在外载荷作用下的总位移解
Figure BDA0003243160580000096
再利用其局部刚体模态阵构造投影矩阵：
Figure BDA0003243160580000097
式中，P(s)为投影矩阵，N为子域边界自由度个数，IN为N阶单位矩阵。
将漂浮结构的总位移解
Figure BDA0003243160580000098
左乘投影矩阵P(s)即可分离出漂浮结构边界的弹性变形位移：
Figure BDA0003243160580000099
式中，
Figure BDA00032431605800000910
漂浮结构边界的弹性变形位移。
由式(16)可知，在外载荷作用下漂浮结构产生弹性变形的过程与一般固定结构不同，其呈现出的是一种广义柔度矩阵，子域边界柔度阵与缩聚后边界刚阵的如下关系式：
Figure BDA00032431605800000911
式中，
Figure BDA00032431605800000912
为子域边界柔度阵。
此外，该刚度矩阵缩聚处理技术同样可用于固定子域，由于固定子域的刚度矩阵一般为可逆矩阵，则可将(17)式写为：
Figure BDA0003243160580000101
式中，N为子域边界自由度个数，IN为N阶单位矩阵。
本发明提出从解方程的角度来求解式(18)，其中
Figure BDA0003243160580000102
为系数矩阵，单位阵IN的每一列可看作一个右端项，局部柔度矩阵
Figure BDA0003243160580000103
的每一列可看作对应右端项的解向量，这样局部柔度阵的计算就可看成多右端项线性方程组的求解问题，如(19)式所示，
Figure BDA0003243160580000104
其中，
Figure BDA0003243160580000105
为局部柔度矩阵
Figure BDA0003243160580000106
的每一列向量，Ii为单位阵IN的每一列向量。
将系数矩阵
Figure BDA0003243160580000107
进行直接分解，从单位阵IN中提取出于该子域边界某一自由度对应的一列，可以看作在相应边界自由度上施加单位载荷；再利用三角回代求出解向量，此解向量为在单位载荷Ii作用下结构所有自由度的响应；进一步利用提取算子从解向量中提取出仅与该子域边界相关的元素，即构成
Figure BDA0003243160580000108
的某一列。重复上述过程，依次计算出该子域所有边界自由度下解向量，即完成
Figure BDA0003243160580000109
的求解。
对于漂浮子域结构来说，求解式(17)过程和固定子域相同，需要注意的是式(19)式右端不是单位矩阵。
步骤6：缩聚FETI算法中的框架方程组。子域分割边界自由度构造的局部刚体模态阵与由等价载荷和框架作用到子域分割边界上的拉格朗日力构成的合外力需要满足正交关系，即子域自平衡条件：
Figure BDA00032431605800001010
式中，
Figure BDA00032431605800001011
为子域局部刚体模态阵的转置，
Figure BDA00032431605800001012
为等价载荷，
Figure BDA00032431605800001013
表示框架作用于子域s边界上的拉格朗日力，
分区框架是一个虚拟的界面连接结构，起到了将分裂的各个子域缝合到一起的作用，即通过分区框架上的耦合节点使相邻子域的界面位移协调一致。对任一子域s，界面位移协调方程表示如下，
Figure BDA00032431605800001014
式中，
Figure BDA0003243160580000111
为界面位移协调算子，是一个Boolean矩阵，其元素值只有0和1，相当于一个装配阵，ugb为分区框架上所有节点位移。
分区框架的另一个作用是实现相邻子域的界面节点力之间的耦合。分区框架自平衡方程可写为：
Figure BDA0003243160580000112
式中，
Figure BDA0003243160580000113
利用缩聚有限元方程的位移通解式(10)、子域自平衡条件式(20)以及框架位移协调条件式(21)和框架载荷平衡条件式(22)，可以形成缩聚FETI算法的框架方程组：
Figure BDA0003243160580000114
式中，
Figure BDA0003243160580000115
步骤7：Dirichlet预处理共轭投影梯度法。对框架方程组式(24)的迭代求解，一般采用共轭投影梯度法，为了加快算法收敛速度，使用Dirichlet矩阵进行预处理：
Figure BDA0003243160580000116
式中，
Figure BDA0003243160580000117
为Dirichlet预处理矩阵。这恰好与式(8)一样，即在缩聚FETI算法中，已经形成了Dirichlet预处理矩阵，这样就能减少重复的计算，提高求解效率。
利用Dirichlet预处理的共轭投影梯度法对式(24)进行求解后得到各子域界面拉格朗日力和漂浮子域边界刚体模态幅值，然后利用式(10)得到各个子域边界位移
Figure BDA0003243160580000122
再根据式(6)便可得到各子域内部自由度位移解。
实施例一
本实施例为两子域板弯模型，平板模型几何长度L＝100mm，高h＝10mm，厚t＝1mm，材料弹性模量E＝70Gpa，泊松比μ＝0.33，其有限元网格采用四边形shell单元，平板一端采用固支约束(节点1和6)。为观察平板位移解产生的显著变化，平板另一端施加1000N载荷(节点5)，如图2所示。
步骤1：将板弯有限元模型简单地划分成2个子域，如图3所示，分别采用FETI算法和缩聚FETI算法两种并行策略组装框架方程组并进行求解。
步骤2：为避免迭代法求解时因控制收敛门槛而带来的数值解的逼近误差，对于该简单模型形成的低阶框架方程组可采用直接法进行求解。
步骤3：两种并行求解策略下各子域弹性位移、刚体位移和总位移产生的子域结构变形示意图分别如图4和图5所示，相应数值结果分别如表1～表3所示，其中，固定子域(子域)刚体位移为0。
由本实施例步骤完成的结果可以得出以下基本结论：
1.在表1中，两种并行策略下，固定子域1的弹性位移一致，而漂浮子域2的弹性位移不同；在表2中，两种并行策略下，固定子域1的刚体位移一致，均为0，而漂浮子域2的刚体位移不同；但在表3中，两种并行策略下，整体有限元模型的总位移保持一致，同时，更高精度的数值计算显示两种并行策略下的总位移之间的误差为1.284901E-11，即表明了本发明提出缩聚FETI方法的有效性。
表1两种并行策略下各子域弹性位移(mm)
Figure BDA0003243160580000121
Figure BDA0003243160580000131
表2两种并行策略下漂浮子域(子域1)刚体位移(mm)
Figure BDA0003243160580000132
表3两种并行策略下模型节点总位移(mm)
Figure BDA0003243160580000133
2.两种并行策略下所计算出的漂浮子域2的“弹性位移”虽然不同，但利用投影分离思想可以发现，两种“弹性位移”之间存在如下关系式，
ue＝[I-R(RTR)-1RT]usje (27)
其中，ue和usje分别为FETI算法和缩聚FETI算法两种并行策略下漂浮子域2的弹性位移解；R为与漂浮子域2所有节点自由度相关的刚体模态阵；I为单位阵，其维数为漂浮子域2节点自由度数；与式(15)表达类似，I-R(RTR)-1RT为相对于子域所有节点自由度的投影矩阵。这表明在缩聚FETI算法中，计算出的漂浮子域的“弹性位移”隐藏着一定的刚体运动。
实施例二
本实施例为一典型翼盒结构的机翼模型，通过网格细化来形成大规模结构有限元待求解问题，并采用4节点集群平台进行并行数值求解。机翼模型内部梁、肋结构的分布及机翼平面形状的主要几何尺寸分别如图6和图7所示。
步骤1：先将机翼几何模型划分为4个较小子域，再利用软件Hypermesh对区域分解后的模型进行网格剖分，整个机翼有限元模型中共包含52593个四边形shell单元和3277个两节点Beam单元，其中，上下蒙皮材料常数设为杨氏模量E＝70GPa，泊松比μ＝0.33，厚度t＝2mm，翼梁和翼肋的材料常数简单设为杨氏模量E＝210GPa，泊松比μ＝0.33，且梁腹板厚度t＝8mm，梁缘条截面积为96mm2，肋腹板厚度t＝4mm，肋缘条截面积为48mm2。
该机翼模型几何分区及各子域有限元模型如图8所示。有限元撕裂合并求解中框架方程组的相关参数信息如表4所示。
表4机翼有限元模型4子域分区参数信息
Figure BDA0003243160580000141
表中，Ne表示有限元模型中单元个数；Np表示有限元模型中节点个数；N表示整个有限元模型待求解自由度个数；Npf表示分区框架上节点个数；Ndf表示分区框架上待求解自由度个数；Nf表示框架方程组阶数；Nsbt为框架方程组经投影操作后需要进行迭代求解的方程组阶数；Nst表示各子域自由度个数；Ns表示各子域待求解自由度个数，即子域有限元方程组阶数；Nb表示各子域边界待求解自由度个数。
步骤2：为了实现对该机翼模型的静力求解，还需设置边界条件。在机翼有限元模型的根部施加固支约束，并在梁结构下缘条上施加均布载荷，即在下缘条有限元网格节点上均施加1N的节点集中力，总载荷为965N。
步骤3：对于该机翼结构4子域分区有限元模型，采用缩聚FETI并行策略组装框架方程组，在PCPG迭代求解过程中采用Dirichlet预处理器，并取搜索方向重正交化处理时正交基的个数分别为100、300和500。取投影残向量h的二范数||h||2和解向量λb增量绝对值最大值|Δλb|max分别作为判断迭代过程是否收敛的两种评定指标，并将收敛门槛值取为右端项hλ二范数的1.0E-9倍，则缩聚FETI-Dirichlet预处理并行策略下PCPG迭代收敛曲线分别如图9所示。
步骤4：利用DSS直接稀疏求解器可解出该机翼整体有限元模型的结构位移响应，并以此为参照，可获得区域分解并行求解时的结构位移响应与DSS数值解之间的相对误差。表5给出了不同正交基个数下PCPG迭代收敛结果以及相应位移响应与DSS数值解之间的相对误差。
表5翼盒结构4子域分区时缩聚FETI并行策略下PCPG迭代求解数值结果
Figure BDA0003243160580000151
由本实施例步骤完成的结果可以得出以下基本结论：
1.由表5可知，缩聚FETI并行求解时的结构位移响应与DSS数值解之间的相对误差在1.0E-4量级，满足工程使用需求。
2.从存储空间需求上来说，采用原始FETI并行策略时，在计算漂浮子域边界柔度矩阵Fb的过程中，需要利用帽子矩阵R(RTR)-1RT来消除其刚度矩阵的奇异性，即需构造K+R(RTR)-1RT，以表4子域2为例，该子域为漂浮子域且包含76302个自由，这就需要构造76302阶近似满阵形式的帽子矩阵，约需要43.38GB存储空间，极大地消耗存储空间和计算时间，这是不经济且不现实的；若采用本发明提出的缩聚FETI并行策略，仅需要构造
Figure BDA0003243160580000152
由于该子域边界自由度个数为1098，构造满阵形式的帽子矩阵时仅需约9.20MB存储空间，可极大节省存储空间，提升计算效率和待求解问题规模。

Claims (5)

1.一种基于区域分解的缩聚FETI工程数值方法，其特征在于，包括如下步骤：
步骤1)：原始结构划分为若干子域结构以及分区框架；
步骤2)：根据步骤1)中若干子域结构生成各子域有限元方程；
步骤3)：对步骤2)中各子域有限元方程进行缩聚，得到缩聚边界刚阵、等价载荷；
步骤4)：根据步骤3)得到的缩聚边界刚阵、等价载荷生成各子域局部刚体模态阵；
步骤5)：利用缩聚边界刚阵、局部刚体模态阵计算各子域局部柔度阵；
步骤6)：根据步骤2)至5)获得各个矩阵构建各子域缩聚FETI算法的框架方程组；
步骤7)：利用Dirichlet预处理共轭投影梯度法求解框架方程组，得到各子域界面拉格朗日力和边界刚体模态幅值；
步骤8)：各子域界面拉格朗日力和边界刚体模态幅值回代，得到各子域位移解。
2.根据权利要求1所述的一种基于区域分解的缩聚FETI工程数值方法，其特征在于，所述步骤5)，具体包括如下步骤：利用漂浮结构的局部刚体模态阵来消除其缩聚边界刚度矩阵的奇异性，从而得到漂浮结构的一种可逆的广义缩聚边界刚度矩阵：
Figure FDA0003243160570000011
式中，
Figure FDA0003243160570000012
为漂浮子结构可逆的广义缩聚边界刚度矩阵，
Figure FDA0003243160570000013
为子域局部刚体模态，
Figure FDA0003243160570000014
为新的“帽子”矩阵，
利用广义缩聚边界刚度矩阵的逆即可得到漂浮结构边界在外载荷作用下的总位移解
Figure FDA0003243160570000015
再利用其局部刚体模态阵构造投影矩阵：
Figure FDA0003243160570000016
式中，P(s)为投影矩阵，N为子域边界自由度个数，IN为N阶单位矩阵，
将漂浮结构的总位移解
Figure FDA0003243160570000017
左乘投影矩阵P(s)即可分离出漂浮结构边界的弹性变形位移：
Figure FDA0003243160570000018
式中，
Figure FDA0003243160570000019
漂浮结构边界的弹性变形位移，
由式(16)可知，在外载荷作用下漂浮结构产生弹性变形的过程与一般固定结构不同，其呈现出的是一种广义柔度矩阵，子域边界柔度阵与缩聚后边界刚阵的如下关系式：
Figure FDA0003243160570000021
式中，
Figure FDA0003243160570000022
为子域边界柔度阵，
该刚度矩阵缩聚处理技术同样可用于固定子域，由于固定子域的刚度矩阵一般为可逆矩阵，则可将(17)式写为：
Figure FDA0003243160570000023
式中，N为子域边界自由度个数，IN为N阶单位矩阵，
提出从解方程的角度来求解式(18)，其中
Figure FDA0003243160570000024
为系数矩阵，单位阵IN的每一列可看作一个右端项，局部柔度矩阵
Figure FDA0003243160570000025
的每一列可看作对应右端项的解向量，这样局部柔度阵的计算就可看成多右端项线性方程组的求解问题，如(19)式所示，
Figure FDA0003243160570000026
其中，
Figure FDA0003243160570000027
为局部柔度矩阵
Figure FDA0003243160570000028
的每一列向量，Ii为单位阵IN的每一列向量，
将系数矩阵
Figure FDA0003243160570000029
进行直接分解，从单位阵IN中提取出于该子域边界某一自由度对应的一列，可以看作在相应边界自由度上施加单位载荷；再利用三角回代求出解向量，此解向量为在单位载荷Ii作用下结构所有自由度的响应；进一步利用提取算子从解向量中提取出仅与该子域边界相关的元素，即构成
Figure FDA00032431605700000210
的某一列，重复上述过程，依次计算出该子域所有边界自由度下解向量，即完成
Figure FDA00032431605700000211
的求解。
3.根据权利要求1所述的一种基于区域分解的缩聚FETI工程数值方法，其特征在于，所述步骤6)，具体包括如下步骤：子域分割边界自由度构造的局部刚体模态阵与由等价载荷和框架作用到子域分割边界上的拉格朗日力构成的合外力需要满足正交关系，即子域自平衡条件：
Figure FDA00032431605700000212
式中，
Figure FDA00032431605700000213
为子域局部刚体模态阵的转置，
Figure FDA00032431605700000214
为等价载荷，
Figure FDA00032431605700000215
表示框架作用于子域s边界上的拉格朗日力，
分区框架是一个虚拟的界面连接结构，起到了将分裂的各个子域缝合到一起的作用，即通过分区框架上的耦合节点使相邻子域的界面位移协调一致，对任一子域s，界面位移协调方程表示如下，
Figure FDA0003243160570000031
式中，
Figure FDA0003243160570000032
为界面位移协调算子，是一个Boolean矩阵，其元素值只有0和1，相当于一个装配阵，ugb为分区框架上所有节点位移，
分区框架的另一个作用是实现相邻子域的界面节点力之间的耦合，分区框架自平衡方程可写为：
Figure FDA0003243160570000033
式中，
Figure FDA0003243160570000034
利用缩聚有限元方程的位移通解式、子域自平衡条件式(20)以及框架位移协调条件式(21)和框架载荷平衡条件式(22)，可以形成缩聚FETI算法的框架方程组：
Figure FDA0003243160570000035
式中，
Figure FDA0003243160570000036
Figure FDA0003243160570000037
4.根据权利要求3所述的一种基于区域分解的缩聚FETI工程数值方法，其特征在于，所述缩聚有限元方程的位移通解式为：
Figure FDA0003243160570000041
式中，
Figure FDA0003243160570000042
为与子域边界自由度相关的局部柔度阵，
Figure FDA0003243160570000043
为等价载荷，
Figure FDA0003243160570000044
表示框架作用于子域s边界上的拉格朗日力，
Figure FDA0003243160570000045
为与子域边界自由度相关的局部刚体模态阵，
Figure FDA0003243160570000046
为子域边界刚体模态幅值。
5.根据权利要求3所述的一种基于区域分解的缩聚FETI工程数值方法，其特征在于，所述步骤7)中，Dirichlet预处理矩阵为：
Figure FDA0003243160570000047
式中，
Figure FDA0003243160570000048
为子域缩聚后边界刚阵，
Figure FDA0003243160570000049
为子域边界-边界自由度对应的刚阵，
Figure FDA00032431605700000410
为子域内部-内部自由度对应的刚阵，
Figure FDA00032431605700000411
为子域内部-边界自由度对应的刚阵，
Figure FDA00032431605700000412
为子域边界-内部自由度对应的刚阵，
利用Dirichlet预处理共轭投影梯度法对式(24)进行求解后得到各子域界面拉格朗日力和漂浮子域边界刚体模态幅值。
CN202111025339.3A 2021-09-02 2021-09-02 一种基于区域分解的缩聚feti工程数值方法 Active CN113779831B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111025339.3A CN113779831B (zh) 2021-09-02 2021-09-02 一种基于区域分解的缩聚feti工程数值方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111025339.3A CN113779831B (zh) 2021-09-02 2021-09-02 一种基于区域分解的缩聚feti工程数值方法

Publications (2)

Publication Number Publication Date
CN113779831A true CN113779831A (zh) 2021-12-10
CN113779831B CN113779831B (zh) 2023-04-18

Family

ID=78840855

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111025339.3A Active CN113779831B (zh) 2021-09-02 2021-09-02 一种基于区域分解的缩聚feti工程数值方法

Country Status (1)

Country Link
CN (1) CN113779831B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114372361A (zh) * 2021-12-23 2022-04-19 北京理工大学 基于bddc区域分解并行算法的粗网格选取方法
CN116090113A (zh) * 2022-10-08 2023-05-09 盐城工业职业技术学院 一种块体结构可变节点悬挂单元网格有限元模型

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0969396A2 (en) * 1998-07-01 2000-01-05 Tadahiko Kawai Numerical analyzing method, element analyzing blocks, and processing apparatus for numerical analysis
CN103324850A (zh) * 2013-06-21 2013-09-25 上海交通大学 基于多文件流的有限元两级分区两次缩聚并行方法
CN104572109A (zh) * 2015-01-19 2015-04-29 上海交通大学 两级分区两次缩聚并行计算系统开发方法及并行计算系统
CN108038277A (zh) * 2017-11-29 2018-05-15 中国空间技术研究院 一种航天器有限元模型的二次缩聚方法
CN110704950A (zh) * 2019-09-27 2020-01-17 西北工业大学 一种消除自由飞配平载荷下飞机变形中刚性位移的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0969396A2 (en) * 1998-07-01 2000-01-05 Tadahiko Kawai Numerical analyzing method, element analyzing blocks, and processing apparatus for numerical analysis
CN103324850A (zh) * 2013-06-21 2013-09-25 上海交通大学 基于多文件流的有限元两级分区两次缩聚并行方法
CN104572109A (zh) * 2015-01-19 2015-04-29 上海交通大学 两级分区两次缩聚并行计算系统开发方法及并行计算系统
CN108038277A (zh) * 2017-11-29 2018-05-15 中国空间技术研究院 一种航天器有限元模型的二次缩聚方法
CN110704950A (zh) * 2019-09-27 2020-01-17 西北工业大学 一种消除自由飞配平载荷下飞机变形中刚性位移的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A. MARKOPOULOS, L. ŘÍHA, T. BRZOBOHATÝ. ET AL.: "Treatment of Singular Matrices in the Hybrid Total FETI Method", 《DOMAIN DECOMPOSITION METHODS IN SCIENCE AND ENGINEERING XXIII. LECTURE NOTES IN COMPUTATIONAL SCIENCE AND ENGINEERING》 *
张晓虎，孙秦: "结构非匹配网格区域分裂的并行数值计算研究", 《西北工业大学学报》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114372361A (zh) * 2021-12-23 2022-04-19 北京理工大学 基于bddc区域分解并行算法的粗网格选取方法
CN114372361B (zh) * 2021-12-23 2024-05-28 北京理工大学 基于bddc区域分解并行算法的粗网格选取方法
CN116090113A (zh) * 2022-10-08 2023-05-09 盐城工业职业技术学院 一种块体结构可变节点悬挂单元网格有限元模型
CN116090113B (zh) * 2022-10-08 2023-09-08 盐城工业职业技术学院 一种块体结构可变节点悬挂单元网格有限元模型

Also Published As

Publication number Publication date
CN113779831B (zh) 2023-04-18

Similar Documents

Publication Publication Date Title
Dunning et al. Coupled aerostructural topology optimization using a level set method for 3D aircraft wings
CN113779831B (zh) 一种基于区域分解的缩聚feti工程数值方法
Farhat et al. The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems
CN108763658B (zh) 基于等几何方法的组合薄壁结构固有频率设计方法
Martins et al. A coupled-adjoint sensitivity analysis method for high-fidelity aero-structural design
Mian et al. Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach
Sun et al. Topology optimization based on level set for a flexible multibody system modeled via ANCF
Gray et al. Geometrically nonlinear high-fidelity aerostructural optimization for highly flexible wings
Yuan et al. An equivalent modeling method for honeycomb sandwich structure based on orthogonal anisotropic solid element
Zhang et al. Efficient monolithic solution algorithm for high-fidelity aerostructural analysis and optimization
CN114756934B (zh) 一种三维多尺度超材料结构优化设计方法
Hui et al. A data-driven CUF-based beam model based on the tree-search algorithm
Dunning et al. Aerostructural level set topology optimization for a common research model wing
Sidorov et al. Nonlocal in time model of material damping in composite structural elements dynamic analysis
CN113505405B (zh) 等效荷载获取方法、基于等效荷载的拓扑优化方法及系统
CN117540590B (zh) 一种约束阻尼板壳结构的建模方法、装置及计算机设备
CN116305589B (zh) 一种直升机桨叶结构降阶分析方法、系统、设备及介质
CN109948253B (zh) 薄板无网格Galerkin结构模态分析的GPU加速方法
Ishihara et al. Computational fluid–structure interaction framework for passive feathering and cambering in flapping insect wings
Mas Colomer et al. Similarity maximization of a scaled aeroelastic flight demonstrator via multidisciplinary optimization
ISHIDA et al. Topology optimization for maximizing linear buckling load based on level set method
CN115292953A (zh) 一种用于分析二维周期性非均质结构的力学仿真分析方法
CN114091124A (zh) 一种基于拓扑优化设计的负热膨胀超材料夹芯板的制备方法
CN110704954B (zh) 变体飞行器柔性蒙皮胞状支撑体的拓扑优化方法及系统
Biswas et al. Global load balancing with parallel mesh adaption on distributed-memory systems

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant