CN113705040B - 结构损伤分析的近场有限元法及在商用软件中的实现方法 - Google Patents

结构损伤分析的近场有限元法及在商用软件中的实现方法 Download PDF

Info

Publication number
CN113705040B
CN113705040B CN202110885638.8A CN202110885638A CN113705040B CN 113705040 B CN113705040 B CN 113705040B CN 202110885638 A CN202110885638 A CN 202110885638A CN 113705040 B CN113705040 B CN 113705040B
Authority
CN
China
Prior art keywords
omega
finite element
near field
unit
bond
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.)
Active
Application number
CN202110885638.8A
Other languages
English (en)
Other versions
CN113705040A (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.)
Dalian University of Technology
Original Assignee
Dalian University of Technology
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 Dalian University of Technology filed Critical Dalian University of Technology
Priority to CN202110885638.8A priority Critical patent/CN113705040B/zh
Publication of CN113705040A publication Critical patent/CN113705040A/zh
Priority to US18/024,780 priority patent/US20230325560A1/en
Priority to PCT/CN2022/100311 priority patent/WO2023011029A1/zh
Application granted granted Critical
Publication of CN113705040B publication Critical patent/CN113705040B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • 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)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及结构损伤分析的近场有限元法及在商用软件中的实现方法，结构损伤分析的近场有限元法包括步骤：(A1)由目标结构得到有限单元；(A2)由有限单元生成近场单元；(A3)计算总体载荷向量；(A4)计算总体刚度矩阵；(A5)求解总体节点位移；(A6)判断当前计算结果是否收敛，若是则进入(A7)；反之返回(A4)；(A7)判断当前施加载荷是否超过预设值，若是则进入(A8)；反之增加载荷后返回(A3)；(A8)输出结果；实现方法包括操作：(B1)由有限单元生成近场单元并输入到软件中；(B2)将(A4)中单元刚度矩阵的公式编写为软件的单元开发子程序；(B3)设定(A6)中函数μ(ξ,t)的更新方式。本发明的方法计算过程简明，与商用软件兼容性好。

Description

