WO2023011029A1 - 结构损伤分析的近场有限元法及在商用软件中的实现方法 - Google Patents
结构损伤分析的近场有限元法及在商用软件中的实现方法 Download PDFInfo
- Publication number
- WO2023011029A1 WO2023011029A1 PCT/CN2022/100311 CN2022100311W WO2023011029A1 WO 2023011029 A1 WO2023011029 A1 WO 2023011029A1 CN 2022100311 W CN2022100311 W CN 2022100311W WO 2023011029 A1 WO2023011029 A1 WO 2023011029A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- finite element
- bond
- field
- calculation
- unit
- Prior art date
Links
- 238000000034 method Methods 0.000 title claims abstract description 82
- 230000006378 damage Effects 0.000 title claims abstract description 54
- 238000004458 analytical method Methods 0.000 title claims abstract description 29
- 238000004364 calculation method Methods 0.000 claims abstract description 84
- 238000006073 displacement reaction Methods 0.000 claims abstract description 83
- 239000011159 matrix material Substances 0.000 claims abstract description 72
- 238000011161 development Methods 0.000 claims abstract description 22
- 230000006870 function Effects 0.000 claims description 41
- 239000000463 material Substances 0.000 claims description 15
- 230000010354 integration Effects 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 4
- 230000008569 process Effects 0.000 abstract description 13
- 238000010586 diagram Methods 0.000 description 13
- 239000002245 particle Substances 0.000 description 5
- 239000007787 solid Substances 0.000 description 5
- 230000003993 interaction Effects 0.000 description 3
- 241000682719 Adina Species 0.000 description 2
- 230000006911 nucleation Effects 0.000 description 2
- 238000010899 nucleation Methods 0.000 description 2
- 238000005381 potential energy Methods 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 230000017105 transposition Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000004141 dimensional analysis Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000002269 spontaneous effect Effects 0.000 description 1
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
- 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
Definitions
- the invention belongs to the field of computational mechanics, and in particular relates to a near-field finite element method for structural damage analysis and an implementation method in commercial software.
- Peridynamics defines "bonds" between non-contacting matter points on an object within a finite distance, through which matter points interact. Based on the peridynamic bond, the three sets of governing equations of solid mechanics are reconstructed: the geometric equation describes the degree of deformation of the bond between two material points. No derivation operation is required in this equation, so the continuity requirement on the deformation of the object is relaxed; its constitutive equation describes the relationship between the force on the bond and the deformation. By defining the bond breaking criterion, it is judged whether the bond is broken or not during the calculation process.
- the gradual increase in the number of interrupted bonds during the calculation process can be used to describe the spontaneous nucleation and growth process of cracks; the balance equation of its force is expressed as an integral equation, so it can be used to describe both continuous deformation and discontinuous deformation, expanding the application of the equation range and solution space.
- Peridynamics can uniformly describe continuous and discontinuous deformation, and can be used to predict the full service process of structures and their materials from deformation, damage, fracture to complete destruction under external loads.
- the reference configuration is discretized into several particles with a certain volume.
- the integral operation in the force balance equation is approximated by summing the interaction forces between a finite number of particles.
- the force balance equation can be discretized directly.
- the calculation accuracy of this discrete summation method is relatively low, and a large number of particles are required for calculation, which greatly increases the calculation amount.
- most of the existing computer-aided engineering (CAE) software is developed based on the finite element method, so the particle method has poor compatibility with these CAE software, which limits the engineering application of the peridynamic model.
- the reference configuration is discretized into several finite units, and the units are connected to each other without overlapping, and the common nodes of adjacent units are called grid nodes.
- the displacement field on the unit can be approximated by interpolation of the displacement values at the grid nodes, and the original problem can be approximated by using the interpolated displacement field.
- the total potential energy of each finite element can be expressed by the unit node displacement, and the element stiffness matrix can be generated, and then the overall stiffness matrix can be obtained by assembling certain rules, and then the overall node displacement can be established.
- the finite element implementation process of the peridynamic model is significantly different from the calculation process of the classical finite element method, and is more complicated.
- Most of the finite element software platforms commonly used in current engineering are designed based on the calculation process of the classical finite element method, which makes it difficult for the finite element realization of the peridynamic model to be compatible with the existing finite element software platforms, which limits the Large-scale engineering applications of models.
- the purpose of the present invention is to solve the deficiencies in the existing numerical realization methods of the peridynamic model, and to provide a new numerical calculation technology—proximate finite element method, by utilizing the ability of the peridynamic model to simulate damage and damage.
- the aim is to implement damage analysis of structures on a commercial finite element software platform. Specifically, by constructing a new type of near-field element, a near-field finite element method based on a commercial finite element software platform for structural damage analysis is constructed.
- (A1) Perform geometric modeling on the target structure ⁇ , and generate a traditional finite element grid according to the mesh size h to obtain m finite elements.
- the value of m is determined by the geometric size of the target structure ⁇ and the adopted grid
- the sub-size h is jointly determined;
- m is the total number of finite elements
- T represents the transpose of the matrix
- F i is the unit load vector of the finite element ⁇ i
- N i (x) is the composition
- the shape function matrix of the finite element ⁇ i of ; x is the composition
- the point in the finite element ⁇ i of ; b(x) is the external body force vector of the external body force field at point x; dV x is the integral microelement;
- [ ] represents the matrix
- n i is the number of nodes of the finite element ⁇ i ;
- N ⁇ (x) is the shape function of the ⁇ th node of the finite element ⁇ i ,
- near field unit the total number of near field unit
- D( ⁇ ) is the micromodulus matrix;
- dV x′ is the integral differential element;
- dV x is the integral differential element;
- [ ] represents the matrix;
- N j (x′) is the composition
- the shape function matrix of the finite element ⁇ j of ; N i (x) is the composition
- step (A6) Judging whether the calculation result converges under the given load at present, if not, then return to step (A4); otherwise, then enter step (A7);
- the method for judging whether the calculation results under the current given load are convergent is as follows: judging the overall node displacement vector d t obtained in this execution step (A5) and the overall node displacement vector d t-1 obtained in the previous execution step (A5) Is it satisfied:
- the method for judging whether the calculation result under the current given load is convergent is as follows: judge whether there is any new broken bond, if there is no new broken bond, it will converge; otherwise, it will not converge;
- the method for judging whether the bond is broken is: if any near-field unit If the elongation s of any bond ⁇ is greater than a given critical elongation s crit (value range is (-1,100]), the bond will be broken, where s crit is related to the material properties of the target structure ⁇ ; On the contrary, the bond is not broken; the calculation formula of s is as follows:
- the method for judging whether the calculation results converge under the current given load is: first judge whether there are new broken keys, if not, then do not update the function ⁇ ( ⁇ ,t); otherwise, update the function ⁇ ( ⁇ , t), and then recalculate the overall stiffness matrix Write down the updated overall stiffness matrix as if (d t is the overall node displacement vector obtained in this execution step (A5), F is the overall load vector, and ⁇ is a given error limit), then it converges; otherwise, it does not converge;
- step (A7) Judging whether the currently applied load value exceeds the preset maximum load value, if not, after increasing the load value, return to step (A3); otherwise, enter step (A8);
- the distance d ji between ⁇ j and ⁇ i is the distance of the closest point between ⁇ j and ⁇ i , 0 ⁇ 100; or, the distance d ji between ⁇ j and ⁇ i is the farthest point between ⁇ j and ⁇ i Point distance, 0 ⁇ 100.
- step (A8) the calculation formula of the equivalent damage d ⁇ (x,t) at any point x in the target structure ⁇ at each calculation step t is as follows:
- H ⁇ (x) is the near-field neighborhood of point x; Represents the union of all ⁇ j ; ⁇ ( ⁇ ,t) is a function of 0 or 1 related to the bond ⁇ and calculation step t; ⁇ crit is the critical fracture energy of the peridynamic bond; dV x ' is the integral microelement; x' is the composition The points in the finite element ⁇ j of ; x is the composition Points within the finite element ⁇ i of ;
- the near-field finite element method for structural damage analysis of the present invention has the same calculation process as the traditional finite element method, and is compatible with commercial finite element software.
- the present invention also provides the implementation method of the near-field finite element method of structural damage analysis as described in any one of the above in commercial software, utilizes the two-dimensional analysis of commercial finite element software (such as ANSYS, ABAQUS, MSC Nastran, MSC Marc, ADINA, etc.)
- step (B2) the near-field unit in step (A4)
- the element stiffness matrix of The calculation formula of is written as the element development subroutine of the commercial finite element software, and embedded into the commercial finite element software to complete the element stiffness matrix of the target structure ⁇
- the element stiffness matrix of The calculated integral formula adopts the numerical integral formula; in forming the near-field unit Integral points are independently selected in ⁇ j and ⁇ i of finite elements, and any pair of integration points in ⁇ j and ⁇ i constitutes a peridynamic bond;
- the global node number of the node with finite element ⁇ j is The global node number of the node with finite element ⁇ i is If the distance d ji between ⁇ j and ⁇ i satisfies 0 ⁇ d ji ⁇ h (0 ⁇ 100, h is the grid division size), then combine ⁇ j and ⁇ i into a near-field unit The global node number of the node is taken as in near field unit the number of nodes.
- (C0) Compile the calculation near-field unit according to the requirements of (B2) and (B3)
- the element stiffness matrix of The unit development subroutine of the unit is embedded into the commercial finite element software for use;
- step (C4) Call the calculated near-field unit written in step (C0)
- the element stiffness matrix of The unit development subroutine for the target structure ⁇ is used for damage analysis, and the commercial finite element software automatically completes the overall stiffness matrix
- step (C5) Commercial finite element software is based on the load information input in step (C3) and the overall stiffness matrix obtained in step (C4) Automatically complete the solution of the overall node displacement vector d described in step (A5);
- step (C6) After determining the judging condition of whether the calculation result converges under the current given load, according to the overall node displacement vector d obtained in the step (C5), judge whether the calculation result under the current given load converges, if yes, Then enter step (C7); otherwise, then return to step (C4);
- step (C7) Judging whether the currently applied load value exceeds the preset maximum load value, if not, after increasing the load value, return to step (C3); otherwise, enter step (C8);
- the near-field finite element method for structural damage analysis disclosed by the present invention has good compatibility with the commercial finite element software platform, and it is easy to embed the near-field dynamics solution code into the software platform without changing the underlying structure and core solution code of the finite element software Medium; it is convenient to further use the peridynamic model for structural damage analysis based on the classical finite element calculation results; the numerical calculation process is simple and clear, and it is convenient for programming or secondary development and application with the help of the finite element software platform.
- Fig. 1 is the schematic diagram of finite element and near-field element
- Fig. 2 is a schematic diagram of generating a near-field unit based on a finite element grid
- Fig. 3 is a schematic diagram of the geometric dimensions and preset maximum load conditions of a double-slotted flat plate
- Fig. 4 is a schematic diagram of an unstructured four-node quadrilateral finite element grid
- the near-field finite element method for structural damage analysis the specific steps are as follows:
- (A1) Perform geometric modeling on the target structure ⁇ , and generate a traditional finite element grid according to the mesh size h to obtain m finite elements.
- the value of m is determined by the geometric size of the target structure ⁇ and the adopted grid
- the sub-size h is jointly determined;
- m is the total number of finite elements
- T represents the transpose of the matrix;
- F i is the unit load vector of finite element ⁇ i ;
- d i is the nodal displacement vector of finite element ⁇ i ;
- d is the overall nodal displacement vector;
- N i (x) is the composition
- the shape function matrix of the finite element ⁇ i of ; x is the composition
- b(x) is the external body force vector of the external body force field at point x;
- dV x is the integral microelement;
- [ ] represents the matrix;
- n i is the number of nodes of the
- near field unit the total number of near field unit
- the transformation matrix of the nodal degrees of freedom ( near field unit The node displacement vector of , d is the overall node displacement vector); T represents the transposition of the matrix; near field unit The element stiffness matrix of ; near field unit The node displacement vector of ; d is the overall node displacement vector; near field unit The shape function difference matrix of ; x′ is the composition
- D( ⁇ ) is the micromodulus matrix;
- dV x′ is the integral differential element;
- dV x is the integral differential element; [ ] represents the matrix;
- N j (x′) is the composition
- the shape function matrix of the finite element ⁇ j of ; N i (x) is the composition
- the shape function matrix of the finite element ⁇ i ; c 0 ( ⁇ ) is the micromodulus coefficient of the bond ⁇ , and the value of c 0 ( ⁇ ) is related to the material properties of the target structure
- step (A6) Judging whether the calculation result converges under the given load at present, if not, then return to step (A4); otherwise, then enter step (A7);
- the method for judging whether the calculation results under the current given load are convergent is as follows: judging the overall node displacement vector d t obtained in this execution step (A5) and the overall node displacement vector d t-1 obtained in the previous execution step (A5) Is it satisfied:
- the method for judging whether the calculation result under the current given load is convergent is as follows: judge whether there is any new broken bond, if there is no new broken bond, it will converge; otherwise, it will not converge;
- the method for judging whether the bond is broken is: if any near-field unit If the elongation s of any bond ⁇ is greater than a given critical elongation s crit (value range is (-1,100]), the bond will be broken, where s crit is related to the material properties of the target structure ⁇ ; On the contrary, the bond is not broken; the calculation formula of s is as follows:
- the method for judging whether the calculation results converge under the current given load is: first judge whether there are new broken keys, if not, then do not update the function ⁇ ( ⁇ ,t); otherwise, update the function ⁇ ( ⁇ , t), and then recalculate the overall stiffness matrix Write down the updated overall stiffness matrix as if (d t is the overall node displacement vector obtained in this execution step (A5), F is the overall load vector, and ⁇ is a given error limit), then it converges; otherwise, it does not converge;
- step (A7) Judging whether the currently applied load value exceeds the preset maximum load value, if not, after increasing the load value, return to step (A3); otherwise, enter step (A8);
- H ⁇ (x) is the near-field neighborhood of point x; Represents the union of all ⁇ j ; ⁇ ( ⁇ ,t) is a function of 0 or 1 related to the bond ⁇ and calculation step t; ⁇ crit is the critical fracture energy of the peridynamic bond; dV x ' is the integral microelement; x' is the composition The points in the finite element ⁇ j of ; x is the composition Points within the finite element ⁇ i of ;
- the implementation method of near-field finite element method for structural damage analysis in commercial software using the secondary development function of commercial finite element software (such as ANSYS, ABAQUS, MSC Nastran, MSC Marc, ADINA, etc.), to execute in commercial finite element software Near-field finite element method for structural damage analysis, including the following operations:
- (B1) Generate near-field element grid data from finite element grid data and input it into commercial finite element software; specifically, the global node number of the node with finite element ⁇ j is The global node number of the node with finite element ⁇ i is If the distance d ji between ⁇ j and ⁇ i satisfies 0 ⁇ d ji ⁇ h (0 ⁇ 100, h is the grid division size), then combine ⁇ j and ⁇ i into a near-field unit The global node number of the node is taken as in near field unit the number of nodes;
- step (B2) the near-field unit in step (A4)
- the element stiffness matrix of The calculation formula of is written as the element development subroutine of the commercial finite element software, and embedded into the commercial finite element software to complete the element stiffness matrix of the target structure ⁇
- the element stiffness matrix of The calculated integral formula adopts the numerical integral formula; in forming the near-field unit Integral points are independently selected in ⁇ j and ⁇ i of finite elements, and any pair of integration points in ⁇ j and ⁇ i constitutes a peridynamic bond;
- (C0) Compile the calculation near-field unit according to the requirements of (B2) and (B3)
- the element stiffness matrix of The unit development subroutine of the unit is embedded into the commercial finite element software for use;
- step (C4) Call the calculated near-field unit written in step (C0)
- the element stiffness matrix of The unit development subroutine for the target structure ⁇ is used for damage analysis, and the commercial finite element software automatically completes the overall stiffness matrix integration;
- step (C5) Commercial finite element software is based on the load information input in step (C3) and the overall stiffness matrix obtained in step (C4) Automatically complete the solution of the overall node displacement vector d described in step (A5);
- step (C6) After determining the judging condition of whether the calculation result converges under the current given load, according to the overall node displacement vector d obtained in the step (C5), judge whether the calculation result under the current given load converges, if yes, Then enter step (C7); otherwise, then return to step (C4);
- step (C7) Judging whether the currently applied load value exceeds the preset maximum load value, if not, after increasing the load value, return to step (C3); otherwise, enter step (C8);
- step (C0) Write and calculate the element stiffness matrix of the 8-node near-field element
- the unit development subroutine is embedded into the commercial finite element software for use; among them, in the unit development subroutine, the update method of the function ⁇ ( ⁇ ,t) is as follows: after step (A5) is completed each time, according to the obtained The overall node displacement vector d, within the unit development subroutine, it is judged whether the peridynamic bond is broken, and the function ⁇ ( ⁇ ,t) is updated;
- step (C4) Call the calculated near-field unit written in step (C0)
- step (C5) Commercial finite element software can be based on the initial load information input in step (C3) and the overall stiffness matrix obtained in step (C4) Automatically complete the solution of the node displacement described in step (A5);
- step (C6) After determining the judging condition of whether the calculation result under the current given load is convergent, according to the overall node displacement vector d obtained in (C5), determine whether the displacement solution result under the current load is convergent, and if it is convergent, enter the step (C7); Otherwise, return to step (C4); Wherein, in the present embodiment, adopt whether newly-increased broken key is used as convergence criterion;
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)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (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之间的距离d
ji满足0≤d
ji≤δh，则将Ω
j和Ω
i组合成一个近场单元
(本发明的近场单元由且仅由两个有限单元组合而成，并且对组成该近场单元的两个有限单元的维数、形状、节点数等特征无特殊要求)，其中，j＝1,2,...,m，i＝1,2,...,m，m为有限单元的总个数，0≤δ≤100；
(A3)基于有限单元网格，计算总体载荷向量F，公式如下：
式中，m为有限单元的总个数；G
i为有限单元Ω
i的节点自由度的转换矩阵(d
i＝G
id，d
i为有限单元Ω
i的节点位移向量，d为总体节点位移向量)；·
T表示矩阵的转置；F
i为有限单元Ω
i的单元载荷向量；N
i(x)为组成
的有限单元Ω
i的形函数矩阵；x为组成
的有限单元Ω
i内的点；b(x)为外体力场在点x处的外体力向量；dV
x为积分微元；[·]表示矩阵；n
i为有限单元Ω
i的节点个数；N
α(x)为有限单元Ω
i的第α个节点的形函数，α＝1,2,...,n
i；
(A4)基于近场单元网格，计算总体刚度矩阵
公式如下：
式中，
为近场单元
的总个数；
为近场单元
的节点自由度的转换矩阵(
为近场单元
的节点位移向量，d为总体节点位移向量)；·
T表示矩阵的转置；
为近场单元
的单元刚度矩阵；
为近场单元
的形函数差矩阵；x′为组成
的有限单元Ω
j内的点；x为组成
的有限单元Ω
i内的点；D(ξ)为微模量矩阵；dV
x′为积分微元；dV
x为积分微元；[·]表示矩阵；N
j(x′)为组成
的有限单元Ω
j的形函数矩阵；N
i(x)为组成
的有限单元Ω
i的形函数矩阵；c
0(ξ)为键ξ的微模量系数，c
0(ξ)的值与目标结构Ω的材料性质有关；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数，其中0表示所述键ξ断裂，1表示所述键ξ未断裂；||·||表示计算向量的长度；ξ＝x′-x为近场动 力学键；ξ
1、ξ
2、ξ
3为键ξ向量的分量；
(A5)求解得到总体节点位移向量d，公式如下：
(A6)判断当前所给定的载荷下计算结果是否收敛，如果否，则返回步骤(A4)；反之，则进入步骤(A7)；
判断当前所给定的载荷下计算结果是否收敛的方法为：判断本次执行步骤(A5)得到的总体节点位移向量d
t与上一次执行步骤(A5)得到的总体节点位移向量d
t-1是否满足：||d
t-d
t-1||/||d
t||≤ε，其中ε为给定的误差限，10
-8h≤ε≤10
-1h(h为网格剖分尺寸)，如果是，则收敛；反之，则不收敛；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：判断有无新增断键，如果无新增断键，则收敛；反之，则不收敛；
判断键是否断裂的方法为：如果任一近场单元
中任一根键ξ的伸长率s大于某一给定的临界伸长率s
crit(取值范围为(-1,100])，则该键断裂，其中s
crit与目标结构Ω的材料性质有关；反之，则该键未断裂；s的计算公式如下：
u
i(x)＝N
i(x)d
i；
u
j(x′)＝N
j(x′)d
j；
式中，||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；x′为组成
的有限单元Ω
j内的点；x为组成
的有限单元Ω
i内的点；u
j(x′)为任一有限单元Ω
j中任一点x′处的位移向量；u
i(x)为任一有限单元Ω
i中任一点x处的位移向量；N
i(x)为组成
的有限单元Ω
i的形函数矩阵；d
i为有限单元Ω
i的节点位移向量；N
j(x′)为组成
的有限单元Ω
j的形函数矩阵；d
j为有限单元Ω
j的节点位移向量；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：先判断有无新增断键，如果无，则不更新函数μ(ξ,t)；反之，则更新函数μ(ξ,t)，再重新计算总体刚度矩阵
记更新后的总体刚度矩阵为
如果
(d
t为本次执行步骤(A5)得到的总体节点位移向量，F为总体载荷向量，ε为给定的误差限)，则收敛；反之，则不收敛；
(A7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后， 返回步骤(A3)；反之，则进入步骤(A8)；
(A8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
作为优选的技术方案：
如上所述的结构损伤分析的近场有限元法，步骤(A2)中，Ω
j和Ω
i之间的距离d
ji为Ω
j和Ω
i的形心距离，δ＝3；或者，Ω
j和Ω
i之间的距离d
ji为Ω
j和Ω
i的最接近点的距离，0≤δ≤100；或者，Ω
j和Ω
i之间的距离d
ji为Ω
j和Ω
i的最远离点的距离，0≤δ≤100。
如上所述的结构损伤分析的近场有限元法，步骤(A8)中，每一计算步t时目标结构Ω中任一点x处的等效损伤d
ξ(x,t)的计算公式如下：
式中，H
δ(x)为点x的近场邻域；
表示所有Ω
j的并集；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数；ω
crit为近场动力学键的临界断裂能；dV
x′为积分微元；x′为组成
的有限单元Ω
j内的点；x为组成
的有限单元Ω
i内的点；||·||表示计算向量的长度；0≤δ≤100；h为网格剖分尺寸；c
0(ξ)为所述键ξ的微模量系数，c
0(ξ)的值与目标结构Ω的材料性质有关；s
crit为给定的临界伸长率；ξ为近场动力学键。
本发明的上述结构损伤分析的近场有限元法与传统有限元法具有相同的计算流程，能够与商用有限元软件兼容。
本发明还提供如上任一项所述的结构损伤分析的近场有限元法在商用软件中的实现方法，利用商用有限元软件(例如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之间的距离d
ji满足0≤d
ji≤δh(0≤δ≤100，h为网格剖分尺寸)，则将Ω
j和Ω
i组合成一个近场单元
的节点的全局节点编号取为
其中
为近场单元
的节点个数。
如上所述的实现方法，具体步骤如下：
(C0)按(B2)和(B3)的要求编写计算近场单元
的单元刚度矩阵
的单元开发子程序，并嵌入到商用有限元软件中以供使用；
(C1)在计算机中对目标结构Ω进行几何建模，并设定网格剖分尺寸h的取值后，按网格剖分尺寸h生成传统的有限单元网格，得到m个有限单元；
(C2)确定距离d
ji的计算方式以及δ的取值后，按照步骤(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)为u
y＝2u
x＝0.35mm时的等效损伤云图；(b)为u
y＝2u
x＝0.5mm时的等效损伤云图；
图6为双开槽平板结构件X方向位移的计算结果：(a)为u
y＝2u
x＝0.35mm时的X方向位移云图；(b)为u
y＝2u
x＝0.5mm时的X方向位移云图；
图7为双开槽平板结构件Y方向位移的计算结果：(a)为u
y＝2u
x＝0.35mm时的Y方向位移云图；(b)为u
y＝2u
x＝0.5mm时的Y方向位移云图。
下面结合具体算例，给出具体实施方式，以进一步阐述本发明。应理解，该实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要 求书所限定的范围。
结构损伤分析的近场有限元法，具体步骤如下：
(A1)对目标结构Ω进行几何建模，并按网格剖分尺寸h生成传统的有限单元网格，得到m个有限单元，m的值由目标结构Ω的几何尺寸和采用的网格剖分尺寸h共同确定；
(A2)对任意两个有限单元Ω
j和Ω
i(如图1所示)，若Ω
j和Ω
i之间的距离d
ji满足0≤d
ji≤δh，则将Ω
j和Ω
i组合成一个近场单元
(如图1所示)，其中，j＝1,2,...,m，i＝1,2,...,m，m为有限单元的总个数，0≤δ≤100；其中，Ω
j和Ω
i之间的距离d
ji为Ω
j和Ω
i的形心距离，δ＝3；或者，Ω
j和Ω
i之间的距离d
ji为Ω
j和Ω
i的最接近点的距离，0≤δ≤100；或者，Ω
j和Ω
i之间的距离d
ji为Ω
j和Ω
i的最远离点的距离，0≤δ≤100；
(A3)基于有限单元网格，计算总体载荷向量F，公式如下：
式中，m为有限单元的总个数；G
i为有限单元Ω
i的节点自由度的转换矩阵(d
i＝G
id，d
i为有限单元Ω
i的节点位移向量，d为总体节点位移向量)；·
T表示矩阵的转置；F
i为有限单元Ω
i的单元载荷向量；d
i为有限单元Ω
i的节点位移向量；d为总体节点位移向量；N
i(x)为组成
的有限单元Ω
i的形函数矩阵；x为组成
的有限单元Ω
i内的点；b(x)为外体力场在点x处的外体力向量；dV
x为积分微元；[·]表示矩阵；n
i为有限单元Ω
i的节点个数；N
α(x)为有限单元Ω
i的第α个节点的形函数，α＝1,2,...,n
i；
(A4)基于近场单元网格，计算总体刚度矩阵
公式如下：
式中，
为近场单元
的总个数；
为近场单元
的节点自由度的转换矩阵(
为近场单元
的节点位移向量，d为总体节点位移向量)；·
T表示矩阵的转置；
为近场单元
的单元刚度矩阵；
为近场单元
的节点位移向量；d为总体节点位移向量；
为近场单元
的形函数差矩阵；x′为组成
的有限单元Ω
j内的点；x为组成
的有限单元Ω
i内的点；D(ξ)为微模量矩阵；dV
x′为积分微元；dV
x为积分微元；[·]表示矩阵；N
j(x′)为组成
的有限单元Ω
j的形函数矩阵；N
i(x)为组成
的有限单元Ω
i的形函数矩阵；c
0(ξ)为所述键ξ的微模量系数，c
0(ξ)的值与目标结构Ω的材料性质有关；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数，其中0表示所述键ξ断裂，1表示所述键ξ未断裂；||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；ξ
1、ξ
2、ξ
3为键ξ向量的分量；
(A5)求解得到总体节点位移向量d，公式如下：
(A6)判断当前所给定的载荷下计算结果是否收敛，如果否，则返回步骤(A4)；反之，则进入步骤(A7)；
判断当前所给定的载荷下计算结果是否收敛的方法为：判断本次执行步骤(A5)得到的总体节点位移向量d
t与上一次执行步骤(A5)得到的总体节点位移向量d
t-1是否满足：||d
t-d
t-1||/||d
t||≤ε，其中ε为给定的误差限，10
-8h≤ε≤10
-1h(h为网格剖分尺寸)，如果是，则收敛；反之，则不收敛；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：判断有无新增断键，如果无新增断键，则收敛；反之，则不收敛；
判断键是否断裂的方法为：如果任一近场单元
中任一根键ξ的伸长率s大于某一给定的临界伸长率s
crit(取值范围为(-1,100])，则该键断裂，其中s
crit与目标结构Ω的材料性质有关；反之，则该键未断裂；s的计算公式如下：
u
i(x)＝N
i(x)d
i；
u
j(x′)＝N
j(x′)d
j；
式中，||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；x′为组成
的有限单元Ω
j内的点；x为组成
的有限单元Ω
i内的点；u
j(x′)为任一有限单元Ω
j中任一点x′处的位移向量；u
i(x)为任一有限单元Ω
i中任一点x处的位移向量；N
i(x)为组成
的有限单元Ω
i的形函数矩阵；d
i为有限单元Ω
i的节点位移向量；N
j(x′)为组成
的有限单元Ω
j的形函数矩阵；d
j为有限单元Ω
j的节点位移向量；
或者，判断当前所给定的载荷下计算结果是否收敛的方法为：先判断有无新增断键，如果无，则不更新函数μ(ξ,t)；反之，则更新函数μ(ξ,t)，再重新计算总体刚度矩阵
记更新后的总体刚度矩阵为
如果
(d
t为本次执行步骤(A5)得到的总体节点位移向量，F为总体载荷向量，ε为给定的误差限)，则收敛；反之，则不收敛；
(A7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(A3)；反之，则进入步骤(A8)；
(A8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图；其中，每一计算步t时目标结构Ω中任一点x处的等效损伤d
ξ(x,t)的计算公式如下：
式中，H
δ(x)为点x的近场邻域；
表示所有Ω
j的并集；μ(ξ,t)为相关于所述键ξ和计算步t的取值为0或1的函数；ω
crit为近场动力学键的临界断裂能；dV
x′为积分微元；x′为组成
的有限单元Ω
j内的点；x为组成
的有限单元Ω
i内的点；||·||表示计算向量的长度；0≤δ≤100；h为网格剖分尺寸；c
0(ξ)为所述键ξ的微模量系数，c
0(ξ)的值与目标结构Ω的材料性质有关；s
crit为给定的临界伸长率；ξ为近场动力学键。
结构损伤分析的近场有限元法在商用软件中的实现方法，利用商用有限元软件(例如ANSYS、ABAQUS、MSC Nastran、MSC Marc、ADINA等)的二次开发功能，在商用有限元软件中执行结构损伤分析的近场有限元法，包括以下操作：
(B1)由有限单元网格数据生成近场单元网格数据，并输入到商用有限元软件中；具体地，设有限单元Ω
j的节点的全局节点编号为
设有限单元Ω
i的节点的全局节点编号为
若Ω
j和Ω
i之间的距离d
ji满足0≤d
ji≤δ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)确定距离d
ji的计算方式以及δ的取值后，按照步骤(A2)和(B1)中所述的方式生成近场单元网格，得到
个近场单元；
(C3)在商用有限元软件中，输入目标结构Ω所受到的载荷信息；
(C4)调用步骤(C0)中编写的计算近场单元
的单元刚度矩阵
的单元开发子程序对目标结构Ω进行损伤分析，商用有限元软件自动完成总体刚度矩阵
的集成；
(C5)商用有限元软件根据步骤(C3)中输入的载荷信息和步骤(C4)中得到的总体刚度矩阵
自动完成步骤(A5)中所述的总体节点位移向量d的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据步骤(C5)中得到的总体节点位移向量d，判断当前所给定的载荷下计算结果是否收敛，如果是，则进入步骤(C7)；反之，则返回步骤(C4)；
(C7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(C3)；反之，则进入步骤(C8)；
(C8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
现结合具体使用过程进行举例说明：
问题设置：考虑一个二维问题，具体是一个具有两条水平开槽的平板(双开槽平板)，其几何尺寸和预设的最大载荷条件如图3所示，预设的最大载荷条件为水平方向位移u
x＝0.25mm，竖直方向u
y＝0.5mm；本使用过程不考虑外力载荷的作用，目标结构的材料的杨氏模量和泊松比分别被设置为E＝30GPa和ν＝1/3；二维时，基于有限单元网格生成近场单元的示意图如图2所示；数值计算时，采用非结构化四节点四边形有限单元网格，具体如图4所示；
结构损伤分析的近场有限元法在商用软件(ANSYS)中的实现方法，具体步骤如下：
(C0)编写计算8节点近场单元的单元刚度矩阵
的单元开发子程序，并嵌入到商用有限元软件中以供使用；其中，单元开发子程序中，函数μ(ξ,t)的更新方式为：每次完成步骤(A5)后，根据得到的总体节点位移向量d，在单元开发子程序内部判断近场动力学键是否断裂，并更新函数μ(ξ,t)；
(C1)在计算机中对目标结构进行几何建模，并按网格剖分尺寸h＝1.2生成有限单元网格，得到m＝27741个有限单元；
(C2)取d
ji为Ω
j和Ω
i的形心距离，并取δ＝3，按照步骤(A2)和步骤(B1)中所述的方式生成近场单元网格，得到
个近场单元；
(C3)在商用有限元软件中，输入目标结构的初始载荷信息和预设的最大载荷信息；其中初始载荷信息为：水平方向位移u
x＝0.025mm，竖直方向u
y＝0.05mm；
(C4)调用步骤(C0)中编写的计算近场单元
的单元刚度矩阵
的单元开发子程序对目标结构进行损伤分析；商用有限元软件能够自动完成总体刚度矩阵
的集成；
(C5)商用有限元软件能够根据步骤(C3)中输入的初始载荷信息和和步骤(C4)中得到的总体刚度矩阵
自动完成步骤(A5)中所述的节点位移的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据(C5)中得到的总体节点位移向量d，判定当前载荷下的位移求解结果是否收敛，如果收敛，则进入步骤(C7)；否则，返回步骤(C4)；其中，本实施例中，采用是否有新增断键作为收敛准则；
(C7)判断当前所施加的载荷值是否超过预设最大载荷值；如果没有超过，返回步骤(C3)，在初始载荷信息的基础上增加载荷值(水平方向位移u
x＝0.025mm，竖直方向u
y＝0.05mm)得到新的初始载荷信息；反之，则进入步骤(C8)；
(C8)根据总体节点位移向量d和断键信息，输出目标结构在每一计算步的位移云图和等效损伤云图；当u
y＝2u
x＝0.35mm时，等效损伤云图如图5(a)所示，X方向位移云图如图6(a)所示，Y方向位移云图如图7(a)所示；当u
y＝2u
x＝0.5mm时，等效损伤云图如图5(b)所示，X方向位移云图如图6(b)所示，Y方向位移云图如图7(b)所示。
Claims (6)
- 结构损伤分析的近场有限元法，其特征在于，包括以下步骤：(A1)对目标结构Ω进行几何建模，并按网格剖分尺寸h生成有限单元网格，得到m个有限单元；(A2)对任意两个有限单元Ω j和Ω i，若Ω j和Ω i之间的距离d ji满足0≤d ji≤δh，则将Ω j和Ω i组合成一个近场单元其中，j＝1,2,…,m，i＝1,2,…,m，0≤δ≤100；
(A3)基于有限单元网格，计算总体载荷向量F，公式如下：
式中，G i为有限单元Ω i的节点自由度的转换矩阵；F i为有限单元Ω i的单元载荷向量；N i(x)为组成的有限单元Ω i的形函数矩阵；x为组成
的有限单元Ω i内的点；b(x)为外体力场在点x处的外体力向量；n i为有限单元Ω i的节点个数；N α(x)为有限单元Ω i的第α个节点的形函数，α＝1,2,…,n i；
(A4)基于近场单元网格，计算总体刚度矩阵公式如下：
式中，为近场单元
的总个数；
为近场单元
的节点自由度的转换矩阵；
为近场单元
的单元刚度矩阵；
为近场单元
的形函数差矩阵；x′为组成
的有限单元Ω j内的点；D(ξ)为微模量矩阵；N j(x′)为组成
的有限单元Ω j的形函数矩阵；c 0(ξ)为键ξ的微模量系数，μ(ξ,t)为相关于所述键ξ和计算步t的取值为 0或1的函数，其中0表示所述键ξ断裂，1表示所述键ξ未断裂；||·||表示计算向量的长度；ξ＝x′-x为近场动力学键；ξ 1、ξ 2、ξ 3为键ξ向量的分量；
(A5)求解得到总体节点位移向量d，公式如下：
(A6)判断当前所给定的载荷下计算结果是否收敛，如果否，则返回步骤(A4)；反之，则进入步骤(A7)；判断当前所给定的载荷下计算结果是否收敛的方法为：判断本次执行步骤(A5)得到的总体节点位移向量d t与上一次执行步骤(A5)得到的总体节点位移向量d t-1是否满足：||d t-d t-1||/||d t||≤ε，10 -8h≤ε≤10 -1h，如果是，则收敛；反之，则不收敛；或者，判断当前所给定的载荷下计算结果是否收敛的方法为：判断有无新增断键，如果无新增断键，则收敛；反之，则不收敛；判断键是否断裂的方法为：如果任一近场单元中任一根键ξ的伸长率s大于某一给定的临界伸长率s crit，则该键断裂，其中s crit与目标结构Ω的材料性质有关；反之，则该键未断裂；s的计算公式如下：
u i(x)＝N i(x)d i；u j(x′)＝N j(x′)d j；式中，u j(x′)为任一有限单元Ω j中任一点x′处的位移向量；u i(x)为任一有限单元Ω i中任一点x处的位移向量；d i为有限单元Ω i的节点位移向量；d j为有限单元Ω j的节点位移向量；或者，判断当前所给定的载荷下计算结果是否收敛的方法为：先判断有无新增断键，如果无，则不更新函数μ(ξ,t)；反之，则更新函数μ(ξ,t)，再重新计算总体刚度矩阵记更新后的总体刚度矩阵为
如果
则收敛；反之，则不收敛；
(A7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(A3)；反之，则进入步骤(A8)；(A8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。 - 根据权利要求1所述的结构损伤分析的近场有限元法，其特征在于，步骤(A2)中，Ω j和Ω i之间的距离d ji为Ω j和Ω i的形心距离，δ＝3；或者，Ω j和Ω i之间的距离d ji为Ω j和Ω i的最接近点的距离，0≤δ≤100；或者，Ω j和Ω i之间的距离d ji为Ω j和Ω i的最远离点的距离，0≤δ≤100。
- 根据权利要求1所述的结构损伤分析的近场有限元法，其特征在于，步骤(A8)中，每一计算步t时目标结构Ω中任一点x处的等效损伤d ξ(x,t)的计算公式如下：
式中，H δ(x)为点x的近场邻域；表示所有Ω j的并集；ω crit为近场动力学键的临界断裂能。
- 如权利要求1～3任一项所述的结构损伤分析的近场有限元法在商用软件中的实现方法，其特征在于，利用商用有限元软件的二次开发功能，在商用有限元软件中执行结构损伤分析的近场有限元法，包括以下操作：(B1)由有限单元网格数据生成近场单元网格数据，并输入到商用有限元软件中；(B2)将步骤(A4)中近场单元的单元刚度矩阵
的计算公式编写为商用有限元软件的单元开发子程序，并嵌入到该商用有限元软件中；子程序中关于近场单元
的单元刚度矩阵
计算的积分公式采用数值积分公式；在组成近场单元
的有限单元Ω j和Ω i内各自独立地选取积分点，并且Ω j和Ω i中的任一对积分点构成一根近场动力学键；
(B3)设定步骤(A6)中函数μ(ξ,t)的更新方式为：每次完成步骤(A5)后，根据得到的总体节点位移向量d，在单元开发子程序内部判断近场动力学键是否断裂，并更新函数μ(ξ,t)。 - 根据权利要求4所述的实现方法，其特征在于，(B1)中，设有限单元Ω j的节点的全局 节点编号为设有限单元Ω i的节点的全局节点编号为
若Ω j和Ω i之间的距离d ji满足0≤d ji≤δh，则将Ω j和Ω i组合成一个近场单元
的节点的全局节点编号取为
其中
为近场单元
的节点个数。
- 根据权利要求4或5所述的实现方法，其特征在于，具体步骤如下：(C0)按(B2)和(B3)的要求编写计算近场单元的单元刚度矩阵
的单元开发子程序，并嵌入到商用有限元软件中以供使用；
(C1)在计算机中对目标结构Ω进行几何建模，并设定网格剖分尺寸h的取值后，按网格剖分尺寸h生成有限单元网格，得到m个有限单元；(C2)确定距离d ji的计算方式以及δ的取值后，按照步骤(A2)和(B1)中所述的方式生成近场单元网格，得到个近场单元；
(C3)在商用有限元软件中，输入目标结构Ω所受到的载荷信息；(C4)调用步骤(C0)中编写的计算近场单元的单元刚度矩阵
的单元开发子程序对目标结构Ω进行损伤分析，商用有限元软件自动完成总体刚度矩阵
的集成；
(C5)商用有限元软件根据步骤(C3)中输入的载荷信息和步骤(C4)中得到的总体刚度矩阵自动完成步骤(A5)中所述的总体节点位移向量d的求解；
(C6)确定当前所给定的载荷下计算结果是否收敛的判断条件后，根据步骤(C5)中得到的总体节点位移向量d，判断当前所给定的载荷下计算结果是否收敛，如果是，则进入步骤(C7)；反之，则返回步骤(C4)；(C7)判断当前所施加的载荷值是否超过预设最大载荷值，如果否，则增加载荷值后，返回步骤(C3)；反之，则进入步骤(C8)；(C8)基于总体节点位移向量d，输出目标结构Ω在每一计算步的位移云图，同时基于断键信息，输出目标结构Ω在每一计算步的等效损伤云图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| 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 |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110885638.8A CN113705040B (zh) | 2021-08-03 | 2021-08-03 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
| CN202110885638.8 | 2021-08-03 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| WO2023011029A1 true WO2023011029A1 (zh) | 2023-02-09 |
Family
ID=78651448
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2022/100311 WO2023011029A1 (zh) | 2021-08-03 | 2022-06-22 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US20230325560A1 (zh) |
| CN (1) | CN113705040B (zh) |
| WO (1) | WO2023011029A1 (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116050027A (zh) * | 2023-03-30 | 2023-05-02 | 陕西空天信息技术有限公司 | 叶轮叶片结构静态分析方法、计算机程序产品和电子设备 |
| CN116364218A (zh) * | 2023-03-30 | 2023-06-30 | 沈阳工业大学 | 基于近场动力学岩石材料率效应本构模型构建与实施方法 |
Families Citing this family (4)
| 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软件的批量修改单元材料方法 |
| CN117688807A (zh) * | 2023-12-01 | 2024-03-12 | 河海大学 | 基于间接力键有限元法的结构变形与损伤分析方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110457790A (zh) * | 2019-07-26 | 2019-11-15 | 顾鑫 | 用于结构变形分析的近场动力学非连续伽辽金有限元方法 |
| CN111339594A (zh) * | 2020-02-26 | 2020-06-26 | 河北工业大学 | 基于dic技术的近场动力学参数实验反演系统及使用方法 |
| CN111814310A (zh) * | 2020-06-11 | 2020-10-23 | 大连理工大学 | 强度准则驱动的近场动力学模型预测结构破坏的方法 |
| CN113705040A (zh) * | 2021-08-03 | 2021-11-26 | 大连理工大学 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
Family Cites Families (4)
| 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 |
| CN111159951B (zh) * | 2019-12-31 | 2023-08-18 | 哈尔滨工业大学（深圳） | 一种基于abaqus有限元与边界元的耦合方法 |
| CN112446163B (zh) * | 2020-11-24 | 2022-12-09 | 西安交通大学 | 基于参数化水平集的能量有限元拓扑优化方法 |
| CN112966421B (zh) * | 2021-03-16 | 2023-12-26 | 昆明理工大学 | 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 |
-
2021
- 2021-08-03 CN CN202110885638.8A patent/CN113705040B/zh active Active
-
2022
- 2022-06-22 US US18/024,780 patent/US20230325560A1/en not_active Abandoned
- 2022-06-22 WO PCT/CN2022/100311 patent/WO2023011029A1/zh unknown
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110457790A (zh) * | 2019-07-26 | 2019-11-15 | 顾鑫 | 用于结构变形分析的近场动力学非连续伽辽金有限元方法 |
| CN111339594A (zh) * | 2020-02-26 | 2020-06-26 | 河北工业大学 | 基于dic技术的近场动力学参数实验反演系统及使用方法 |
| CN111814310A (zh) * | 2020-06-11 | 2020-10-23 | 大连理工大学 | 强度准则驱动的近场动力学模型预测结构破坏的方法 |
| CN113705040A (zh) * | 2021-08-03 | 2021-11-26 | 大连理工大学 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
Non-Patent Citations (2)
| Title |
|---|
| HAN FEI; LUBINEAU GILLES; AZDOUD YAN; ASKARI ABE: "A morphing approach to couple state-based peridynamics with classical continuum mechanics", COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING, NORTH-HOLLAND, AMSTERDAM, NL, vol. 301, 4 January 2016 (2016-01-04), AMSTERDAM, NL , pages 336 - 358, XP029411294, ISSN: 0045-7825, DOI: 10.1016/j.cma.2015.12.024 * |
| SHEN SHANGKUN, YANG ZIHAO, HAN FEI, CUI JUNZHI, ZHANG JIEQIONG: "Peridynamic modeling with energy-based surface correction for fracture simulation of random porous materials", THEORETICAL AND APPLIED FRACTURE MECHANICS, ELSEVIER BV, NL, vol. 114, 1 August 2021 (2021-08-01), NL , pages 102987, XP093031372, ISSN: 0167-8442, DOI: 10.1016/j.tafmec.2021.102987 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116050027A (zh) * | 2023-03-30 | 2023-05-02 | 陕西空天信息技术有限公司 | 叶轮叶片结构静态分析方法、计算机程序产品和电子设备 |
| CN116364218A (zh) * | 2023-03-30 | 2023-06-30 | 沈阳工业大学 | 基于近场动力学岩石材料率效应本构模型构建与实施方法 |
| CN116050027B (zh) * | 2023-03-30 | 2023-07-07 | 陕西空天信息技术有限公司 | 叶轮叶片结构静态分析方法、存储介质和电子设备 |
| CN116364218B (zh) * | 2023-03-30 | 2023-09-08 | 沈阳工业大学 | 基于近场动力学岩石材料率效应本构模型构建与实施方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| US20230325560A1 (en) | 2023-10-12 |
| CN113705040A (zh) | 2021-11-26 |
| CN113705040B (zh) | 2024-03-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2023011029A1 (zh) | 结构损伤分析的近场有限元法及在商用软件中的实现方法 | |
| Zhang et al. | Structural topology optimization through explicit boundary evolution | |
| Tang | Moving mesh methods for computational fluid dynamics | |
| Benson et al. | A generalized finite element formulation for arbitrary basis functions: from isogeometric analysis to XFEM | |
| Rannou et al. | A local multigrid X‐FEM strategy for 3‐D crack propagation | |
| Tian et al. | An efficient hybrid method for multibody dynamics simulation based on absolute nodal coordinate formulation | |
| WO2012026383A1 (ja) | 計算用データ生成装置、計算用データ生成方法及び計算用データ生成プログラム | |
| Cui et al. | An ABAQUS implementation of the cell-based smoothed finite element method (CS-FEM) | |
| Lv et al. | A matrix-free implicit unstructured multigrid finite volume method for simulating structural dynamics and fluid–structure interaction | |
| Wang et al. | From computer-aided design (CAD) toward human-aided design (HAD): an isogeometric topology optimization approach | |
| Liang et al. | Extended material point method for the three‐dimensional crack problems | |
| Li et al. | Phase field crack model with diffuse description for fracture problem and implementation in engineering applications | |
| Yang et al. | Isogeometric double-objective shape optimization of free-form surface structures with Kirchhoff–Love shell theory | |
| Shao | Comparison of BESO and SIMP to do structural topology optimization in discrete digital design, and then combine them into a hybrid method | |
| Rong et al. | Hybrid finite element transfer matrix method and its parallel solution for fast calculation of large-scale structural eigenproblem | |
| Yang et al. | Morphogenesis of free-form surfaces by an effective approach based on isogeometric analysis and particle swarm optimization | |
| CN107688691A (zh) | 一种木结构用齿板连接节点受力性能的数值模拟方法 | |
| Sotiropoulos et al. | High performance topology optimization computing platform | |
| Su et al. | Computational morphogenesis of free‐form grid structures with Voronoi diagram | |
| Fan et al. | Cone-complementary manifold method for stability and failure analysis of jointed/fractured rock masses | |
| Wang et al. | Extended Eulerian SPH and its realization of FVM | |
| Peng et al. | Quadtree-polygonal smoothed finite element method for adaptive brittle fracture problems | |
| Xu et al. | An integrated method of CAD, CAE and multi-objective optimization | |
| Gao et al. | Overview of Advanced Numerical Methods Classified by Operation Dimensions | |
| Gou et al. | Automatic modelling of heterotypic tunnel structures via NURBS-based meshfree method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| 121 | Ep: the epo has been informed by wipo that ep was designated in this application |
Ref document number: 22851745 Country of ref document: EP Kind code of ref document: A1 |
|
| NENP | Non-entry into the national phase |
Ref country code: DE |