CN113722965B - 一种基于积分-广义有限差分数值离散算子的断裂模拟方法 - Google Patents

一种基于积分-广义有限差分数值离散算子的断裂模拟方法 Download PDF

Info

Publication number
CN113722965B
CN113722965B CN202111044179.7A CN202111044179A CN113722965B CN 113722965 B CN113722965 B CN 113722965B CN 202111044179 A CN202111044179 A CN 202111044179A CN 113722965 B CN113722965 B CN 113722965B
Authority
CN
China
Prior art keywords
discrete
point
integral
discrete point
points
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
CN202111044179.7A
Other languages
English (en)
Other versions
CN113722965A (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.)
Wuhan Institute of Rock and Soil Mechanics of CAS
Original Assignee
Wuhan Institute of Rock and Soil Mechanics of CAS
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 Wuhan Institute of Rock and Soil Mechanics of CAS filed Critical Wuhan Institute of Rock and Soil Mechanics of CAS
Priority to CN202111044179.7A priority Critical patent/CN113722965B/zh
Publication of CN113722965A publication Critical patent/CN113722965A/zh
Application granted granted Critical
Publication of CN113722965B publication Critical patent/CN113722965B/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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential 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
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • 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)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Business, Economics & Management (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Strategic Management (AREA)
  • Quality & Reliability (AREA)
  • Game Theory and Decision Science (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Development Economics (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明属于数值模拟方法、裂缝模拟技术领域，具体涉及一种基于积分‑广义有限差分数值离散算子的断裂模拟方法；对研究域进行空间离散，利用数值离散算子构建各点场量值及其导数的近似表达，后代入等效积分弱形式平衡方程，结合边界条件，组装并求解以离散点位移场为未知量的整体方程组，通过节点生死模拟裂纹扩展行为，并重复以上步骤，直至所得断裂参数不再满足裂纹扩展准则；本发明计算原理简单、易实现，通过节点生死可方便地模拟裂纹向任意方向、任意长度的扩展问题，应用于不同扩展准则下各类型、各维度的裂纹扩展研究，且无需引入断裂力学先验解，方程物理意义明确；本发明实现了有无网格的二象性，可为裂模模拟提供新发展思路。

Description

一种基于积分-广义有限差分数值离散算子的断裂模拟方法
技术领域
本发明属于数值模拟方法、裂缝模拟技术领域，具体涉及一种基于积分-广义有限差分数值离散算子的断裂模拟方法。
背景技术
在工程技术领域，许多工程材料、构件及系统的破坏与裂纹衍生、扩展、贯通是密切相关的，理论分析、试验预测、数值模拟是开展裂纹方面工作研究的主要手段，其中，数值模拟对各种问题适用性强、应用范围广，这些特点使得它在裂纹及其扩展问题中的研究被充分重视。数值方法是一种计算机算法，其基础是数值离散化格式，裂纹模拟是建立在离散格式基础上的，良好的离散格式有利于实现断裂的模拟。然而，断裂模拟涉及到复杂的非线性以及非连续数值模拟问题，一直是研究的难点和热点。
人们已经发展了众多的数值离散及裂纹模拟方法，常用的裂纹数值模拟方法有：扩展有限单元法、无网格伽辽金方法、单位分解有限元法以及流形元法等。其中，扩展有限元法的核心思想是用扩充带有不连续性质的形函数代表区域内的间断，不连续场描述独立于网格边界，而可以方便地模拟裂纹任意方向、任意长度的扩展问题，还可模拟带有孔洞和夹杂的非均质材料，但仍是在标准有限元的框架内解决问题；无网格伽辽金方法使用高斯积分方法对求解域中的背景方法进行积分，采用移动最小二乘技术来构造形函数，但同时也丧失了过点插值的性质，而不便于施加本质边界条件。由于裂纹问题的高度非线性和复杂性，尚没有一种类型数值方法能够有效解决其扩展的各类问题，主要存在以下不足：一是连续-非连续刻画的困难，二是对网格的依赖性，三是边界条件的准确施加。有网格和无网格的方法各有优势和劣势，尽管将二者机械结合的方法也已存在，但如何自然地将二者融合、发挥二者优势和作用尚未见到实质的突破。连续-非连续问题的重要性和复杂性，呼唤相关研究人员丰富、发展新的数值分析思路，提出新的解决方案。
发明内容
针对现有方法存在的问题，本发明从底层离散格式着手，提供了一种基于积分-广义有限差分数值离散算子的断裂模拟方法，本发明的基本思路是：联合泰勒级数展开技术和移动最小二乘技术，建立研究域内各离散点处场变量及其导数值的离散近似格式，又称之为离散算子，并将其代入平衡方程的等效积分弱形式中，结合边界条件，求解研究域离散点的位移场量值，进一步得到断裂力学参数，根据裂纹扩展准则执行节点生死算法，模拟裂纹扩展行为。本发明所述数值离散算子克服了对网格的依赖性，具有计算原理简单、易于实现的特点，高计算效率、高精度，灵活性与应用性强，同时，因平衡方程采用等效积分弱形式，降低了对场函数本身连续性的要求；本发明所述断裂模拟方法，是所述积分-广义有限差分数值离散算子在裂纹问题中的扩展应用，可方便地引入边界条件，且针对裂纹问题，无需先验引入基于断裂力学解的增强函数。本发明可随时切换有网格和无网格，实现有无网格的二象性，这可能为进一步解决现有数值方法在模拟裂纹方面的不足提供新方向。
为达到上述目的，本发明所采用的技术方案如下：
一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于，使用积分-广义有限差分方法离散并求解线弹性力学问题，通过断裂准则对裂纹行为进行预测，并通过节点生死实现裂纹及其扩展过程的追踪。
一种基于积分-广义有限差分数值离散算子的断裂模拟方法，具体包括如下步骤：
步骤一：执行空间离散算法，将研究域空间离散化，获得研究域空间内一系列离散点的编号及其坐标；
步骤二：利用步骤一离散点对边界线条进行逆时针描述，形成闭合回路，限定各边界线条的边界条件；
步骤三：确定平衡方程的等效积分弱形式；然后对步骤一所得的每个离散点均按照步骤四-步骤八操作；
步骤四：对离散点执行积分域搜索算法，确定并记录每个离散点的积分域信息；
步骤五：根据步骤一离散点坐标信息与步骤二逆时针边界描述，执行数值离散算法将离散点的场变量值及其导数近似表达为关于局部离散点场变量的关系式；
步骤六：赋予步骤五场变量为位移的物理含义，根据步骤二边界条件，执行自然边界条件引入算法，得到引入自然边界条件的平衡方程等效积分弱形式；
步骤七：将步骤五离散点处导数关于局部离散点位移场变量的关系式代入步骤六引入自然边界条件的平衡方程等效积分弱形式，扩充形成关于所有离散点处位移场变量值的线性方程；
步骤八：根据步骤二边界条件，执行本质边界条件引入算法，得到引入本质边界条件的线性方程；
步骤九：组装步骤八引入本质边界条件的线性方程得到线性方程组；
步骤十：执行计算，得到所有离散点的位移场变量值；
步骤十一：断裂力学参数计算；
步骤十二：裂纹扩展准则判断：根据步骤十一断裂力学参数，
若离散点的位移满足裂纹扩展准则，执行步骤十三，并将满足扩展准则的离散点重复步骤三至步骤十二，每一次重复为一次计算步；
若离散点的位移不满足裂纹扩展准则，说明该状态下该离散点所在的研究区域已达到变形稳定阶段，裂纹不再扩展，跳至步骤十四；
步骤十三：节点生死算法，离散点集合和边界条件的更新，裂纹扩展行为预测；
步骤十四：计算结果的可视化与后处理：将每一计算步下的步骤十离散点的位移场变量值、步骤十一离散点应力场变量值、应变场变量值、断裂力学参数，及直至达到变形稳定阶段时的裂纹扩展路径可视化或进行其他后处理，以满足人员特别需要。
所述步骤一是通过调用外源软件或采用内置算法进行的；
其中通过调用外源软件具体包括如下步骤：
步骤1，向外源软件内输入研究域的边界信息；
步骤2，利用外源软件完成研究域的网格剖分；
步骤3，输出网格点的编号及坐标信息，得到的网格点编号及坐标信息即为研究域内离散点的编号及其坐标(xi,yi),i＝1,2,...,N；
其中采用内置算法具体包括如下步骤：
步骤1，根据研究域的边界信息，确定边界的下限、上限、左限、右限，即ymin、ymax、xmin、xmax，围成最大整合区域；
步骤2，确定最大整合区域内x、y轴等分数量Nx、Ny，并将最大整合区域等分，得到
Figure GDA0004246385180000032
个交点；
步骤3，剔除位于研究域之外的交点，剩余交点即为研究域内的离散点，对其进行编号并记录坐标信息(xi,yi),i＝1,2,...,N。
所述步骤三确定的平衡方程的等效积分弱形式，针对忽略体积力的线弹性力学问题中，具体如下：
Figure GDA0004246385180000031
式中，Ln是围绕微元面积A的逆时针闭合回路中线条的数量，LSj是第j(j＝1,2,...,Ln)条边界线条的长度，sp,j、ep,j分别代表第j条边界线条的起始点和终点,(nx,j,ny,j)是第j条边界线条的外矢量方向，
Figure GDA0004246385180000041
Figure GDA0004246385180000042
式中，针对平面应变问题有
Figure GDA0004246385180000043
对平面应变问题E'＝E，ν'＝ν，E为弹性模量，v为泊松比，w1和w2分别为x、y方向上的位移。
所述步骤四是通过单元覆盖方法或随机覆盖方法进行的；
其中单元覆盖方法为根据步骤一离散点的编号和坐标信息，生成Delaunay三角形单元，包含离散点的Delaunay三角形单元即为该离散点的积分域；
其中随机覆盖方法具体包括以下步骤：
步骤1，根据研究域类型确定积分域的类型，当研究域为长条形时，选择长方形或椭圆形为积分域类型，当研究域为圆形时，选择圆形或正方形作为积分域类型，再确定积分域的形状参数；
步骤2，根据步骤一离散点的编号和坐标信息，筛选出距离离散点小于或者等于步骤1积分域形状参数的积分域内的局部离散点；
步骤3，根据步骤2积分域内的局部离散点，其中属于步骤一研究域内的离散点的积分域作为该离散点的积分域。
所述步骤五中目标离散点i的场量值及其导数近似表达为关于局部离散点场变量的关系式如下：
G＝DT·ST·W·S
H＝DT·ST·W
Figure GDA0004246385180000044
其中：
Figure GDA0004246385180000045
|△x|max＝max{|△x1|,|△x2|,...,|△xn|}，|△y|max＝max{|△y1|,|△y2|,...,|△yn|}，T代表转置，
Figure GDA0004246385180000051
W＝diag([W1 W2 … Wn])，
Figure GDA0004246385180000052
目标离散点i的周围有n个局部离散点，(xj,yj)是第j个局部离散点的坐标，j＝1、2、……n，△xj＝xj-x0，△yj＝yj-y0，u0、
Figure GDA0004246385180000053
为目标离散点i(x0,y0)处的场变量和各阶导数。
所述步骤六的具体过程如下：
步骤6.1，判断目标离散点i的积分域是否包含步骤二中边界线条，若是，执行步骤6.2；若否，跳至步骤七；
步骤6.2，进一步判断步骤6.1中包含的边界线条是否为自然边界条件，若是，执行步骤6.3；若否，跳至步骤七；
步骤6.3，将步骤三得到平衡方程的等效积分弱形式中关于该边界线条的表达移动到等号右侧，得到引入自然边界条件的平衡方程等效积分弱形式。
所述步骤七将步骤五目标离散点i处导数关于局部离散点位移场变量的关系式代入步骤6.3引入自然边界条件的平衡方程的等效积分弱形式，若是由步骤6.1或步骤6.2直接调至本步骤的则将步骤五目标离散点i处导数关于局部离散点位移场变量的关系式代入步骤三的平衡方程的等效积分弱形式，扩充形成关于所有离散点处位移场变量值的线性方程；
其中扩充形成关于所有离散点处位移场变量值的线性方程为：
Figure GDA0004246385180000054
其中：Pi，Qi为经扩充之后大小为1×2N的系数矩阵，Ri＝[Ri,P Ri,Q]T为扩充后大小为2×1的右边项矩阵，U＝[u1 1 u1 2 u2 1 u2 2 … … uN 1 uN 2]T为扩充后大小为2N×1的未知量矩阵。
所述步骤八的具体过程包括如下步骤：
步骤8.1，对于目标离散点i，判断其是否包含步骤二中边界线条，若是，执行步骤8.2，若否，跳至步骤9；
步骤8.2，进一步判断步骤8.1所包含边界线条是否为本质边界条件，若是，执行步骤8.3；若否，跳至步骤9；
步骤8.3，若在目标离散点i处施加数值为DB(DB≠0)、方向为x方向的本质边界条件,则令步骤七所得线性方程中的Pi,i＝||Pi T||1×1025，Ri,P＝DB×||Pi T||1×1025，|| ||1为一阶范数；若在目标离散点i处施加数值为DB(DB≠0)、方向为y方向的本质边界条件,则令步骤七所得线性方程中的Qi,i＝||Qi T||1×1025，Ri,Q＝DB×||Ri T||1×1025；若在目标离散点i处施加数值为0、方向为x方向的本质边界条件,则令步骤七所得线性方程中的Pi,i＝||Pi T||1×1025；若在目标离散点i处施加数值为0、方向为y方向的本质边界条件,则令步骤七所得线性方程中的Qi,i＝||Qi T||1×1025；其中，Pi,i、Qi,i分别为系数矩阵Pi、Qi的第i个元素，上角标T为矩阵中的转置符号；
所述步骤九中组装步骤八引入本质边界条件的线性方程得到线性方程组即将步骤八系数矩阵、右边项矩阵装配成整体矩阵；若是由步骤8.1或步骤8.2直接跳转过来的则组装步骤七的所有离散点处位移场变量值的线性方程得到线性方程组；
组装后，有大小为2N×2N的整体系数矩阵K＝[P1 T Q1 T P2 T Q2 T … … PN T QN T]T，有大小为2N×1的右边项矩阵Z＝[R1 T R2 T … RN T]T；故有以所有离散点处位移场变量为未知的线性方程组：K·U＝Z。
所述步骤十一断裂力学参数计算，具体包括以下步骤：
步骤1，提取步骤五目标离散点i处导数关于局部离散点场变量的近似表达式；
步骤2，将目标离散点i代入场变量的具体含义，即x方向位移场或y方向位移场，根据步骤十计算得到的离散点位移场变量值，得到目标离散点i处x方向位移场、y方向位移场的导数近似值；
步骤3，将步骤2位移场的导数近似值代入弹性力学表达式，得到目标离散点i的应力、应变场量值；
步骤4，根据断裂力学知识，计算得到断裂力学参数：应力强度因子、应变能密度、最大环向应力、最大能量释放率；
所述步骤十三具体包括以下步骤：
步骤1，根据步骤十二裂纹扩展准则和步骤十一断裂力学参数，确定裂纹起裂点；若有多个裂纹起裂点，则对每个裂纹起裂点执行步骤2至步骤5；
步骤2，根据裂纹扩展准则，计算裂纹起裂点的裂纹扩展方向；
步骤3，记录距离裂纹起裂点最近的离散点与裂纹起裂点之间的距离，记为L0，在距离起裂点
Figure GDA0004246385180000071
方向为步骤2裂纹扩展方向的地方生成新的离散点，若距离新生离散点最近的离散点与新生离散点之间的距离不超过L0×10-5，则认为该离散点与新生离散点重合；若超过L0×10-5，则认为无重合；
步骤4，根据步骤2裂纹扩展方向和步骤3裂纹起裂点、新生离散点，对应生成裂纹节点对，并对系列新生离散点进行编号，并记录坐标信息，纳入步骤一离散点编号及坐标信息的集合中去，完成离散点集合的更新；
步骤5，将裂纹起裂点与新生离散点连线组成的线条纳入步骤二逆时针边界描述及边界条件中去，完成边界条件的更新；同时，记录裂纹起裂点与新生离散点的连线情况，此为裂纹扩展路径。
本发明的有益效果：
1.本发明所述数值离散算子可配合各类数值分析方法使用，以满足高精度后处理的需求，亦可单独使用，为偏微分方程求解提供新思路，灵活性与应用性强，可进一步扩展应用至连续、非连续问题中去；
2.本发明所述断裂模拟方法通过节点生死算法，可方便地模拟裂纹向任意方向、任意长度的扩展问题，可应用于不同扩展准则下各类型、各维度的裂纹扩展研究；
3.本发明所述断裂模拟方法中各方程有明确物理意义：块体平衡，故可进一步扩展应用至流固耦合问题中去，应用范围广，并可对发明中平衡方程做调整，适用于各种体积力情况中去；
4.本发明所述方法可随时切换有网格和无网格，实现有网格和无网格的二象性，可为解决现有算法模拟裂纹不足提供新方向，发展前景广阔。
附图说明
图1为本发明中圆形积分域定义方法示意图；
图2为本发明关于裂纹的描述示意图；
图3为实施例1目标离散点与局部离散点集合的坐标关系示意图；
图4为实施例2运算域示意图；
图5为实施例3模型尺寸及边界条件示意图；
图6为实施例3离散点分布示意图；
图7为实施例3积分域搜索定义示例；
图8为实施例3初始模型直接计算结果——X方向位移场云图；
图9为实施例3初始模型直接计算结果——Y方向位移场云图；
图10为实施例3初始模型间接结果——σxx应力场云图；
图11为实施例3初始模型间接结果——σyy应力场云图；
图12为实施例3初始模型间接结果——σxy应力场云图；
图13为实施例3裂纹扩展路径预测结果。
具体实施方式
为使本发明实施的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行更加详细的描述。
需要说明的是：下面通过参考附图描述的实施例都是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。所描述的实施例是本发明一部分实施例，而不是全部的实施例，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。基于本发明中的实施例，相关领域技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
具体地，一种基于积分-广义有限差分数值离散算子的断裂模拟方法，包括以下步骤：
步骤1，执行空间离散算法，将研究域空间离散化，获得研究域空间内一系列离散点的编号及其坐标(xi,yi),i＝1,2,...,N；
空间离散算法——目标是在研究域内布点，该步骤可通过两种方式完成：
方式1，调用外源软件，如Gmsh、Ansys、Abaqus、HyperMesh软件进行离散：
步骤1.1.1，向外源软件内输入研究域的边界信息；
步骤1.1.2，利用外源软件完成研究域的网格剖分；
步骤1.1.3，输出网格点的编号及坐标信息，得到网格点的编号及坐标信息即为研究域内离散点的编号及其坐标(xi,yi),i＝1,2,...,N；
方式2，采用内置算法，简单布点：以二维图形为例，
步骤1.2.1，根据研究域的边界信息，确定边界的下限、上限、左限、右限，即ymin、ymax、xmin、xmax，围成最大整合区域；
步骤1.2.2，确定最大整合区域内x、y轴等分数量Nx、Ny，并将最大整合区域等分，得到
Figure GDA0004246385180000095
个交点(即等分线以及边界的下限、上限、左限、右限之间的全部交点)；
步骤1.2.3，剔除位于研究域之外的交点，剩余交点即为研究域内的离散点，对其(N个点)进行编号并记录坐标信息(xi,yi),i＝1,2,...,N；
步骤2，利用步骤1的离散点对边界线条进行逆时针描述，形成闭合回路，限定各边界线条的边界条件；
逆时针描述——将位于边界上的离散点进行逆时针连线，以离散点的编号对该逆时针闭合回路进行描述。
所述边界条件分为自然边界条件、本质边界条件、混合边界条件，根据所研究问题具体确定。
步骤3，确定的是平衡方程的等效积分弱形式；
根据弹性力学，有忽略体积力的高阶平衡方程(当体积力不为零时可做类似推导)：
Figure GDA0004246385180000091
其中，平面应变问题
Figure GDA0004246385180000092
平面应力问题
Figure GDA0004246385180000093
其中E为弹性模量，v为泊松比，A为微元体的面积，w1和w2分别为x、y方向上的位移；
根据格林公式(式(4))，实现对高阶平衡方程的降阶，最终形式如式(5)所示。
Figure GDA0004246385180000094
式中，L是围成区域Ω的闭合逆时针回路，P、Q是任意一个在区域Ω内连续的函数；
Figure GDA0004246385180000101
式中，
Figure GDA0004246385180000102
是围成式(1)中微元面积A的逆时针闭合回路，(nx,ny)是闭合回路/>
Figure GDA0004246385180000103
中各线条的外矢量方向。
至此，面积分(式(1))转化为线积分(式(5))。进一步地，假设式(5)中①②③④项在各线条上为线性分布，则线积分可以转化为关于线条始末两点的积分形式，如式(6)所示。
Figure GDA0004246385180000104
其中，LS是线条S的长度，sp、ep分别代表线条S的起始点和终点。
将式(6)代入式(5)中，得到平衡方程的等效积分弱形式(式(7))：
Figure GDA0004246385180000105
式中，Ln是围绕微元面积A的逆时针闭合回路中线条的数量，LSj是第j(j＝1,2,...,Ln)条边界线条的长度，sp,j、ep,j分别代表第j条边界线条的起始点和终点,(nx,j,ny,j)是第j条边界线条的外矢量方向。
步骤4，对步骤1所得的每个离散点i(i＝1,2,...,N)均按照步骤4-步骤8操作(即循环离散点)：对离散点i(i＝1,2,...,N)执行积分域搜索算法，记录每个离散点的积分域信息；
循环——对所有离散点执行至“循环结束”的步骤内容。
积分域搜索算法——目标是确定每个离散点的积分域，该步骤可通过两种方式完成：
方式1，单元覆盖方法——根据生成单元信息输出所需积分域：根据步骤1离散点的编号和坐标信息，生成Delaunay三角形单元，包含离散点i的Delaunay三角形单元即为该离散点的积分域；
方式2，随机覆盖方法——随机生成合适大小的积分域并输出，具体包括以下步骤：
步骤4.2.1，根据研究域类型确定积分域的类型，当研究域为长条形时，选择长方形或椭圆形为积分域类型，当研究域为圆形时，选择圆形或正方形作为积分域类型，进一步再确定积分域的形状参数，参考附图1，以圆形积分域形成为例；
步骤4.2.2，根据步骤1离散点的编号和坐标信息，筛选出距离离散点i小于或者等于步骤4.2.1积分域形状参数的积分域内的局部离散点；
步骤4.2.3，根据步骤4.2.2积分域内的局部离散点，其中属于步骤1研究域内的离散点i的积分域作为该离散点的积分域(附图1中灰色部分)；并进行记录；
步骤5，根据步骤1离散点坐标信息与步骤2逆时针边界描述，执行数值离散算法，将目标离散点i的场量值及其导数近似表达为关于局部离散点(n个点)场变量的关系式，具体内容如下；
所述数值离散算法，具体包括以下步骤：
步骤5.1，将目标离散点i的场变量值在局部离散点(n个点)处进行泰勒展开，并根据用户目的确定泰勒展开式截断项数，写成矩阵形式，得到式(11)、式(13)的具体形式；
步骤5.2，引入移动最小二乘方法，选择合适的加权函数类型，根据目标离散点与局部离散点坐标关系建立加权矩阵(即构建式(12))；
步骤5.3，根据步骤5.1截断项数确定比例矩阵内容(式(16))；
步骤5.4，根据步骤5.1矩阵形式、步骤5.2加权矩阵、步骤5.3比例矩阵，构建目标离散点i的场量值及其导数关于局部离散点场量的近似表达式(式(15)、式(17))。
具体如下：设目标离散点i的周围有n个局部离散点，有泰勒完全展开式(式(8))：
Figure GDA0004246385180000111
其中，△xj＝xj-x0，△yj＝yj-y0，(xj,yj)是第j个局部离散点的坐标，u0、
Figure GDA0004246385180000112
Figure GDA0004246385180000113
为目标离散点i(x0,y0)处的场变量和各阶导数，分别简写为u0、u,x|0、u,y|0、u0,xy|0、u,xx|0、u,yy|0；ui为第j个局部离散点处的场变量真值；式(8)中未知量包括目标离散点的场变量及各阶导数。在本发明中，所述场变量的物理意义包括位移、应力、应变、断裂力学参数，在后续步骤中，根据不同的需求赋予其不同的实际意义；
假设式(8)保留了前K项，为使第j个局部离散点处场变量的精确值
Figure GDA0004246385180000121
与通过式(8)计算得到的近似值/>
Figure GDA0004246385180000122
(公式(8)中是带有省略号的，是无穷尽的，所以可以说是真值，等号成立；当实际求解过程中，只能写出有穷尽项，所以严格来说真值约等于右边，也就是说，通过公式8，保留前K项，计算得到的是近似值/>
Figure GDA0004246385180000123
)之间的总体偏差/>
Figure GDA0004246385180000124
最小，采用移动最小二乘法求解，要求：
Figure GDA0004246385180000125
其中，Wj为移动最小二乘方法中第j个局部离散点的不为零的加权值，通过用户确定的加权函数类型确定，进一步地，式(9)可写成式(10)的矩阵形式：
Figure GDA0004246385180000126
其中，
Figure GDA0004246385180000127
W＝diag([W1 W2 … Wn]) (12)
Figure GDA0004246385180000128
Figure GDA0004246385180000131
为便于书写，令G＝ST·W·S，H＝ST·W，代入式(10)，有
Figure GDA0004246385180000132
为避免矩阵G的病态与奇异，引入比例矩阵
Figure GDA0004246385180000133
其中，|△x|max＝max{|△x1|,|△x2|,...,|△xn|}，|△y|max＝max{|△y1|,|△y2|,...,|△yn|}，比例矩阵修正后，有如下表达式：
Figure GDA0004246385180000134
/>
从而能够确定任意目标离散点i的场变量u0及其各阶导数关于局部离散点(n个点)场变量的近似表达式；
步骤6，赋予步骤5所有离散点场变量为位移的物理含义，根据步骤2边界条件，执行自然边界条件引入算法，得到引入自然边界条件的平衡方程等效积分弱形式；(平衡方程等效积分弱形式在步骤3中就已经得到了一个初步的表达式，在后续步骤中，引入自然边界条件或是本质边界条件，都属于细微调整)
自然边界条件引入——引入方式同有限单元法，具体包括以下步骤：
步骤6.1，判断目标离散点i的积分域是否包含步骤2中边界线条，即：对于目标离散点i，根据步骤4积分域信息，判断其是否包含步骤2中边界线条，若是，执行步骤6.2；若否，跳至步骤7；
步骤6.2，进一步判断步骤6.1中包含的边界线条是否为自然边界条件，若是，执行步骤6.3；若否，跳至步骤7；
步骤6.3，将步骤3得到的式(7)中关于该边界线条的表达移动到等号右侧，得到引入自然边界条件的平衡方程等效积分弱形式；
步骤7，将步骤5目标离散点i处导数关于局部离散点(n个点)位移场变量的关系式代入步骤6.3引入自然边界条件的平衡方程的等效积分弱形式，若是由步骤6.1或步骤6.2直接跳至本步骤的则将步骤5目标离散点i处导数关于局部离散点(n个点)位移场变量的关系式代入步骤3的平衡方程的等效积分弱形式，扩充形成关于所有离散点(N个点)处位移场变量值的线性方程；
扩充——通过添加零元素，将以n个局部离散点位移为未知量的表达式，写成以N个离散点位移为未知量的线性方程，具体过程如下：
假设以n个局部离散点位移场变量为未知的表达式有如下形式：
Figure GDA0004246385180000141
其中，上标1、2分别代表x方向、y方向，Nj是第j(＝1,2,...,n)个局部离散点在所有离散点中的编号，即
Figure GDA0004246385180000142
代表编号为Nj的离散点在x方向、y方向上的位移未知量，Pi、Qi、Ri,P、Ri,Q为将步骤5目标离散点i处导数关于局部离散点(n个点)场变量的关系式代入步骤3平衡方程的等效积分弱形式所得到的线性方程系数；经扩充之后，有大小为1×2N的系数矩阵Pi，Qi，第1行、第2Nj-1列的元素分别为/>
Figure GDA0004246385180000143
第1行、第2Nj列的元素分别为
Figure GDA0004246385180000144
Figure GDA0004246385180000145
而其他元素为零；有大小为2×1的右边项矩阵Ri＝[Ri,P Ri,Q]T；有大小为2N×1的未知量矩阵U＝[u1 1 u1 2 u2 1 u2 2 … … uN 1 uN 2]T，即式(18)可以写成式(19)的形式：
Figure GDA0004246385180000146
步骤8，根据步骤2边界条件，对步骤7得到的线性方程进行修正，执行本质边界条件引入算法，得到引入本质边界条件的线性方程；
本质边界条件引入——通过乘大数法引入，具体包括以下步骤：
步骤8.1，对于目标离散点i，判断其是否包含步骤2中边界线条，若是，执行步骤8.2，若否，跳至步骤9；
步骤8.2，进一步判断步骤8.1所包含边界线条是否为本质边界条件，若是，执行步骤8.3；若否，跳至步骤9；
步骤8.3，若在目标离散点i处施加数值为DB(DB≠0)、方向为x方向的本质边界条件,则令步骤7所得线性方程中的Pi,i＝||Pi T||1×1025，Ri,P＝DB×||Pi T||1×1025，|| ||1为一阶范数；若在目标离散点i处施加数值为DB(DB≠0)、方向为y方向的本质边界条件,则令步骤7所得线性方程中的Qi,i＝||Qi T||1×1025，Ri,Q＝DB×||Ri T||1×1025；若在目标离散点i处施加数值为0、方向为x方向的本质边界条件,则令步骤7所得线性方程中的Pi,i＝||Pi T||1×1025；若在目标离散点i处施加数值为0、方向为y方向的本质边界条件,则令步骤7所得线性方程中的Qi,i＝||Qi T||1×1025；其中，Pi,i、Qi,i分别为系数矩阵Pi、Qi的第i个元素；Pi T、Qi T、Ri T分别是式(19)中Pi、Qi、Ri的转置，即上角标T为矩阵中的转置符号。
步骤9，循环结束，组装步骤8引入本质边界条件的线性方程得到线性方程组，即将步骤8系数矩阵、右边项矩阵装配成整体矩阵；若是由步骤8.1或步骤8.2直接跳转过来的则组装步骤七的所有离散点处位移场变量值的线性方程得到线性方程组，
组装后，有大小为2N×2N的整体系数矩阵K＝[P1 T Q1 T P2 T Q2 T … … PN T QN T]T，有大小为2N×1的右边项矩阵Z＝[R1 T R2 T … RN T]T；
故有以所有离散点(N个点)处位移场变量为未知的线性方程组(式(20))：
K·U＝Z (20)
步骤10，按照下式计算，得到所有离散点(N个点)的位移场变量值；
根据式(20)，有
U＝K-1·Z (21)
则可以计算得到离散点的位移场量值；
步骤11，断裂力学参数计算，具体包括以下步骤：
步骤11.1，循环离散点，提取步骤5目标离散点i处导数关于局部离散点场变量的近似表达式；
步骤11.2，将目标离散点i代入场变量的具体含义，如x方向位移场、y方向位移场，及步骤10离散点位移场变量值，得到目标离散点i处x方向位移场、y方向位移场的导数近似值；
步骤11.3，将步骤11.2位移场的导数近似值代入弹性力学表达式，得到目标离散点i的应力、应变场量值；
步骤11.4，进一步地，根据断裂力学知识，计算得到断裂力学参数：应力强度因子、应变能密度、最大环向应力、最大能量释放率；
步骤12，裂纹扩展准则判断：根据步骤11断裂力学参数，若存在离散点的断裂力学参数满足裂纹扩展准则，执行步骤13，并重复步骤3至步骤12，每一次重复为一次计算步；若研究域内所有离散点的断裂力学参数不满足裂纹扩展准则，说明该状态下研究域已达到变形稳定阶段，裂纹不再扩展，跳至步骤14；
步骤13，节点生死算法，离散点集合和边界条件的更新，裂纹扩展行为预测，具体包括以下步骤：
步骤13.1，根据步骤12裂纹扩展准则和步骤11断裂力学参数，确定裂纹起裂点；若有多个裂纹起裂点，则对每个裂纹起裂点执行步骤13.2至步骤13.5；
步骤13.2，根据裂纹扩展准则，计算裂纹起裂点的裂纹扩展方向；(根据不同的断裂力学参数和裂纹扩展准则，有不同的计算公式，具体到某一断裂力学参数和裂纹扩展准则，其计算方法为本领域人员公知的内容)
步骤13.3，记录距离裂纹起裂点最近的离散点与裂纹起裂点之间的距离，记为L0，在距离起裂点
Figure GDA0004246385180000161
方向为步骤13.2裂纹扩展方向的地方生成新的离散点，若距离新生离散点最近的离散点与新生离散点之间的距离不超过L0×10-5，则认为该离散点与新生离散点重合；若超过L0×10-5，则认为无重合；
步骤13.4，如附图2所示，根据步骤13.2裂纹扩展方向和步骤13.3裂纹起裂点、新生离散点，对应生成裂纹节点对，并对系列新生离散点进行编号，并记录坐标信息，纳入步骤1离散点编号及坐标信息的集合中去，完成离散点集合的更新；(裂纹节点对，参考附图2，是将新生离散点垂直于裂纹线方向偏移一定距离得到的)
步骤13.5，将裂纹起裂点与新生离散点连线组成的线条纳入步骤2逆时针边界描述及边界条件中去，完成边界条件的更新；同时，记录裂纹起裂点与新生离散点的连线情况，此为裂纹扩展路径；
步骤14，计算结果的可视化与后处理；
通过步骤12、步骤13，得到每一计算步下的裂纹扩展情况，从而实现对裂纹扩展路径的追踪；
在步骤14将每一计算步下的步骤10离散点的位移场量值、步骤11离散点应力场量值、应变场量值、断裂力学参数，及直至达到变形稳定阶段时的裂纹扩展路径可视化或进行其他后处理，以满足人员特别需要。
本发明所提供的一种基于积分-广义有限差分数值离散算子的断裂模拟方法，具有如下特点：
1.本发明所述数值离散算子联合泰勒级数展开技术和移动最小二乘技术，获得离散算子的近似表达，克服了对网格的依赖性，计算原理简单，易于实现，具有高精度、高计算效率的特点；
2.本发明结合格林公式或变分原理，实现高阶微分算子的降阶，从而降低了对场函数本身连续性的要求；
3.本发明所述断裂模拟方法无需先验引入基于断裂力学解的增强函数；
4.本发明所述断裂模拟方法可方便地通过乘大数法或罚函数法引入自然边界条件，克服了移动最小二乘近似的缺点，本质边界条件引入方式同有限单元法，即：边界条件的引入，方便快捷且精准；
5.本发明所述断裂模拟方法通过节点生死实现裂纹及其扩展过程的追踪，简单易施行；通过节点对来描述裂纹，将裂纹面分为上下两部分，以节点对垂直于裂纹线方向的相对位移描述裂纹开度，以节点对平行于裂纹线方向相对位移描述上下裂纹面的相互错动，可模拟张拉型、剪切型、混合型裂纹扩展行为；
6.本发明所述断裂模拟方法以等效积分弱形式平衡方程为基础组装整体方程组，物理意义为积分域即块体的平衡，而非传统有限元的点的平衡。
为更好地解释本发明的可实施性，做如下实施例：
实施例1：单点一阶导数的数值离散格式
假定目标点坐标为(1.85,1.675)。其局部离散点集合为(0,0)、(2,0.5)、(2.25,2)、(1,2)、(2,0)、(1.5,2.5)、(0,1)、(0,-1)，共8个，如附图3所示。确定截断项为6项，故有泰勒展开式
Figure GDA0004246385180000171
将目标点坐标与周围散点坐标代入式(11)中，得到泰勒展开式的矩阵系数
Figure GDA0004246385180000172
Figure GDA0004246385180000173
选择四次样条曲线(式(23)，
Figure GDA0004246385180000174
△dmax＝α·max{△d1,△d2,...,△d8}，本实施例中取α＝1.2)为加权函数，代入目标点与周围离散点坐标，根据式(12)得到加权矩阵
W＝diag{[0.0162 0.8377 0.9928 0.9368 0.5305 0.9933 0.9975 0.6432]}
并根据式(16)有比例矩阵
D＝diag{[1 0.5405 0.5970 0.2922 0.3564 0.3227]}
将以上矩阵代入式(17)中，即可得到G、H的数学表达式，进一步地，有如下目标点一阶导数的离散近似格式：
Figure GDA0004246385180000181
若假设场分布形式为u＝x2+y，即
Figure GDA0004246385180000182
代入式(24)中，得到
Figure GDA0004246385180000183
与解析解(分别为3.70和1)一致。
实施例2：单点二阶微分算子的数值离散格式
计算模型同实施例1。
对于在散点的二阶微分算子
Figure GDA0004246385180000184
用它在运算域R0的平均值来近似，并运用格林公式有：
Figure GDA0004246385180000185
假设
Figure GDA0004246385180000186
在运算域R0每一条围线上为线性分布，则可进一步写成
Figure GDA0004246385180000187
其中，m为运算域R0的围线数量，sp和ep代表第i(i＝1,2,…,m)条围线的起点与终点。
附图4给出了实施例2的运算域示意图。同实施例1所述步骤，构建各围线起点、终点的离散近似格式，后代入式(26)中得到目标点二阶微分算子的离散表达式：
Figure GDA0004246385180000191
若假设场分布形式为u＝x2+y，即
Figure GDA0004246385180000192
代入式(27)中，得到
Figure GDA0004246385180000193
/>
与解析解(分别为2和0)一致。
实施例3：有限宽板中边缘裂纹受单向拉伸分布载荷的作用
假定该板为平面应力状态，模型尺寸及边界条件如图5所示。首先利用空间离散算法，对研究域进行空间上的离散，并对研究域的边界进行逆时针闭合回路表达，离散点分布情况见图6。假定材料参数为：弹性模型E＝1000，泊松比v＝0.3。
之后根据离散点的坐标和边界描述，选择加权函数类型为四次样条曲线(式(23)，其中
Figure GDA0004246385180000194
△dmax＝α·max{△d1,△d2,...,△dm}，本实施例中取α＝1.2)，循环离散点，根据式(8)-式(17)构建各离散点场量值及其导数关于局部离散点场量的近似表达式。
根据离散点的坐标和边界描述，构建Delaunay三角形背景网格，为积分域搜索算法做准备工作。图7给出了一个积分域定义的示例。
将构建的离散近似表达式，结合积分域，代入平衡方程的等效积分弱形式中去，可以得到一个形如
Ax＝b (28)的以离散点处位移场量值为未知量的线性方程组。应注意，过程中如碰到施加自然边界条件的边界线条，则将其移动到等号右侧，以此实现自然边界条件的引入。
之后，对需要施加本质边界条件的离散点，通过乘大数法对线性方程组进行处理，精准施加本质边界条件。
然后就可以进行求解运算，式(28)所得到的直接结果为各离散点位移场量值，如图8-9所示。根据弹性力学的本构方程和前述构建的离散近似表达式，可得到离散点处应力场量值，如图10-12所示。
在本实施例中，以最大周向应力为裂纹扩展准则，即：裂纹将沿着σθθmax所对应的θ的方向扩展；设定每计算步裂纹扩展长度为0.05。最终得到的裂纹扩展路径如图13所示。

Claims (6)

1.一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于，使用积分-广义有限差分方法离散并求解线弹性力学问题，通过断裂准则对裂纹行为进行预测，并通过节点生死实现裂纹及其扩展过程的追踪；
具体包括如下步骤：
步骤一：执行空间离散算法，将研究域空间离散化，获得研究域空间内一系列离散点的编号及其坐标；
步骤二：利用步骤一离散点对边界线条进行逆时针描述，形成闭合回路，限定各边界线条的边界条件；
步骤三：确定平衡方程的等效积分弱形式；然后对步骤一所得的每个离散点均按照步骤四-步骤八操作；
步骤四：对离散点执行积分域搜索算法，确定并记录每个离散点的积分域信息；
步骤五：根据步骤一离散点坐标信息与步骤二逆时针边界描述，执行数值离散算法将离散点的场变量值及其导数近似表达为关于局部离散点场变量的关系式；
步骤六：赋予步骤五场变量为位移的物理含义，根据步骤二边界条件，执行自然边界条件引入算法，得到引入自然边界条件的平衡方程等效积分弱形式；
步骤七：将步骤五离散点处导数关于局部离散点位移场变量的关系式代入步骤六引入自然边界条件的平衡方程等效积分弱形式，扩充形成关于所有离散点处位移场变量值的线性方程；
步骤八：根据步骤二边界条件，执行本质边界条件引入算法，得到引入本质边界条件的线性方程；
步骤九：组装步骤八引入本质边界条件的线性方程得到线性方程组；
步骤十：执行计算，得到所有离散点的位移场变量值；
步骤十一：断裂力学参数计算；
步骤十二：裂纹扩展准则判断：根据步骤十一断裂力学参数，
若离散点的位移满足裂纹扩展准则，执行步骤十三，并将满足扩展准则的离散点重复步骤三至步骤十二，每一次重复为一次计算步；
若离散点的位移不满足裂纹扩展准则，说明该离散点所在的研究区域已达到变形稳定阶段，裂纹不再扩展，跳至步骤十四；
步骤十三：节点生死算法，离散点集合和边界条件的更新，裂纹扩展行为预测；
步骤十四：计算结果的可视化与后处理：将每一计算步下的步骤十离散点的位移场变量值、步骤十一离散点应力场变量值、应变场变量值、断裂力学参数，及直至达到变形稳定阶段时的裂纹扩展路径可视化或进行其他后处理，以满足人员特别需要；
所述步骤一是通过调用外源软件或采用内置算法进行的；
其中通过调用外源软件具体包括如下步骤：
步骤1，向外源软件内输入研究域的边界信息；
步骤2，利用外源软件完成研究域的网格剖分；
步骤3，输出网格点的编号及坐标信息，得到的网格点编号及坐标信息即为研究域内离散点的编号及其坐标(xi,yi),i＝1,2,...,N；
其中采用内置算法具体包括如下步骤：
步骤1，根据研究域的边界信息，确定边界的下限、上限、左限、右限，即ymin、ymax、xmin、xmax，围成最大整合区域；
步骤2，确定最大整合区域内x、y轴等分数量Nx、Ny，并将最大整合区域等分，得到
Figure FDA0004246385170000021
个交点；
步骤3，剔除位于研究域之外的交点，剩余交点即为研究域内的离散点，对其进行编号并记录坐标信息(xi,yi),i＝1,2,...,N；
所述步骤三确定的平衡方程的等效积分弱形式，针对忽略体积力的线弹性力学问题中，具体如下：
Figure FDA0004246385170000031
式中，Ln是围绕微元面积A的逆时针闭合回路中线条的数量，LSj是第j(j＝1,2,...,Ln)条边界线条的长度，sp,j、ep,j分别代表第j条边界线条的起始点和终点,(nx,j,ny,j)是第j条边界线条的外矢量方向，
Figure FDA0004246385170000032
Figure FDA0004246385170000033
式中，针对平面应变问题有/>
Figure FDA0004246385170000034
对平面应变问题E'＝E，ν'＝ν，E为弹性模量，v为泊松比，w1和w2分别为x、y方向上的位移；
所述步骤十一断裂力学参数计算，具体包括以下步骤：
步骤1，提取步骤五目标离散点i处导数关于局部离散点场变量的近似表达式；
步骤2，将目标离散点i代入场变量的具体含义，即x方向位移场或y方向位移场，根据步骤十计算得到的离散点位移场变量值，得到目标离散点i处x方向位移场、y方向位移场的导数近似值；
步骤3，将步骤2位移场的导数近似值代入弹性力学表达式，得到目标离散点i的应力、应变场量值；
步骤4，根据断裂力学知识，计算得到断裂力学参数：应力强度因子、应变能密度、最大环向应力、最大能量释放率；
所述步骤十三具体包括以下步骤：
步骤1，根据步骤十二裂纹扩展准则和步骤十一断裂力学参数，确定裂纹起裂点；若有多个裂纹起裂点，则对每个裂纹起裂点执行步骤2至步骤5；
步骤2，根据裂纹扩展准则，计算裂纹起裂点的裂纹扩展方向；
步骤3，记录距离裂纹起裂点最近的离散点与裂纹起裂点之间的距离，记为L0，在距离起裂点
Figure FDA0004246385170000041
方向为步骤2裂纹扩展方向的地方生成新的离散点，若距离新生离散点最近的离散点与新生离散点之间的距离不超过L0×10-5，则认为该离散点与新生离散点重合；若超过L0×10-5，则认为无重合；
步骤4，根据步骤2裂纹扩展方向和步骤3裂纹起裂点、新生离散点，对应生成裂纹节点对，并对系列新生离散点进行编号，并记录坐标信息，纳入步骤一离散点编号及坐标信息的集合中去，完成离散点集合的更新；
步骤5，将裂纹起裂点与新生离散点连线组成的线条纳入步骤二逆时针边界描述及边界条件中去，完成边界条件的更新；同时，记录裂纹起裂点与新生离散点的连线情况，此为裂纹扩展路径。
2.根据权利要求1所述的一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于所述步骤四是通过单元覆盖方法或随机覆盖方法进行的；
其中，单元覆盖方法为根据步骤一离散点的编号和坐标信息，生成Delaunay三角形单元，包含离散点的Delaunay三角形单元即为该离散点的积分域；
其中，随机覆盖方法具体包括以下步骤：
步骤1，根据研究域类型确定积分域的类型，当研究域为长条形时，选择长方形或椭圆形为积分域类型，当研究域为圆形时，选择圆形或正方形作为积分域类型，再确定积分域的形状参数；
步骤2，根据步骤一离散点的编号和坐标信息，筛选出距离离散点小于或者等于步骤1积分域形状参数的积分域内的局部离散点；
步骤3，根据步骤2积分域内的局部离散点，其中属于步骤一研究域内的离散点的积分域作为该离散点的积分域。
3.根据权利要求2所述的一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于所述步骤五中目标离散点i的场量值及其导数近似表达为关于局部离散点场变量的关系式如下：
G＝DT·ST·W·S
H＝DT·ST·W
Figure FDA0004246385170000051
其中：
Figure FDA0004246385170000052
|△x|max＝max{|△x1|,|△x2|,...,|△xn|}，|△y|max＝max{|△y1|,|△y2|,...,|△yn|}，T代表转置，
Figure FDA0004246385170000053
W＝diag([W1 W2 … Wn])，
Figure FDA0004246385170000054
目标离散点i的周围有n个局部离散点，(xj,yj)是第j个局部离散点的坐标，j＝1、2、……n，△xj＝xj-x0，△yj＝yj-y0，u0、
Figure FDA0004246385170000055
Figure FDA0004246385170000056
为目标离散点i(x0,y0)处的场变量和各阶导数。
4.根据权利要求3所述的一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于所述步骤六的具体过程如下：
步骤6.1，判断目标离散点i的积分域是否包含步骤二中边界线条，若是，执行步骤6.2；若否，跳至步骤七；
步骤6.2，进一步判断步骤6.1中包含的边界线条是否为自然边界条件，若是，执行步骤6.3；若否，跳至步骤七；
步骤6.3，将步骤三得到平衡方程的等效积分弱形式中关于该边界线条的表达移动到等号右侧，得到引入自然边界条件的平衡方程等效积分弱形式。
5.根据权利要求4所述的一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于所述步骤七将步骤五目标离散点i处导数关于局部离散点位移场变量的关系式代入步骤6.3引入自然边界条件的平衡方程的等效积分弱形式，若是由步骤6.1或步骤6.2直接调至本步骤的则将步骤五目标离散点i处导数关于局部离散点位移场变量的关系式代入步骤三的平衡方程的等效积分弱形式，扩充形成关于所有离散点处位移场变量值的线性方程；
其中扩充形成关于所有离散点处位移场变量值的线性方程为：
Figure FDA0004246385170000061
其中：Pi，Qi为经扩充之后大小为1×2N的系数矩阵，Ri＝[Ri,P Ri,Q]T为扩充后大小为2×1的右边项矩阵，U＝[u1 1 u1 2 u2 1 u2 2 … … uN 1 uN 2]T为扩充后大小为2N×1的未知量矩阵。
6.根据权利要求5所述的一种基于积分-广义有限差分数值离散算子的断裂模拟方法，其特征在于所述步骤八的具体过程包括如下步骤：
步骤8.1，对于目标离散点i，判断其是否包含步骤二中边界线条，若是，执行步骤8.2，若否，跳至步骤九；
步骤8.2，进一步判断步骤8.1所包含边界线条是否为本质边界条件，若是，执行步骤8.3；若否，跳至步骤九；
步骤8.3，若在目标离散点i处施加数值为DB(DB≠0)、方向为x方向的本质边界条件,则令步骤七所得线性方程中的Pi,i＝||Pi T||1×1025，Ri,P＝DB×||Pi T||1×1025，||||1为一阶范数；若在目标离散点i处施加数值为DB(DB≠0)、方向为y方向的本质边界条件,则令步骤七所得线性方程中的Qi,i＝||Qi T||1×1025，Ri,Q＝DB×||Ri T||1×1025；若在目标离散点i处施加数值为0、方向为x方向的本质边界条件,则令步骤七所得线性方程中的Pi,i＝||Pi T||1×1025；若在目标离散点i处施加数值为0、方向为y方向的本质边界条件,则令步骤七所得线性方程中的Qi,i＝||Qi T||1×1025；其中，Pi,i、Qi,i分别为系数矩阵Pi、Qi的第i个元素，上角标T为矩阵中的转置符号；
所述步骤九中组装步骤八引入本质边界条件的线性方程得到线性方程组即将步骤八系数矩阵、右边项矩阵装配成整体矩阵；若是由步骤8.1或步骤8.2直接跳转过来的则组装步骤七的所有离散点处位移场变量值的线性方程得到线性方程组；
组装后，有大小为2N×2N的整体系数矩阵K＝[P1 T Q1 T P2 T Q2 T … … PN T QN T]T，有大小为2N×1的右边项矩阵Z＝[R1 T R2 T … RN T]T；故有以所有离散点处位移场变量为未知的线性方程组：K·U＝Z。
CN202111044179.7A 2021-09-07 2021-09-07 一种基于积分-广义有限差分数值离散算子的断裂模拟方法 Active CN113722965B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111044179.7A CN113722965B (zh) 2021-09-07 2021-09-07 一种基于积分-广义有限差分数值离散算子的断裂模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111044179.7A CN113722965B (zh) 2021-09-07 2021-09-07 一种基于积分-广义有限差分数值离散算子的断裂模拟方法

Publications (2)

Publication Number Publication Date
CN113722965A CN113722965A (zh) 2021-11-30
CN113722965B true CN113722965B (zh) 2023-06-27

Family

ID=78682187

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111044179.7A Active CN113722965B (zh) 2021-09-07 2021-09-07 一种基于积分-广义有限差分数值离散算子的断裂模拟方法

Country Status (1)

Country Link
CN (1) CN113722965B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115292990B (zh) * 2022-07-18 2023-04-11 南方科技大学 一种连续-非连续耦合的二维固体破裂模拟方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108197358A (zh) * 2017-12-20 2018-06-22 北京石油化工学院 一种高效快速模拟水力压裂的方法
CN110555229A (zh) * 2019-07-12 2019-12-10 北京航空航天大学 一种无网格固体力学仿真方法、电子设备及存储介质
CN110779795A (zh) * 2019-11-04 2020-02-11 中国石油大学(华东) 裂缝性储层地质力学建模网格单元大小确定方法
CN111460568A (zh) * 2020-04-22 2020-07-28 山西省河道与水库技术中心 一种混凝土重力坝运行期裂缝扩展判别方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7822585B2 (en) * 2005-05-03 2010-10-26 Hrl Laboratories, Llc System and a method for numerical simulation of wave propagation in homogeneous media

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108197358A (zh) * 2017-12-20 2018-06-22 北京石油化工学院 一种高效快速模拟水力压裂的方法
CN110555229A (zh) * 2019-07-12 2019-12-10 北京航空航天大学 一种无网格固体力学仿真方法、电子设备及存储介质
CN110779795A (zh) * 2019-11-04 2020-02-11 中国石油大学(华东) 裂缝性储层地质力学建模网格单元大小确定方法
CN111460568A (zh) * 2020-04-22 2020-07-28 山西省河道与水库技术中心 一种混凝土重力坝运行期裂缝扩展判别方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
基于单位分解的扩展径向点插值无网格法;马文涛 等;岩土力学;第33卷(第12期);全文 *
岩土体损伤局部化及其灾变渐进破坏数值模拟研究;张小强;中国博士学位论文全文数据库 基础科学辑;全文 *
水泥土桩加固边坡变形破坏机理与稳定性研究;许胜才;中国博士学位论文全文数据库 工程科技Ⅱ辑;全文 *
考虑块体内耦合作用及锚杆塑性效应的非连续变形分析方法;王文;中国博士学位论文全文数据库 工程科技Ⅱ辑;全文 *
裂纹扩展的无网格数值模拟方法;娄路亮 等;航空材料学报;第21卷(第3期);全文 *

Also Published As

Publication number Publication date
CN113722965A (zh) 2021-11-30

Similar Documents

Publication Publication Date Title
US11429765B2 (en) Meshless method for solid mechanics simulation, electronic device, and storage medium
Oliver et al. Reduced order modeling strategies for computational multiscale fracture
JP7058902B2 (ja) ハイブリッド繊維複合材料の板巻きシェル構造に対する高速協調最適化方法
Bathe et al. The finite element method with overlapping elements–a new paradigm for CAD driven simulations
CN116151084B (zh) 基于结构网格的模拟方法、装置、终端设备及存储介质
CN109684740B (zh) 一种基于混合网格及时间步长的电磁学多尺度计算方法
CN113722965B (zh) 一种基于积分-广义有限差分数值离散算子的断裂模拟方法
Callens et al. Local explicit interval fields for non-stationary uncertainty modelling in finite element models
Daxini et al. Numerical shape optimization based on meshless method and stochastic optimization technique
Jiang et al. Explicit layout optimization of complex rib-reinforced thin-walled structures via computational conformal mapping (CCM)
CN105404751B (zh) 基于热-力-电磁场网络统一的实现方法
Li et al. Analysis of elasto-plastic thin-shell structures using layered plastic modeling and absolute nodal coordinate formulation
Lin et al. A simple non-conforming isogeometric formulation with superior accuracy for free vibration analysis of thin beams and plates
CN112231863B (zh) 太阳翼电池阵基板建模方法、装置、设备及存储介质
US7773814B2 (en) Multi-scale analysis device
CN105160130B (zh) 一种基于三维图像的有限差分法预测材料热导率的方法
CN115630542A (zh) 一种薄壁加筋结构的加筋布局优化方法
Remmers et al. Isogeometric analysis for modelling of failure in advanced composite materials
Ge et al. Blending isogeometric and Lagrangian elements in three-dimensional analysis
CN114297877A (zh) 杆结构超材料结构多工况仿真自动化系统及方法
Vescovini Analysis of variable stiffness panels with complex geometries using R-functions
Ma Application of Building Information Model (BIM) in the Design of Marine Architectural Structures
Cheng et al. Re-analysis method for inversion of block matrix based on change threshold
Georgiou Interactive Structural Analysis and Form-Finding
Leidinger et al. Explicit isogeometric b-rep analysis on trimmed nurbs-based multi-patch cad models in ls-dyna

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