结构损伤分析的近场有限元法及在商用软件中的实现方法
技术领域
本发明属于计算力学领域，具体涉及一种结构损伤分析的近场有限元法及在商用软件中的实现方法。
背景技术
结构及其材料在外部载荷作用下的破坏问题是固体力学领域的百年难题，具有重要的工程应用价值，在航空航天、机械制造、土木水利等重大工程领域都不可避免地遇到此类问题。由于破坏问题涉及到材料的缺陷、损伤和断裂等非连续变形，所以基于连续性假设的经典连续介质力学模型已不再适用。上世纪20年代兴起的线弹性断裂力学和上世纪50年代兴起的连续损伤力学都是固体力学家们对材料破坏过程的有益探索，且收获颇丰。然而，在传统的线弹性断裂力学中，必须人为地引入初始裂纹。这种方式并未考虑裂纹成核的物理机理，因此在预测裂纹萌生时具有局限性。在连续损伤力学中，介质被假设是连续分布的。人们通过定义损伤变量来描述介质的破坏程度，但这种唯象的处理与“断裂是一种非连续变形”的客观事实不符。
对于结构及其材料的破坏问题，目前国际上新兴的近场动力学研究另辟蹊径，以非局部相互作用的思想为基础建立起一套新型的固体力学理论体系。近场动力学在物体上有限距离内的非接触物质点之间定义了“键”，物质点通过键进行相互作用。基于近场动力学的键，固体力学的三组控制方程被重新构建：其几何方程描述了两物质点之间键的变形程度。在该方程中不需要求导运算，因此放宽了对物体变形的连续性要求；其本构方程描述了键上力与变形之间的关系。通过定义键的断裂准则，判断计算过程中键断裂与否。计算过程中断键数量逐渐增加可以用于描述裂纹自发地成核和扩展过程；其力的平衡方程表现为积分方程，因此既可用于描述连续变形也可用于描述非连续变形，扩大了方程的应用范围和求解空间。近场动力学能够统一描述连续和非连续变形，能够用于预测结构及其材料在外部载荷作用下经变形、损伤、断裂到完全破坏的全使役过程。
鉴于近场动力学的上述优点，结合当今快速发展的计算机技术，如何在计算机上利用近场动力学模型高效地实现结构件的损伤分析具有重要价值。针对近场动力学模型的数值实现，目前有两种较为常见的数值计算方法：
(1)基于离散粒子的无网格方法
在该方法中，参考构型被离散为若干个具有一定体积的粒子。并且，采用对有限个粒子之间相互作用力进行求和的方式近似力的平衡方程中的积分运算。由此，可以直接对力平衡方程进行离散。这种离散求和的方式计算精度相对较低，需要采用大量的粒子进行计算，从而极大地增加了计算量。此外，在计算固体力学领域，现有的计算机辅助工程(CAE)软件多是基于有限单元方法开发的，因此粒子类方法与这些CAE软件兼容性差，限制了近场动力学模型的工程应用。
(2)基于连续单元的有限元法
在该方法中，参考构型被离散为若干个有限单元，单元之间相互连接且互不重叠，相邻单元的公共节点被称为网格节点。基于插值技术，单元上的位移场可以由网格节点处的位移值近似插值表示，并利用插值出的位移场对原问题进行近似表示。在经典的有限元法中，可以由单元节点位移表示出每个有限单元的单元总势能，并生成单元刚度矩阵，然后通过一定的规则进行组装即可得到总体刚度矩阵，进而建立关于总体节点位移的代数平衡方程。对于近场动力学模型，由于物质点间相互作用的非局部性，无法针对每个有限单元表示出其总势能。因此近场动力学模型的有限元实现过程与经典有限元法的计算流程差别显著，且更为复杂。当前工程中常用的有限元软件平台大多是基于经典有限元法的计算流程设计的，这导致近场动力学模型的有限元实现难以与已有的有限元软件平台兼容，限制了近场动力学模型的大规模工程应用。
发明内容
本发明的目的是解决现有的近场动力学模型数值实现方法存在的不足，利用近场动力学模型对损伤破坏模拟的能力，提供一种新的数值计算技术—近场有限元法，用于在商用有限元软件平台实现结构的损伤分析。具体是通过构造新型的近场单元，从而构建基于商用有限元软件平台实现结构损伤分析的近场有限元法。
为达到上述目的，本发明采用的技术方案如下：
结构损伤分析的近场有限元法，包括以下步骤：
(A1)对目标结构Ω进行几何建模，并按网格剖分尺寸h生成传统的有限单元网格，得到m个有限单元，m的值由目标结构Ω的几何尺寸和采用的网格剖分尺寸h共同确定；
(A2)对任意两个有限单元Ωj和Ωi，若Ωj和Ωi之间的距离dji满足0≤dji≤δh，则将Ωj和Ωi组合成一个近场单元(本发明的近场单元由且仅由两个有限单元组合而成，并且对组成该近场单元的两个有限单元的维数、形状、节点数等特征无特殊要求)，其中，j＝1,2,…,m，i＝1,2,…,m，m为有限单元的总个数，0≤δ≤100；
(A3)基于有限单元网格，计算总体载荷向量F，公式如下：
式中，m为有限单元的总个数；Gi为有限单元Ωi的节点自由度的转换矩阵(di＝Gid，di为有限单元Ωi的节点位移向量，d为总体节点位移向量)；·T表示矩阵的转置；Fi为有限单元Ωi的单元载荷向量；Ni(x)为组成的有限单元Ωi的形函数矩阵；x为组成/>的有限单元Ωi内的点；b(x)为外体力场在点x处的外体力向量；dVx为积分微元；[·]表示矩阵；ni为有限单元Ωi的节点个数；Nα(x)为有限单元Ωi的第α个节点的形函数，α＝1,2,…,ni；
(A4)基于近场单元网格，计算总体刚度矩阵公式如下：
式中，为近场单元/>的总个数；/>为近场单元/>的节点自由度的转换矩阵(为近场单元/>的节点位移向量，d为总体节点位移向量)；·T表示矩阵的转置；/>为近场单元/>的单元刚度矩阵；/>为近场单元/>的形函数差矩阵；x′为组成的有限单元Ωj内的点；x为组成/>的有限单元Ωi内的点；D(ξ)为微模量矩阵；dVx′为积分微元；dVx为积分微元；[·]表示矩阵；Nj(x′)为组成/>的有限单元Ωj的形函数矩阵；Ni(x)为组成/>的有限单元Ωi的形函数矩阵；c0(ξ)为键ξ的微模量系数，c0(ξ)的值与目标结构Ω的材料性质有关；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数，其中0表示所述键ξ断裂，1表示所述键ξ未断裂；||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；ξ1、ξ2、ξ3为键ξ向量的分量；
(A5)求解得到总体节点位移向量d，公式如下：
(A6)判断当前所给定的载荷下计算结果是否收敛，如果否，则返回步骤(A4)；反之，则进入步骤(A7)；
判断当前所给定的载荷下计算结果是否收敛的方法为：判断本次执行步骤(A5)得到的总体节点位移向量dt与上一次执行步骤(A5)得到的总体节点位移向量dt-1是否满足：||dt-dt-1||/||dt||≤ε，其中ε为给定的误差限，10-8h≤ε≤10-1h(h为网格剖分尺寸)，如果是，则收敛；反之，则不收敛；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：判断有无新增断键，如果无新增断键，则收敛；反之，则不收敛；
判断键是否断裂的方法为：如果任一近场单元中任一根键ξ的伸长率s大于某一给定的临界伸长率scrit(取值范围为(-1,100])，则该键断裂，其中scrit与目标结构Ω的材料性质有关；反之，则该键未断裂；s的计算公式如下：
ui(x)＝Ni(x)di；
uj(x′)＝Nj(x′)dj；
式中，||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；x′为组成的有限单元Ωj内的点；x为组成/>的有限单元Ωi内的点；uj(x′)为任一有限单元Ωj中任一点x′处的位移向量；ui(x)为任一有限单元Ωi中任一点x处的位移向量；Ni(x)为组成/>的有限单元Ωi的形函数矩阵；di为有限单元Ωi的节点位移向量；Nj(x′)为组成/>的有限单元Ωj的形函数矩阵；dj为有限单元Ωj的节点位移向量；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：先判断有无新增断键，如果无，则不更新函数μ(ξ,t)；反之，则更新函数μ(ξ,t)，再重新计算总体刚度矩阵记更新后的总体刚度矩阵为/>如果/>(dt为本次执行步骤(A5)得到的总体节点位移向量，F为总体载荷向量，ε为给定的误差限)，则收敛；反之，则不收敛；
(A7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(A3)；反之，则进入步骤(A8)；
(A8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
作为优选的技术方案：
如上所述的结构损伤分析的近场有限元法，步骤(A2)中，Ωj和Ωi之间的距离dji为Ωj和Ωi的形心距离，δ＝3；或者，Ωj和Ωi之间的距离dji为Ωj和Ωi的最接近点的距离，0≤δ≤100；或者，Ωj和Ωi之间的距离dji为Ωj和Ωi的最远离点的距离，0≤δ≤100。
如上所述的结构损伤分析的近场有限元法，步骤(A8)中，每一计算步t时目标结构Ω中任一点x处的等效损伤dξ(x,t)的计算公式如下：
式中，Hδ(x)为点x的近场邻域；表示所有Ωj的并集；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数；ωcrit为近场动力学键的临界断裂能；dVx′为积分微元；x′为组成/>的有限单元Ωj内的点；x为组成/>的有限单元Ωi内的点；||·||表示计算向量的长度；0≤δ≤100；h为网格剖分尺寸；c0(ξ)为所述键ξ的微模量系数，c0(ξ)的值与目标结构Ω的材料性质有关；scrit为给定的临界伸长率；ξ为近场动力学键。
本发明的上述结构损伤分析的近场有限元法与传统有限元法具有相同的计算流程，能够与商用有限元软件兼容。
本发明还提供如上任一项所述的结构损伤分析的近场有限元法在商用软件中的实现方法，利用商用有限元软件(例如ANSYS、ABAQUS、MSC Nastran、MSC Marc、ADINA等)的二次开发功能，在商用有限元软件中执行结构损伤分析的近场有限元法，包括以下操作：
(B1)由有限单元网格数据生成近场单元网格数据，并输入到商用有限元软件中；
(B2)将步骤(A4)中近场单元的单元刚度矩阵/>的计算公式编写为商用有限元软件的单元开发子程序，并嵌入到该商用有限元软件中用于完成目标结构Ω的单元刚度矩阵/>的计算；子程序中关于近场单元/>的单元刚度矩阵/>计算的积分公式采用数值积分公式；在组成近场单元/>的有限单元Ωj和Ωi内各自独立地选取积分点，并且Ωj和Ωi中的任一对积分点构成一根近场动力学键；
(B3)设定步骤(A6)中描述近场动力学键是否断裂的函数μ(ξ,t)的更新方式为：每次完成步骤(A5)后，根据得到的总体节点位移向量d，在单元开发子程序内部判断近场动力学键是否断裂，并更新函数μ(ξ,t)。
作为优选的技术方案：
如上所述的实现方法，(B1)中，设有限单元Ωj的节点的全局节点编号为设有限单元Ωi的节点的全局节点编号为/>若Ωj和Ωi之间的距离dji满足0≤dji≤δh(0≤δ≤100，h为网格剖分尺寸)，则将Ωj和Ωi组合成一个近场单元/>的节点的全局节点编号取为其中/>为近场单元/>的节点个数。
如上所述的实现方法，具体步骤如下：
(C0)按(B2)和(B3)的要求编写计算近场单元的单元刚度矩阵/>的单元开发子程序，并嵌入到商用有限元软件中以供使用；
(C1)在计算机中对目标结构Ω进行几何建模，并设定网格剖分尺寸h的取值后，按网格剖分尺寸h生成传统的有限单元网格，得到m个有限单元；
(C2)确定距离dji的计算方式以及δ的取值后，按照步骤(A2)和(B1)中所述的方式生成近场单元网格，得到个近场单元；
(C3)在商用有限元软件中，输入目标结构Ω所受到的载荷信息；
(C4)调用步骤(C0)中编写的计算近场单元的单元刚度矩阵/>的单元开发子程序对目标结构Ω进行损伤分析，商用有限元软件自动完成总体刚度矩阵/>的集成；例如：ANSYS中使用FORTRAN语言编写UserElem子程序，通过编译功能，将子程序嵌入到ANSYS软件中，调用此编译过的ANSYS软件执行计算；ABAQUS中使用FORTRAN语言编写UEL子程序，使用时直接调用这段子程序，ABAQUS软件在执行计算时能够自动完成编译并使用用户编写的子程序进行计算；
(C5)商用有限元软件根据步骤(C3)中输入的载荷信息和步骤(C4)中得到的总体刚度矩阵自动完成步骤(A5)中所述的总体节点位移向量d的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据步骤(C5)中得到的总体节点位移向量d，判断当前所给定的载荷下计算结果是否收敛，如果是，则进入步骤(C7)；反之，则返回步骤(C4)；
(C7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(C3)；反之，则进入步骤(C8)；
(C8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
有益效果
本发明公开的结构损伤分析的近场有限元法，与商用有限元软件平台兼容性好，在不改变有限元软件底层架构与核心求解代码的前提下易于将近场动力学求解代码嵌入到软件平台中；便于在基于经典有限元计算结果的基础上进一步利用近场动力学模型进行结构损伤破坏分析；数值计算过程简单明确，方便编程实现或借助有限元软件平台做二次开发应用。
附图说明
图1为有限单元和近场单元的示意图；
图2为基于有限单元网格生成近场单元的示意图；
图3为双开槽平板的几何尺寸与预设的最大载荷条件示意图；
图4为非结构化四节点四边形有限元网格示意图；
图5为双开槽平板结构件破坏的计算结果：(a)为uy＝2ux＝0.35mm时的等效损伤云图；(b)为uy＝2ux＝0.5mm时的等效损伤云图；
图6为双开槽平板结构件X方向位移的计算结果：(a)为uy＝2ux＝0.35mm时的X方向位移云图；(b)为uy＝2ux＝0.5mm时的X方向位移云图；
图7为双开槽平板结构件Y方向位移的计算结果：(a)为uy＝2ux＝0.35mm时的Y方向位移云图；(b)为uy＝2ux＝0.5mm时的Y方向位移云图。
具体实施方式
下面结合具体算例，给出具体实施方式，以进一步阐述本发明。应理解，该实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
结构损伤分析的近场有限元法，具体步骤如下：
(A1)对目标结构Ω进行几何建模，并按网格剖分尺寸h生成传统的有限单元网格，得到m个有限单元，m的值由目标结构Ω的几何尺寸和采用的网格剖分尺寸h共同确定；
(A2)对任意两个有限单元Ωj和Ωi(如图1所示)，若Ωj和Ωi之间的距离dji满足0≤dji≤δh，则将Ωj和Ωi组合成一个近场单元(如图1所示)，其中，j＝1,2,…,m，i＝1,2,…,m，m为有限单元的总个数，0≤δ≤100；其中，Ωj和Ωi之间的距离dji为Ωj和Ωi的形心距离，δ＝3；或者，Ωj和Ωi之间的距离dji为Ωj和Ωi的最接近点的距离，0≤δ≤100；或者，Ωj和Ωi之间的距离dji为Ωj和Ωi的最远离点的距离，0≤δ≤100；
(A3)基于有限单元网格，计算总体载荷向量F，公式如下：
式中，m为有限单元的总个数；Gi为有限单元Ωi的节点自由度的转换矩阵(di＝Gid，di为有限单元Ωi的节点位移向量，d为总体节点位移向量)；·T表示矩阵的转置；Fi为有限单元Ωi的单元载荷向量；di为有限单元Ωi的节点位移向量；d为总体节点位移向量；Ni(x)为组成的有限单元Ωi的形函数矩阵；x为组成/>的有限单元Ωi内的点；b(x)为外体力场在点x处的外体力向量；dVx为积分微元；[·]表示矩阵；ni为有限单元Ωi的节点个数；Nα(x)为有限单元Ωi的第α个节点的形函数，α＝1,2,…,ni；
(A4)基于近场单元网格，计算总体刚度矩阵公式如下：
式中，为近场单元/>的总个数；/>为近场单元/>的节点自由度的转换矩阵(为近场单元/>的节点位移向量，d为总体节点位移向量)；·T表示矩阵的转置；/>为近场单元/>的单元刚度矩阵；/>为近场单元/>的节点位移向量；d为总体节点位移向量；/>为近场单元/>的形函数差矩阵；x′为组成/>的有限单元Ωj内的点；x为组成/>的有限单元Ωi内的点；D(ξ)为微模量矩阵；dVx′为积分微元；dVx为积分微元；[·]表示矩阵；Nj(x′)为组成/>的有限单元Ωj的形函数矩阵；Ni(x)为组成/>的有限单元Ωi的形函数矩阵；c0(ξ)为所述键ξ的微模量系数，c0(ξ)的值与目标结构Ω的材料性质有关；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数，其中0表示所述键ξ断裂，1表示所述键ξ未断裂；|·|表示计算向量的长度；ξ＝x′-x为近场动力学键；ξ1、ξ2、ξ3为键ξ向量的分量；
(A5)求解得到总体节点位移向量d，公式如下：
(A6)判断当前所给定的载荷下计算结果是否收敛，如果否，则返回步骤(A4)；反之，则进入步骤(A7)；
判断当前所给定的载荷下计算结果是否收敛的方法为：判断本次执行步骤(A5)得到的总体节点位移向量dt与上一次执行步骤(A5)得到的总体节点位移向量dt-1是否满足：||dt-dt-1||/||dt||≤ε，其中ε为给定的误差限，10-8h≤ε≤10-1h(h为网格剖分尺寸)，如果是，则收敛；反之，则不收敛；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：判断有无新增断键，如果无新增断键，则收敛；反之，则不收敛；
判断键是否断裂的方法为：如果任一近场单元中任一根键ξ的伸长率s大于某一给定的临界伸长率scrit(取值范围为(-1,100])，则该键断裂，其中scrit与目标结构Ω的材料性质有关；反之，则该键未断裂；s的计算公式如下：
ui(x)＝Ni(x)di；
uj(x′)＝Nj(x′)dj；
式中，||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；x′为组成的有限单元Ωj内的点；x为组成/>的有限单元Ωi内的点；uj(x′)为任一有限单元Ωj中任一点x′处的位移向量；ui(x)为任一有限单元Ωi中任一点x处的位移向量；Ni(x)为组成/>的有限单元Ωi的形函数矩阵；di为有限单元Ωi的节点位移向量；Nj(x′)为组成/>的有限单元Ωj的形函数矩阵；dj为有限单元Ωj的节点位移向量；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：先判断有无新增断键，如果无，则不更新函数μ(ξ,t)；反之，则更新函数μ(ξ,t)，再重新计算总体刚度矩阵记更新后的总体刚度矩阵为/>如果/>(dt为本次执行步骤(A5)得到的总体节点位移向量，F为总体载荷向量，ε为给定的误差限)，则收敛；反之，则不收敛；
(A7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(A3)；反之，则进入步骤(A8)；
(A8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图；其中，每一计算步t时目标结构Ω中任一点x处的等效损伤dξ(x,t)的计算公式如下：
式中，Hδ(x)为点x的近场邻域；表示所有Ωj的并集；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数；ωcrit为近场动力学键的临界断裂能；dVx′为积分微元；x′为组成/>的有限单元Ωj内的点；x为组成/>的有限单元Ωi内的点；||·||表示计算向量的长度；0≤δ≤100；h为网格剖分尺寸；c0(ξ)为所述键ξ的微模量系数，c0(ξ)的值与目标结构Ω的材料性质有关；scrit为给定的临界伸长率；ξ为近场动力学键。
结构损伤分析的近场有限元法在商用软件中的实现方法，利用商用有限元软件(例如ANSYS、ABAQUS、MSC Nastran、MSC Marc、ADINA等)的二次开发功能，在商用有限元软件中执行结构损伤分析的近场有限元法，包括以下操作：
(B1)由有限单元网格数据生成近场单元网格数据，并输入到商用有限元软件中；具体地，设有限单元Ωj的节点的全局节点编号为设有限单元Ωi的节点的全局节点编号为/>若Ωj和Ωi之间的距离dji满足0≤dji≤δh(0≤δ≤100，h为网格剖分尺寸)，则将Ωj和Ωi组合成一个近场单元/>的节点的全局节点编号取为/>其中/>为近场单元/>的节点个数；
(B2)将步骤(A4)中近场单元的单元刚度矩阵/>的计算公式编写为商用有限元软件的单元开发子程序，并嵌入到该商用有限元软件中用于完成目标结构Ω的单元刚度矩阵/>的计算；子程序中关于近场单元/>的单元刚度矩阵/>计算的积分公式采用数值积分公式；在组成近场单元/>的有限单元Ωj和Ωi内各自独立地选取积分点，并且Ωj和Ωi中的任一对积分点构成一根近场动力学键；
(B3)设定步骤(A6)中描述近场动力学键是否断裂的函数μ(ξ,t)的更新方式为：每次完成步骤(A5)后，根据得到的总体节点位移向量d，在单元开发子程序内部判断近场动力学键是否断裂，并更新函数μ(ξ,t)。
结构损伤分析的近场有限元法在商用软件中的实现方法的具体步骤如下：
(C0)按(B2)和(B3)的要求编写计算近场单元的单元刚度矩阵/>的单元开发子程序，并嵌入到商用有限元软件中以供使用；
(C1)在计算机中对目标结构Ω进行几何建模，并设定网格剖分尺寸h的取值后，按网格剖分尺寸h生成传统的有限单元网格，得到m个有限单元；
(C2)确定距离dji的计算方式以及δ的取值后，按照步骤(A2)和(B1)中所述的方式生成近场单元网格，得到个近场单元；
(C3)在商用有限元软件中，输入目标结构Ω所受到的载荷信息；
(C4)调用步骤(C0)中编写的计算近场单元的单元刚度矩阵/>的单元开发子程序对目标结构Ω进行损伤分析，商用有限元软件自动完成总体刚度矩阵/>的集成；
(C5)商用有限元软件根据步骤(C3)中输入的载荷信息和步骤(C4)中得到的总体刚度矩阵K，自动完成步骤(A5)中所述的总体节点位移向量d的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据步骤(C5)中得到的总体节点位移向量d，判断当前所给定的载荷下计算结果是否收敛，如果是，则进入步骤(C7)；反之，则返回步骤(C4)；
(C7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(C3)；反之，则进入步骤(C8)；
(C8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
现结合具体使用过程进行举例说明：
问题设置：考虑一个二维问题，具体是一个具有两条水平开槽的平板(双开槽平板)，其几何尺寸和预设的最大载荷条件如图3所示，预设的最大载荷条件为水平方向位移ux＝0.25mm，竖直方向uy＝0.5mm；本使用过程不考虑外力载荷的作用，目标结构的材料的杨氏模量和泊松比分别被设置为E＝30GPa和v＝1/3；二维时，基于有限单元网格生成近场单元的示意图如图2所示；数值计算时，采用非结构化四节点四边形有限单元网格，具体如图4所示；
结构损伤分析的近场有限元法在商用软件(ANSYS)中的实现方法，具体步骤如下：
(C0)编写计算8节点近场单元的单元刚度矩阵的单元开发子程序，并嵌入到商用有限元软件中以供使用；其中，单元开发子程序中，函数μ(ξ,t)的更新方式为：每次完成步骤(A5)后，根据得到的总体节点位移向量d，在单元开发子程序内部判断近场动力学键是否断裂，并更新函数μ(ξ,t)；
(C1)在计算机中对目标结构进行几何建模，并按网格剖分尺寸h＝1.2生成有限单元网格，得到m＝27741个有限单元；
(C2)取dji为Ωj和Ωi的形心距离，并取δ＝3，按照步骤(A2)和步骤(B1)中所述的方式生成近场单元网格，得到m＝1085727个近场单元；
(C3)在商用有限元软件中，输入目标结构的初始载荷信息和预设的最大载荷信息；其中初始载荷信息为：水平方向位移ux＝0.025mm，竖直方向uy＝0.05mm；
(C4)调用步骤(C0)中编写的计算近场单元的单元刚度矩阵/>的单元开发子程序对目标结构进行损伤分析；商用有限元软件能够自动完成总体刚度矩阵/>的集成；
(C5)商用有限元软件能够根据步骤(C3)中输入的初始载荷信息和和步骤(C4)中得到的总体刚度矩阵自动完成步骤(A5)中所述的节点位移的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据(C5)中得到的总体节点位移向量d，判定当前载荷下的位移求解结果是否收敛，如果收敛，则进入步骤(C7)；否则，返回步骤(C4)；其中，本实施例中，采用是否有新增断键作为收敛准则；
(C7)判断当前所施加的载荷值是否超过预设最大载荷值；如果没有超过，返回步骤(C3)，在初始载荷信息的基础上增加载荷值(水平方向位移ux＝0.025mm，竖直方向uy＝0.05mm)得到新的初始载荷信息；反之，则进入步骤(C8)；
(C8)根据总体节点位移向量d和断键信息，输出目标结构在每一计算步的位移云图和等效损伤云图；当uy＝2ux＝0.35mm时，等效损伤云图如图5(a)所示，X方向位移云图如图6(a)所示，Y方向位移云图如图7(a)所示；当uy＝2ux＝0.5mm时，等效损伤云图如图5(b)所示，X方向位移云图如图6(b)所示，Y方向位移云图如图7(b)所示。

Claims (6)

1.结构损伤分析的近场有限元法，其特征在于，包括以下步骤：
(A1)对目标结构Ω进行几何建模，并按网格剖分尺寸h生成有限单元网格，得到m个有限单元；
(A2)对任意两个有限单元Ωj和Ωi，若Ωj和Ωi之间的距离dji满足0≤dji≤δh，则将Ωj和Ωi组合成一个近场单元其中，j＝1，2，…，m，i＝1，2，…，m，0≤δ≤100；
(A3)基于有限单元网格，计算总体载荷向量F，公式如下：
式中，Gi为有限单元Ωi的节点自由度的转换矩阵；Fi为有限单元Ωi的单元载荷向量；dVx为积分微元；Ni(x)为组成的有限单元Ωi的形函数矩阵；x为组成/>的有限单元Ωi内的点；b(x)为外体力场在点x处的外体力向量；ni为有限单元Ωi的节点个数；Nα(x)为有限单元Ωi的第α个节点的形函数，α＝1，2，…，ni；
(A4)基于近场单元网格，计算总体刚度矩阵公式如下：
式中，为近场单元/>的总个数；/>为近场单元/>的节点自由度的转换矩阵；/>为近场单元/>的单元刚度矩阵；dVx′为积分微元；/>为近场单元/>的形函数差矩阵；x′为组成/>的有限单元Ωj内的点；D(ξ)为微模量矩阵；Nj(x′)为组成/>的有限单元Ωj的形函数矩阵；c0(ξ)为键ξ的微模量系数，μ(ξ，t)为相关于所述键ξ和计算步t的取值为0或1的函数，其中0表示所述键ξ断裂，1表示所述键ξ未断裂；||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；ξ1、ξ2、ξ3为键ξ向量的分量；
(A5)求解得到总体节点位移向量d，公式如下：
(A6)判断当前所给定的载荷下计算结果是否收敛，如果否，则返回步骤(A4)；反之，则进入步骤(A7)；
判断当前所给定的载荷下计算结果是否收敛的方法为：判断本次执行步骤(A5)得到的总体节点位移向量dt与上一次执行步骤(A5)得到的总体节点位移向量dt-1是否满足：||dt-dt-1||/||dt||≤ε，10-8h≤ε≤10-1h，如果是，则收敛；反之，则不收敛；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：判断有无新增断键，如果无新增断键，则收敛；反之，则不收敛；
判断键是否断裂的方法为：如果任一近场单元中任一根键ξ的伸长率s大于某一给定的临界伸长率scrit，则该键断裂，其中scrit与目标结构Ω的材料性质有关；反之，则该键未断裂；s的计算公式如下：
ui(x)＝Ni(x)di；
uj(x′)＝Nj(x′)dj；
式中，uj(x′)为任一有限单元Ωj中任一点x′处的位移向量；ui(x)为任一有限单元Ωi中任一点x处的位移向量；di为有限单元Ωi的节点位移向量；dj为有限单元Ωj的节点位移向量；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：先判断有无新增断键，如果无，则不更新函数μ(ξ，t)；反之，则更新函数μ(ξ，t)，再重新计算总体刚度矩阵记更新后的总体刚度矩阵为/>如果/>则收敛；反之，则不收敛；
(A7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(A3)；反之，则进入步骤(A8)；
(A8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
2.根据权利要求1所述的结构损伤分析的近场有限元法，其特征在于，步骤(A2)中，Ωj和Ωi之间的距离dji为Ωj和Ωi的形心距离，δ＝3；或者，Ωj和Ωi之间的距离dji为Ωj和Ωi的最接近点的距离，o≤δ≤100；或者，Ωj和Ωi之间的距离dji为Ωj和Ωi的最远离点的距离，o≤δ≤100。
3.根据权利要求1所述的结构损伤分析的近场有限元法，其特征在于，步骤(A8)中，每一计算步t时目标结构Ω中任一点x处的等效损伤dξ(x，t)的计算公式如下：
式中，Hδ(x)为点x的近场邻域；表示所有Ωj的并集；ωcrit为近场动力学键的临界断裂能。
4.如权利要求1～3任一项所述的结构损伤分析的近场有限元法在商用软件中的实现方法，其特征在于，利用商用有限元软件的二次开发功能，在商用有限元软件中执行结构损伤分析的近场有限元法，包括以下操作：
(B1)由有限单元网格数据生成近场单元网格数据，并输入到商用有限元软件中；
(B2)将步骤(A4)中近场单元的单元刚度矩阵/>的计算公式编写为商用有限元软件的单元开发子程序，并嵌入到该商用有限元软件中；子程序中关于近场单元/>的单元刚度矩阵/>计算的积分公式采用数值积分公式；在组成近场单元/>的有限单元Ωj和Ωi内各自独立地选取积分点，并且Ωj和Ωi中的任一对积分点构成一根近场动力学键；
(B3)设定步骤(A6)中函数μ(ξ，t)的更新方式为：每次完成步骤(A5)后，根据得到的总体节点位移向量d，在单元开发子程序内部判断近场动力学键是否断裂，并更新函数μ(ξ，t)。
5.根据权利要求4所述的实现方法，其特征在于，(B1)中，设有限单元Ωj的节点的全局节点编号为设有限单元Ωi的节点的全局节点编号为若Ωj和Ωi之间的距离dji满足0≤dji≤δh，则将Ωj和Ωi组合成一个近场单元/> 的节点的全局节点编号取为其中/>为近场单元/>的节点个数。
6.根据权利要求4或5所述的实现方法，其特征在于，具体步骤如下：
(C0)按(B2)和(B3)的要求编写计算近场单元的单元刚度矩阵/>的单元开发子程序，并嵌入到商用有限元软件中以供使用；
(C1)在计算机中对目标结构Ω进行几何建模，并设定网格剖分尺寸h的取值后，按网格剖分尺寸h生成有限单元网格，得到m个有限单元；
(C2)确定距离dji的计算方式以及δ的取值后，按照步骤(A2)和(B1)中所述的方式生成近场单元网格，得到个近场单元；
(C3)在商用有限元软件中，输入目标结构Ω所受到的载荷信息；
(C4)调用步骤(C0)中编写的计算近场单元的单元刚度矩阵/>的单元开发子程序对目标结构Ω进行损伤分析，商用有限元软件自动完成总体刚度矩阵/>的集成；
(C5)商用有限元软件根据步骤(C3)中输入的载荷信息和步骤(C4)中得到的总体刚度矩阵自动完成步骤(A5)中所述的总体节点位移向量d的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据步骤(C5)中得到的总体节点位移向量d，判断当前所给定的载荷下计算结果是否收敛，如果是，则进入步骤(C7)；反之，则返回步骤(C4)；
(C7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(C3)；反之，则进入步骤(C8)；
(C8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
CN202110885638.8A 2021-08-03 2021-08-03 结构损伤分析的近场有限元法及在商用软件中的实现方法 Active CN113705040B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN202110885638.8A CN113705040B (zh) 2021-08-03 2021-08-03 结构损伤分析的近场有限元法及在商用软件中的实现方法
US18/024,780 US20230325560A1 (en) 2021-08-03 2022-06-22 A peridynamics-based finite element method (perifem) for the analysis of structural damage and its implementation in commercial software
PCT/CN2022/100311 WO2023011029A1 (zh) 2021-08-03 2022-06-22 结构损伤分析的近场有限元法及在商用软件中的实现方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110885638.8A CN113705040B (zh) 2021-08-03 2021-08-03 结构损伤分析的近场有限元法及在商用软件中的实现方法

Publications (2)

Publication Number Publication Date
CN113705040A CN113705040A (zh) 2021-11-26
CN113705040B true CN113705040B (zh) 2024-03-22

Family

ID=78651448

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110885638.8A Active CN113705040B (zh) 2021-08-03 2021-08-03 结构损伤分析的近场有限元法及在商用软件中的实现方法

Country Status (3)

Country Link
US (1) US20230325560A1 (zh)
CN (1) CN113705040B (zh)
WO (1) WO2023011029A1 (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113705040B (zh) * 2021-08-03 2024-03-22 大连理工大学 结构损伤分析的近场有限元法及在商用软件中的实现方法
CN114491944B (zh) * 2021-12-23 2022-10-25 贵州大学 一种桥梁损伤分析方法、分析报警系统及装置
CN115828585B (zh) * 2022-11-30 2023-07-28 四川大学 基于GiD-Code_bright软件的批量修改单元材料方法
CN116364218B (zh) * 2023-03-30 2023-09-08 沈阳工业大学 基于近场动力学岩石材料率效应本构模型构建与实施方法
CN116050027B (zh) * 2023-03-30 2023-07-07 陕西空天信息技术有限公司 叶轮叶片结构静态分析方法、存储介质和电子设备
CN117688807A (zh) * 2023-12-01 2024-03-12 河海大学 基于间接力键有限元法的结构变形与损伤分析方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110457790A (zh) * 2019-07-26 2019-11-15 顾鑫 用于结构变形分析的近场动力学非连续伽辽金有限元方法
CN111159951A (zh) * 2019-12-31 2020-05-15 哈尔滨工业大学（深圳） 一种基于abaqus有限元与边界元的耦合方法
CN112446163A (zh) * 2020-11-24 2021-03-05 西安交通大学 基于参数化水平集的能量有限元拓扑优化方法
CN112966421A (zh) * 2021-03-16 2021-06-15 昆明理工大学 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9405867B2 (en) * 2012-06-07 2016-08-02 Dassault Systemes Simulia Corp. Hydraulic fracture simulation with an extended finite element method
CN111339594B (zh) * 2020-02-26 2022-12-16 河北工业大学 基于dic技术的近场动力学参数实验反演系统及使用方法
CN111814310B (zh) * 2020-06-11 2022-09-13 大连理工大学 强度准则驱动的近场动力学模型预测结构破坏的方法
CN113705040B (zh) * 2021-08-03 2024-03-22 大连理工大学 结构损伤分析的近场有限元法及在商用软件中的实现方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110457790A (zh) * 2019-07-26 2019-11-15 顾鑫 用于结构变形分析的近场动力学非连续伽辽金有限元方法
CN111159951A (zh) * 2019-12-31 2020-05-15 哈尔滨工业大学（深圳） 一种基于abaqus有限元与边界元的耦合方法
CN112446163A (zh) * 2020-11-24 2021-03-05 西安交通大学 基于参数化水平集的能量有限元拓扑优化方法
CN112966421A (zh) * 2021-03-16 2021-06-15 昆明理工大学 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
破片冲击作用下舰船复合材料结构损伤的近场动力学模拟;杨娜娜;赵天佑;陈志鹏;武国勋;姚熊亮;;爆炸与冲击(第02期);67-77 *

Also Published As

Publication number Publication date
US20230325560A1 (en) 2023-10-12
CN113705040A (zh) 2021-11-26
WO2023011029A1 (zh) 2023-02-09

Similar Documents

Publication Publication Date Title
CN113705040B (zh) 结构损伤分析的近场有限元法及在商用软件中的实现方法
Rannou et al. A local multigrid X‐FEM strategy for 3‐D crack propagation
CN111859766A (zh) 可变计算域的拉格朗日积分点有限元数值仿真系统及方法
Wang et al. Simulations of the dynamics and interaction between a floating structure and a near-field explosion bubble
CN111125963A (zh) 基于拉格朗日积分点有限元的数值仿真系统及方法
Jiang Nonlinear thermomechanical analysis of structures using OpenSees
Wang et al. Reduced order modeling with local enrichment for the nonlinear geometric response of a cracked panel
Zhang et al. Investigations on the hydroelastic slamming of deformable wedges by using the smoothed particle element method
Burczyński et al. Intelligent optimal design of spatial structures
CN113761760A (zh) 工程尺度岩体破裂全过程模拟的pd-fem数值计算方法及系统
US20230297737A1 (en) Modelling an object, determination of load capacity, improvement of the design and generating component, and system
Su et al. Computational morphogenesis of free‐form grid structures with Voronoi diagram
Chen et al. Contact analysis within the bi-potential framework using cell-based smoothed finite element method
Averseng et al. Interactive dynamic design and analysis of tensegrity systems
Xu et al. An integrated method of CAD, CAE and multi-objective optimization
Chu Evolutionary structural optimization method for systems with stiffness and displacement constraints
CN112989661A (zh) 一种联合拓扑优化与形状优化的水下结构设计方法
CN114036689B (zh) 一种基于迭代的构件强度应力优化方法
Jäger et al. Simulation of masonry in ANSYS and LS-DYNA the features and challenges
Song et al. Geometrically nonlinear analysis of Reissner–Mindlin plates using multi-patch isogeometric analysis based on Nitsche’s method
Tae et al. Computational differential algebraic equation framework and multi spatial and time discretizations preserving consistent second-order accuracy: nonlinear dynamics
WO2024103241A1 (zh) 一种基于物质点法的软体机器人仿真方法
Liu et al. Gradient Smoothing Methods with Programming: Applications to Fluids and Landslides
Zhu et al. Development of a Data-driven Turbulence Model for 3d Thermal Stratification Simulation during Reactor Transients
CN117688807A (zh) 基于间接力键有限元法的结构变形与损伤分析方法

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