CN114970260A - 一种用于模拟复合材料破坏的格构相场方法 - Google Patents

一种用于模拟复合材料破坏的格构相场方法 Download PDF

Info

Publication number
CN114970260A
CN114970260A CN202210554806.XA CN202210554806A CN114970260A CN 114970260 A CN114970260 A CN 114970260A CN 202210554806 A CN202210554806 A CN 202210554806A CN 114970260 A CN114970260 A CN 114970260A
Authority
CN
China
Prior art keywords
phase field
unit
lattice
composite material
calculation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210554806.XA
Other languages
English (en)
Other versions
CN114970260B (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 University WHU
Original Assignee
Wuhan University WHU
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 University WHU filed Critical Wuhan University WHU
Priority to CN202210554806.XA priority Critical patent/CN114970260B/zh
Publication of CN114970260A publication Critical patent/CN114970260A/zh
Application granted granted Critical
Publication of CN114970260B publication Critical patent/CN114970260B/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/26—Composites

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Operations Research (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Computer Graphics (AREA)
  • Databases & Information Systems (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开了一种用于模拟复合材料破坏的格构相场方法，包括步骤：通过CT扫描技术进行试样结构三维重构，根据重构图像进行三维格构建模，对材料属性参数进行离散化等效，建立适用于格构模型的相场控制方程，对三维格构模型赋予等效的材料参数与边界条件，采用高效迭代算法求解控制方程得到应力场和相场分布，根据相场值判断试样的损伤程度，实现复合材料在荷载作用下裂纹扩展过程、裂纹扩展方向与裂纹扩展路径的精准模拟。本发明为复合材料的数值模拟提供了有效途径，并很大程度上解决了传统相场模型求解裂纹过程中的的计算效率问题；大幅度提高了对裂纹扩展过程的模拟精度，进而解决了脆性复合材料破坏模拟过程中的一系列难题。

Description

一种用于模拟复合材料破坏的格构相场方法
技术领域
本发明属于计算力学领域，具体的说，本发明涉及一种用于模拟复合材料破坏的格构相场方法。
背景技术
相场模型作为基于有限元的新兴弥散方法，在模拟脆性材料断裂破坏时具有着众多优势，一方面避免了传统断裂求解方法所面临的裂纹尖端奇异性问题，另一方面极大程度上消除了裂纹路径的网格敏感性。而由于相场模型的控制方程具有高度的非线性，因此在求解的过程中需要大量的时间耗费，尤其是三维问题，往往由于巨大的计算量令人们在对其应用上望而却步，如何提高相场模型的计算效率，一直以来是科学研究者们探讨的问题，另一方面，有限元法等数值方法在空间内的连续离散对复合材料的模拟并不十分友好，材料的高度非均质性会为网格划分带来困难，寻找适用于复合材料的空间离散方法同样是一个科学难题。
相场模型一直致力于对脆性材料断裂过程的精确预测模拟，而其自身的非线性特性造成计算耗费较大，随着分析问题的维度增加，若材料有限元等连续离散方法，计算量会呈指数性增加，并且造成以下问题：(1)求解速度过慢，难以在有限的时间内得到计算结果；(2)占用微机内存大，对电脑配置要求高；(3)连续性离散限制了对复合材料进行数值模拟的能力。
发明内容
本发明主要解决相场模型在模拟裂纹扩展的过程中的计算耗费大以及连续性离散难以较好地对复合材料进行模拟的问题，为岩石、混凝土一类的脆性复合材料的断裂过程提供一种高效模拟方法。
本发明的上述技术问题主要是通过下述的技术方案得以解决的：
一种用于模拟复合材料破坏的格构相场方法，包括以下步骤：
步骤S1.通过CT扫描成像技术，对复合材料试件的内部结构进行扫描，获得能够表现复合材料各组分分布情况的灰度图像；
步骤S2.对灰度图像进行识别和处理，并根据灰度值对复合材料的组分分布进行重构，获得各组成材料的真实分布数据；
步骤S3.根据图像重构得到的数据，将空间连续的试件采用非连续的格构单元结构进行离散；
步骤S4.建立材料属性参数等效方法，消除采用非连续的格构单元结构离散所带来的误差；
步骤S5.建立适用于格构相场模型的控制方程，控制方程形式如下：
Figure BDA0003654475660000021
其中：第一个公式为控制方程中的弹性方程，第二个公式为控制方程中的裂纹演化方程，Ω表示弹性体的求解域，
Figure BDA0003654475660000022
为Ω的边界，σ为应力张量，u为位移向量，δ为变分符号，δu即表示虚位移，
Figure BDA0003654475660000023
为微分符号，
Figure BDA0003654475660000024
s∈[0,1]为相场变量，相场变量的值可以表现材料的破坏程度s＝0代表材料完好无损，s＝1代表材料完全破坏，
Figure BDA0003654475660000025
和
Figure BDA0003654475660000026
分别为体力和力边界条件，q和Q分别为裂纹演化方程的扩散项和源项，位移向量u和相场变量s即为计算过程中的待求解未知量，位移向量u表征物体变形，相场变量s表征物体破坏程度；
步骤S6.对三维格构模型赋予等效后的材料参数与边界条件，该方法需要的复合材料各组分材料参数有：(1)弹性模量，(2)泊松比，(3)密度，(4)单轴抗拉强度，(5)能量释放率；
步骤S7.采用适用于格构相场方法的高效迭代算法求解控制方程得到位移场场和相场分布，当求得的结果满足设定的收敛条件，则认定该次求解是精确的，开始进行下一荷载步的计算；
步骤S8.根据相场值判断试样的损伤程度，如某位置的相场值达到1，则认定此位移已经发生断裂破坏；
步骤S9.直至施加的边界条件已经全部加载完毕，或者材料发生整体断裂失去承受荷载的能力，计算结束。
进一步地，所述步骤S2中获得各组成材料的真实分布数据包括砂石、水泥砂浆、界面过渡区以及孔洞的分布位置和形状轮廓。
进一步地，所述步骤S3中采用一维桁架单元组成体心格构模块，并以体心格构模块作为基本离散单位建立三维离散体系，进而对连续体进行离散，一个体心格构模块由9个节点和20个桁架单元组成，体心格构模块中的桁架单元可以根据长度分为两类，我们称其为横向单元和斜向单元，它们的长度计算方法为：
Figure BDA0003654475660000031
其中：L为体心格构模块的边长，Ll和Ld分别为横向单元和斜向单元的长度。
进一步地，所述步骤S4中需要进行等效的参数包括桁架单元的横截面积、材料的密度以及材料的能量释放率，其计算方法分别为：
Figure BDA0003654475660000032
Figure BDA0003654475660000033
Figure BDA0003654475660000034
φ＝(9+8ζ)/(18+24ζ) (30)
ζ＝9ν/(4-8ν) (31)
Figure BDA0003654475660000035
其中:Al和Ad分别为横向单元和斜向单元的等效横截面积，ρl和ρd分别为横向单元和斜向单元的等效密度，ρ为等效前的密度，Gl和Gd分别为横向单元和斜向单元的等效能量释放率，Gl和Gd的值相等并可用表达符号Geq代替，Gc为等效前而能量释放率，ν为泊松比，φ、ζ及cA为过渡变量，用来导出泊松比和各等效参数之间的关系，没有实际物理意义。
进一步地，所述步骤S5中，对于其中的弹性方程，采用一维桁架单元的位移场来表达三维连续问题的位移场，即采用不连续离散方法来模拟连续介质的力学响应过程，横向单元和斜向单元的应力张量σ和体力
Figure BDA0003654475660000041
用如下方法计算：
Figure BDA0003654475660000042
Figure BDA0003654475660000043
其中：
Figure BDA0003654475660000044
和
Figure BDA0003654475660000045
分别为横向单元和斜向单元的体力表达形式，g为重力加速度，横向单元和斜向单元的应力向量计算方法相同，均用公式(9)计算，EA为弹性模量，
Figure BDA0003654475660000046
为退化前的应力张量，ε为应变张量，ω(s)为退化函数，其表达形式为：
Figure BDA0003654475660000047
其中：l0为相场模型的裂纹弥散化宽度，ft为单轴抗拉强度，c0是一个尺度参数，表示为：
Figure BDA0003654475660000048
其中：α(s)为裂纹几何函数，用如下公式计算：
Figure BDA0003654475660000049
进一步地，所述步骤S5中控制方程中扩散项q和源项Q采用如下公式表达：
Figure BDA00036544756600000410
Figure BDA00036544756600000411
其中：
Figure BDA00036544756600000412
为裂纹驱动力，并且由如下公式计算：
Figure BDA00036544756600000413
其中：
Figure BDA0003654475660000051
为第n个荷载步的历史场，并且服从Rankine准则
Figure BDA0003654475660000052
Macaulay括号定义为<x>＝max(x,0)，T为第n个荷载步所在的时间，
Figure BDA0003654475660000053
为设定的用于判断单元是否开始发生破坏的历史场临界值。
进一步地，建立了一种适用于格构相场方法的高效迭代求解算法，将控制方程即公式(1)改写为：
Figure BDA0003654475660000054
并进一步简写为：
Figure BDA0003654475660000055
其中r:＝(ru,rs)表示控制方程的总残余，
Figure BDA0003654475660000056
表示为待求解未知量向量，进而方程求解的刚度矩阵可以表示为如下位移场和相场的耦合形式：
Figure BDA0003654475660000057
其中：K为总刚度矩阵，其各分量分别表示为：
Figure BDA0003654475660000058
其中：
Figure BDA0003654475660000059
为位移场的形函数矩阵，I为单位向量，
Figure BDA00036544756600000510
为位移场的协调矩阵，
Figure BDA00036544756600000511
为形函数
Figure BDA00036544756600000512
在桁架单元轴向方向的偏导，
Figure BDA00036544756600000513
为相场的形函数，
Figure BDA00036544756600000514
为相场的协调矩阵，实际上Kus和Ksu相对其它两相较小，其数值对计算结果影响微小，为加快计算和收敛速度，设定Kus＝Ksu＝0恒定成立，因此刚度矩阵可以简化为：
Figure BDA00036544756600000515
为了体现桁架单元在三维空间中的方向，刚度矩阵分量Kuu需要进一步乘以方向转换矩阵，得到最终形式的刚度矩阵：
Figure BDA0003654475660000061
其中：
Figure BDA0003654475660000062
为最终形式的刚度矩阵，
Figure BDA0003654475660000063
为方向转换矩阵，其分量用如下公式计算：
Figure BDA0003654475660000064
其中：(x1,y1,z1)和(x2,y2,z2)分别为桁架单元的两个节点的空间坐标，对于横向单元Le＝ll，对于斜向单元Le＝ld，进而通过计算刚度矩阵可以得到每一迭代步的待求解未知量变化值，即
Figure BDA0003654475660000065
当某一迭代步求得的待求解未知量变化值小于设定的收敛阙值时，即认为计算收敛，开始下一步计算。
因此，与现有技术相比，本发明具有如下优点：
1.格构离散结构的采用和高效迭代求解算法的联合使用，能够大幅度提高计算效率，在很大程度上解决了基于连续性离散的传统相场模型求解所面临的的计算耗费大的问题，一方面，用一维单元对三维空间进行离散，增强了刚度矩阵的稀疏性，并减少了单元内插值的复杂度；另一方面，对刚度矩阵进行了优化，提高了迭代收敛效率，经验证，在自由度数量相同的情况下，该方法相较于传统相场模型求解过程较少计算时间约50％左右。
2.模拟复合材料破坏的格构相场方法将相场模型与格构离散结构有机结合，解决了相场模型在模拟三维问题时由于计算耗费难以施展的问题；桁架单元的使用，提高了刚度矩阵的稀疏度，减少了内存占用，降低了问题求解对计算机配置的要求；同时由于非连续离散的特性，使材料的非均匀特性更加容易体现；在保证了计算精度的同时，效率也将有质的飞跃，使相场模型的应用范围也得到大幅度的扩展；本发明实施过程简单，实用性强，适用范围广，甚至可以向其它领域推广，可在ABAQUS、COMSOL、ANSYS等商业有限元软件中构建，也可通过其它开源有限元软件中或自行编程实施。
附图说明
附图1是本发明双线性自适应相场方法实施技术路线示意图；
附图2是本发明基于桁架单元的体心格构模块示意图及其对三维空间的离散效果示意图；
附图3是本发明处理后的混凝土试样CT扫描图像及各组分分布图；
附图4是本发明格构离散建模及求解裂纹路径示意图；
附图5是本发明混凝土试样开裂过程示意图；
附图6是本发明混凝土内部裂纹切片示意图；
附图7是本发明格构离散与有限元离散计算效率对比图。
具体实施方式
下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。
参阅图1，本发明提供了一种用于模拟复合材料破坏的格构相场方法，以混凝土试样的受拉破坏过程模拟为例，试样为尺寸2.5mm×2.5mm×2.5mm的立方体块，本案例将混凝土视为由砂石、水泥砂浆、界面过渡区(ITZ)以及孔洞组成的结构。本实施例将材料基于桁架单元的体心格构模块对混凝土模型进行离散，体心格构模块的示意图如图2所示。
步骤S1.通过CT扫描成像技术，对复合材料试件的内部结构进行扫描，获得能够表现复合材料各组分分布情况的灰度图像；
步骤S2.对灰度图像进行识别和处理，并对复合材料的组分分布进行重构，获得各组成材料的真实分布数据，包括砂石、水泥砂浆、界面过渡区以及孔洞的分布位置和形状轮廓，重构后得到的混凝土内部结构及各组分的分布情况如图3所示，其中图3中(a)为3D CT图像，(b)为孔洞分布图像，(c)为孔界面过渡区分布图像，(d)为水泥砂浆分布图像，(e)为砂石分布图像。
步骤S3.根据图像重构得到的数据，将空间连续的试件采用非连续的格构单元结构进行离散，得到的格构离散模型如图4(a)所示，图4(b)为裂纹路径示意图；
步骤S4.建立材料属性参数等效方法，消除采用非连续的格构单元结构离散所带来的误差；
步骤S5.建立适用于格构相场模型的控制方程，控制方程形式如下：
Figure BDA0003654475660000081
其中：第一个公式为控制方程中的弹性方程，第二个公式为控制方程中的裂纹演化方程，Ω表示弹性体的求解域，
Figure BDA0003654475660000082
为Ω的边界，σ为应力张量，u为位移向量，δ为变分符号，δu即表示虚位移，
Figure BDA0003654475660000083
为微分符号，
Figure BDA0003654475660000084
s∈[0,1]为相场变量，相场变量的值可以表现材料的破坏程度s＝0代表材料完好无损，s＝1代表材料完全破坏，
Figure BDA0003654475660000085
和
Figure BDA0003654475660000086
分别为体力和力边界条件，q和Q分别为裂纹演化方程的扩散项和源项，位移向量u和相场变量s即为计算过程中的待求解未知量，位移向量u表征物体变形，相场变量s表征物体破坏程度；
步骤S6.对三维格构模型赋予等效后的材料属性参数与边界条件，施加的边界条件为试件底面竖直方向位移固定为0，顶面设定一个竖直向上的位移荷载，大小为u*＝0.03mm，分成400个荷载步逐渐施加。材料属性参数等效方法如下公式(2)-(8),该方法需要的复合材料各组分材料参数有：(1)弹性模量，(2)泊松比，(3)密度，(4)单轴抗拉强度，(5)能量释放率。复合材料中由于各组分物体性质的不同，赋予的各材料属性数值也不同。
采用一维桁架单元组成体心格构模块，并以体心格构模块作为基本离散单位建立三维离散体系，进而对连续体进行离散，一个体心格构模块由9个节点和20个桁架单元组成，由于此结构的非连续特性，能够使非均质物体，尤其是复合材料的建模更加方便，体心格构模块中的桁架单元可以根据长度分为两类，我们称其为横向单元和斜向单元，它们的长度计算方法为：
Figure BDA0003654475660000091
其中：L为体心格构模块的边长，Ll和Ld分别为横向单元和斜向单元的长度。
桁架单元的横截面积、材料的密度以及材料的能量释放率，其计算方法分别为：
Figure BDA0003654475660000092
Figure BDA0003654475660000093
Figure BDA0003654475660000094
φ＝(9+8ζ)/(18+24ζ) (54)
ζ＝9ν/(4-8ν) (55)
Figure BDA0003654475660000095
其中:Al和Ad分别为横向单元和斜向单元的等效横截面积，ρl和ρd分别为横向单元和斜向单元的等效密度，ρ为等效前的密度，Gl和Gd分别为横向单元和斜向单元的等效能量释放率，Gl和Gd的值相等并可用表达符号Geq代替，Gc为等效前而能量释放率，ν为泊松比，φ、ζ及cA为过渡变量，用来导出泊松比和各等效参数之间的关系，没有实际物理意义。
混凝土各组成成分的参数选取如表1所示。边界条件设定为底面各方向位移固定为0，顶面施加向上的位移荷载，总施加位移为0.03mm，分100个相等的子步施加。
表1混凝土各成分的材料属性
Figure BDA0003654475660000096
Figure BDA0003654475660000101
步骤S7.采用适用于格构相场方法的高效迭代算法求解控制方程得到位移场场和相场分布，当求得的结果满足设定的收敛条件，则认定该次求解是精确的，开始进行下一荷载步的计算；
为得到材料在加载过程中的破坏过程，创建了适用于格构结构离散的相场模型控制方程，对于其中的弹性方程，不同于传统的连续性相场方法，本发明采用一维桁架单元的位移场来表达三维连续问题的位移场，即采用不连续离散方法来模拟连续介质的力学响应过程，横向单元和斜向单元的应力张量σ和体力
Figure BDA0003654475660000102
用如下方法计算：
Figure BDA0003654475660000103
Figure BDA0003654475660000104
其中：
Figure BDA0003654475660000105
和
Figure BDA0003654475660000106
分别为横向单元和斜向单元的体力表达形式，g为重力加速度，横向单元和斜向单元的应力向量计算方法相同，均用公式(9)计算，EA为弹性模量，
Figure BDA0003654475660000107
为退化前的应力张量，ε为应变张量，ω(s)为退化函数，其表达形式为：
Figure BDA0003654475660000108
其中：l0为相场模型的裂纹弥散化宽度，ft为单轴抗拉强度，c0是一个尺度参数，表示为：
Figure BDA0003654475660000109
其中：α(s)为裂纹几何函数，用如下公式计算：
Figure BDA00036544756600001010
不同于传统的连续性相场方法，本发明控制方程中的裂纹演化方程需要用一维单元的破坏过程来反映三维连续物体的破坏过程，其中的扩散项q和源项Q采用如下公式表达：
Figure BDA0003654475660000111
Figure BDA0003654475660000112
其中：
Figure BDA0003654475660000113
为裂纹驱动力，并且由如下公式计算：
Figure BDA0003654475660000114
其中：
Figure BDA0003654475660000115
为第n个荷载步的历史场，并且服从Rankine准则
Figure BDA0003654475660000116
Macaulay括号定义为<x>＝max(x,0)，T为第n个荷载步所在的时间，
Figure BDA0003654475660000117
为设定的用于判断单元是否开始发生破坏的历史场临界值。
建立了一种适用于格构相场方法的高效迭代求解算法，将控制方程即公式(1)改写为：
Figure BDA0003654475660000118
并进一步简写为：
Figure BDA0003654475660000119
其中r:＝(ru,rs)表示控制方程的总残余，
Figure BDA00036544756600001110
表示为待求解未知量向量，进而方程求解的刚度矩阵可以表示为如下位移场和相场的耦合形式：
Figure BDA00036544756600001111
其中：K为总刚度矩阵，其各分量分别表示为：
Figure BDA00036544756600001112
其中：
Figure BDA00036544756600001113
为位移场的形函数矩阵，I为单位向量，
Figure BDA0003654475660000121
为位移场的协调矩阵，
Figure BDA0003654475660000122
为形函数
Figure BDA0003654475660000123
在桁架单元轴向方向的偏导，
Figure BDA0003654475660000124
为相场的形函数，
Figure BDA0003654475660000125
为相场的协调矩阵，实际上Kus和Ksu相对其它两相较小，其数值对计算结果影响微小，为加快计算和收敛速度，本发明设定Kus＝Ksu＝0恒定成立，能够在减少计算量的同时保证计算的精度，因此刚度矩阵可以简化为：
Figure BDA0003654475660000126
为了体现桁架单元在三维空间中的方向，刚度矩阵分量Kuu需要进一步乘以方向转换矩阵，得到最终形式的刚度矩阵：
Figure BDA0003654475660000127
其中：
Figure BDA0003654475660000128
为最终形式的刚度矩阵，
Figure BDA0003654475660000129
为方向转换矩阵，其分量用如下公式计算：
Figure BDA00036544756600001210
其中：(x1,y1,z1)和(x2,y2,z2)分别为桁架单元的两个节点的空间坐标，对于横向单元Le＝ll，对于斜向单元Le＝ld，进而通过计算刚度矩阵可以得到每一迭代步的待求解未知量变化值，即
Figure BDA00036544756600001211
步骤S8.根据相场值判断试样的损伤程度，如某位置的相场值达到1，则认定此位移已经发生断裂破坏；
步骤S9.直至施加的边界条件已经全部加载完毕，或者材料发生整体断裂失去承受荷载的能力，计算结束。计算所得到的最终裂纹路径如图4-图6所示，其中图5展示了混凝土在加载过程中的逐渐开裂过程，(a)为u*＝0.0035mm，(b)为u*＝0.0045mm，(c)为u*＝0.03mm时开裂示意图。如图6所示为混凝土内部裂纹切片示意图。
为了进一步展示本发明的能力，本案例通过改变离散化的自由度多次进行计算，来对比本发明的格构离散(LPFM)和传统有限元连续离散(PFM)的计算效率，如图7所示，结果显示，本发明能够节省计算时间约50％左右。
综上所述，本发明是一种高效、稳定、实施流程清晰且适应性强的计算方法，该发明可为复合材料在的裂纹扩展模拟问题的求解提供理论与技术支撑，促进有限元及相场模型对岩石、混凝土、地聚物等脆性复合材料的准确模拟。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种用于模拟复合材料破坏的格构相场方法，其特征在于，包括以下步骤：
步骤S1.通过CT扫描成像技术，对复合材料试件的内部结构进行扫描，获得能够表现复合材料各组分分布情况的灰度图像；
步骤S2.对灰度图像进行识别和处理，并根据灰度值对复合材料的组分分布进行重构，获得各组成材料的真实分布数据；
步骤S3.根据图像重构得到的数据，将空间连续的试件采用非连续的格构单元结构进行离散；
步骤S4.建立材料属性参数等效方法，消除采用非连续的格构单元结构离散所带来的误差；
步骤S5.建立适用于格构相场模型的控制方程，控制方程形式如下：
Figure FDA0003654475650000011
其中：第一个公式为控制方程中的弹性方程，第二个公式为控制方程中的裂纹演化方程，Ω表示弹性体的求解域，
Figure FDA0003654475650000012
为Ω的边界，σ为应力张量，u为位移向量，δ为变分符号，δu即表示虚位移，
Figure FDA0003654475650000013
为微分符号，
Figure FDA0003654475650000014
s∈[0,1]为相场变量，相场变量的值可以表现材料的破坏程度s＝0代表材料完好无损，s＝1代表材料完全破坏，
Figure FDA0003654475650000015
和
Figure FDA0003654475650000016
分别为体力和力边界条件，q和Q分别为裂纹演化方程的扩散项和源项，位移向量u和相场变量s即为计算过程中的待求解未知量，位移向量u表征物体变形，相场变量s表征物体破坏程度；
步骤S6.对三维格构模型赋予等效后的材料参数与边界条件，该方法需要的复合材料各组分材料参数有：(1)弹性模量，(2)泊松比，(3)密度，(4)单轴抗拉强度，(5)能量释放率；
步骤S7.采用适用于格构相场方法的高效迭代算法求解控制方程得到位移场场和相场分布，当求得的结果满足设定的收敛条件，则认定该次求解是精确的，开始进行下一荷载步的计算；
步骤S8.根据相场值判断试样的损伤程度，如某位置的相场值达到1，则认定此位移已经发生断裂破坏；
步骤S9.直至施加的边界条件已经全部加载完毕，或者材料发生整体断裂失去承受荷载的能力，计算结束。
2.根据权利要求1所述的一种用于模拟复合材料破坏的格构相场方法，其特征性在于：所述步骤S2中获得各组成材料的真实分布数据包括砂石、水泥砂浆、界面过渡区以及孔洞的分布位置和形状轮廓。
3.根据权利要求1所述的一种用于模拟复合材料破坏的格构相场方法，其特征性在于：所述步骤S3中采用一维桁架单元组成体心格构模块，并以体心格构模块作为基本离散单位建立三维离散体系，进而对连续体进行离散，一个体心格构模块由9个节点和20个桁架单元组成，体心格构模块中的桁架单元可以根据长度分为两类，我们称其为横向单元和斜向单元，它们的长度计算方法为：
Figure FDA0003654475650000021
其中：L为体心格构模块的边长，Ll和Ld分别为横向单元和斜向单元的长度。
4.根据权利要求1所述的一种用于模拟复合材料破坏的格构相场方法，其特征性在于：所述步骤S4中需要进行等效的参数包括桁架单元的横截面积、材料的密度以及材料的能量释放率，其计算方法分别为：
Figure FDA0003654475650000022
Figure FDA0003654475650000023
Figure FDA0003654475650000024
φ＝(9+8ζ)/(18+24ζ) (6)
ζ＝9ν/(4-8ν) (7)
Figure FDA0003654475650000031
其中:Al和Ad分别为横向单元和斜向单元的等效横截面积，ρl和ρd分别为横向单元和斜向单元的等效密度，ρ为等效前的密度，Gl和Gd分别为横向单元和斜向单元的等效能量释放率，Gl和Gd的值相等并可用表达符号Geq代替，Gc为等效前而能量释放率，ν为泊松比，φ、ζ及cA为过渡变量，用来导出泊松比和各等效参数之间的关系，没有实际物理意义。
5.根据权利要求1所述的一种用于模拟复合材料破坏的格构相场方法，其特征性在于：所述步骤S5中，对于其中的弹性方程，采用一维桁架单元的位移场来表达三维连续问题的位移场，即采用不连续离散方法来模拟连续介质的力学响应过程，横向单元和斜向单元的应力张量σ和体力
Figure FDA0003654475650000032
用如下方法计算：
Figure FDA0003654475650000033
Figure FDA0003654475650000034
其中：
Figure FDA0003654475650000035
和
Figure FDA0003654475650000036
分别为横向单元和斜向单元的体力表达形式，g为重力加速度，横向单元和斜向单元的应力向量计算方法相同，均用公式(9)计算，EA为弹性模量，
Figure FDA0003654475650000037
为退化前的应力张量，ε为应变张量，ω(s)为退化函数，其表达形式为：
Figure FDA0003654475650000038
其中：l0为相场模型的裂纹弥散化宽度，ft为单轴抗拉强度，c0是一个尺度参数，表示为：
Figure FDA0003654475650000039
其中：α(s)为裂纹几何函数，用如下公式计算：
Figure FDA00036544756500000310
6.根据权利要求1所述的一种用于模拟复合材料破坏的格构相场方法，其特征性在于：所述步骤S5中控制方程中扩散项q和源项Q采用如下公式表达：
Figure FDA0003654475650000041
Figure FDA0003654475650000042
其中：
Figure FDA0003654475650000043
为裂纹驱动力，并且由如下公式计算：
Figure FDA0003654475650000044
其中：
Figure FDA0003654475650000045
为第n个荷载步的历史场，并且服从Rankine准则
Figure FDA0003654475650000046
Macaulay括号定义为<x>＝max(x,0)，T为第n个荷载步所在的时间，
Figure FDA0003654475650000047
为设定的用于判断单元是否开始发生破坏的历史场临界值。
7.根据权利要求1所述的一种用于模拟复合材料破坏的格构相场方法，其特征性在于：建立了一种适用于格构相场方法的高效迭代求解算法，将控制方程即公式(1)改写为：
Figure FDA0003654475650000048
并进一步简写为：
Figure FDA0003654475650000049
其中r:＝(ru,rs)表示控制方程的总残余，
Figure FDA00036544756500000410
表示为待求解未知量向量，进而方程求解的刚度矩阵可以表示为如下位移场和相场的耦合形式：
Figure FDA00036544756500000411
其中：K为总刚度矩阵，其各分量分别表示为：
Figure FDA0003654475650000051
其中：
Figure FDA0003654475650000052
为位移场的形函数矩阵，I为单位向量，
Figure FDA0003654475650000053
为位移场的协调矩阵，
Figure FDA0003654475650000054
为形函数
Figure FDA0003654475650000055
在桁架单元轴向方向的偏导，
Figure FDA0003654475650000056
为相场的形函数，
Figure FDA0003654475650000057
为相场的协调矩阵，实际上Kus和Ksu相对其它两相较小，其数值对计算结果影响微小，为加快计算和收敛速度，设定Kus＝Ksu＝0恒定成立，因此刚度矩阵可以简化为：
Figure FDA0003654475650000058
为了体现桁架单元在三维空间中的方向，刚度矩阵分量Kuu需要进一步乘以方向转换矩阵，得到最终形式的刚度矩阵：
Figure FDA0003654475650000059
其中：
Figure FDA00036544756500000510
为最终形式的刚度矩阵，
Figure FDA00036544756500000511
为方向转换矩阵，其分量用如下公式计算：
Figure FDA00036544756500000512
其中：(x1,y1,z1)和(x2,y2,z2)分别为桁架单元的两个节点的空间坐标，对于横向单元Le＝ll，对于斜向单元Le＝ld，进而通过计算刚度矩阵可以得到每一迭代步的待求解未知量变化值，即
Figure FDA00036544756500000513
当某一迭代步求得的待求解未知量变化值小于设定的收敛阙值时，即认为计算收敛，开始下一步计算。
CN202210554806.XA 2022-05-20 2022-05-20 一种用于模拟复合材料破坏的格构相场方法 Active CN114970260B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210554806.XA CN114970260B (zh) 2022-05-20 2022-05-20 一种用于模拟复合材料破坏的格构相场方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210554806.XA CN114970260B (zh) 2022-05-20 2022-05-20 一种用于模拟复合材料破坏的格构相场方法

Publications (2)

Publication Number Publication Date
CN114970260A true CN114970260A (zh) 2022-08-30
CN114970260B CN114970260B (zh) 2024-04-02

Family

ID=82985886

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210554806.XA Active CN114970260B (zh) 2022-05-20 2022-05-20 一种用于模拟复合材料破坏的格构相场方法

Country Status (1)

Country Link
CN (1) CN114970260B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115544834A (zh) * 2022-09-30 2022-12-30 东南大学 基于相场模型的混凝土材料损伤演变模拟方法
CN115618691A (zh) * 2022-11-10 2023-01-17 东南大学 基于纤维增强复合材料各向异性损伤破裂的相场分析方法
CN117494486A (zh) * 2024-01-03 2024-02-02 南通泰胜蓝岛海洋工程有限公司 一种局部集中载荷作用下的组合梁精细化应力位移分析方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6813592B1 (en) * 1999-06-20 2004-11-02 Mtu Aero Engines Gmbh Method for crack propagation simulation
US20160299046A1 (en) * 2014-09-25 2016-10-13 East China University Of Science And Technology A method of measurement and determination on fracture toughness of structural materials at high temperature
WO2019238451A1 (en) * 2018-06-13 2019-12-19 Danmarks Tekniske Universitet A method and a system for modelling and simulating a fractured geological structure
CN111177901A (zh) * 2019-12-18 2020-05-19 常州工学院 基于Mohr-Coulomb准则的围岩破坏危险性数值模拟评价方法
CN112036060A (zh) * 2020-08-03 2020-12-04 武汉大学 一种用于模拟脆性材料破坏的双线性自适应相场方法
CN112051142A (zh) * 2020-08-03 2020-12-08 武汉大学 一种用于模拟脆性材料不同破坏模式的普适性相场方法
CN113094946A (zh) * 2021-03-23 2021-07-09 武汉大学 一种用于模拟材料开裂的相场模型局部化自适应算法
US20210340857A1 (en) * 2020-04-29 2021-11-04 Saudi Arabian Oil Company Method for automated crack detection and analysis using ultrasound images
US20220075911A1 (en) * 2020-06-11 2022-03-10 Dalian University Of Technology Method for predicting structural failure by strength-criterion-driven peridynamic model

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6813592B1 (en) * 1999-06-20 2004-11-02 Mtu Aero Engines Gmbh Method for crack propagation simulation
US20160299046A1 (en) * 2014-09-25 2016-10-13 East China University Of Science And Technology A method of measurement and determination on fracture toughness of structural materials at high temperature
WO2019238451A1 (en) * 2018-06-13 2019-12-19 Danmarks Tekniske Universitet A method and a system for modelling and simulating a fractured geological structure
CN111177901A (zh) * 2019-12-18 2020-05-19 常州工学院 基于Mohr-Coulomb准则的围岩破坏危险性数值模拟评价方法
US20210340857A1 (en) * 2020-04-29 2021-11-04 Saudi Arabian Oil Company Method for automated crack detection and analysis using ultrasound images
US20220075911A1 (en) * 2020-06-11 2022-03-10 Dalian University Of Technology Method for predicting structural failure by strength-criterion-driven peridynamic model
CN112036060A (zh) * 2020-08-03 2020-12-04 武汉大学 一种用于模拟脆性材料破坏的双线性自适应相场方法
CN112051142A (zh) * 2020-08-03 2020-12-08 武汉大学 一种用于模拟脆性材料不同破坏模式的普适性相场方法
CN113094946A (zh) * 2021-03-23 2021-07-09 武汉大学 一种用于模拟材料开裂的相场模型局部化自适应算法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
QIANG YUE: "A phase-field lattice model (PFLM) for fracture problem", COMPOSITE STRUCTURE, 31 August 2023 (2023-08-31), pages 1 - 17 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115544834A (zh) * 2022-09-30 2022-12-30 东南大学 基于相场模型的混凝土材料损伤演变模拟方法
CN115544834B (zh) * 2022-09-30 2023-11-07 东南大学 基于相场模型的混凝土材料损伤演变模拟方法
CN115618691A (zh) * 2022-11-10 2023-01-17 东南大学 基于纤维增强复合材料各向异性损伤破裂的相场分析方法
CN115618691B (zh) * 2022-11-10 2024-01-26 东南大学 基于纤维增强复合材料各向异性损伤破裂的相场分析方法
CN117494486A (zh) * 2024-01-03 2024-02-02 南通泰胜蓝岛海洋工程有限公司 一种局部集中载荷作用下的组合梁精细化应力位移分析方法
CN117494486B (zh) * 2024-01-03 2024-04-02 南通泰胜蓝岛海洋工程有限公司 一种局部集中载荷作用下的组合梁精细化应力位移分析方法

Also Published As

Publication number Publication date
CN114970260B (zh) 2024-04-02

Similar Documents

Publication Publication Date Title
CN114970260A (zh) 一种用于模拟复合材料破坏的格构相场方法
Wei et al. A study on X-FEM in continuum structural optimization using a level set model
Gao et al. Topology optimization of micro-structured materials featured with the specific mechanical properties
Nguyen et al. Two-and three-dimensional isogeometric cohesive elements for composite delamination analysis
Simone et al. From continuous to discontinuous failure in a gradient-enhanced continuum damage model
Gu et al. Multi-inclusions modeling by adaptive XIGA based on LR B-splines and multiple level sets
Düster et al. Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method
Ullah et al. An adaptive finite element/meshless coupled method based on local maximum entropy shape functions for linear and nonlinear problems
Tanaka et al. Study on crack propagation simulation of surface crack in welded joint structure
Shen et al. Multivariate uncertainty analysis of fracture problems through model order reduction accelerated SBFEM
Zhao et al. Multiscale topology optimization using feature-driven method
Yu et al. Error-controlled adaptive LR B-splines XIGA for assessment of fracture parameters in through-cracked Mindlin-Reissner plates
Evangelista Jr et al. A global–local strategy with the generalized finite element framework for continuum damage models
CN114756934B (zh) 一种三维多尺度超材料结构优化设计方法
March et al. Evaluation of computational homogenization methods for the prediction of mechanical properties of additively manufactured metal parts
Schiftner et al. Statics-sensitive layout of planar quadrilateral meshes
Chiappa et al. Improvement of 2D finite element analysis stress results by radial basis functions and balance equations
Patni et al. Efficient modelling of beam-like structures with general non-prismatic, curved geometry
Oztoprak et al. Two-scale analysis of spaceframes with complex additive manufactured nodes
Krätzig et al. Onbest'shell models–From classical shells, degenerated and multi-layered concepts to 3D
O’Hara et al. A two-scale generalized finite element method for fatigue crack propagation simulations utilizing a fixed, coarse hexahedral mesh
Abdi Evolutionary topology optimization of continuum structures using X-FEM and isovalues of structural performance
Wang et al. Two-dimensional mixed mode crack simulation using the material point method
Qu et al. Development of a fully automatic damage simulation framework for quasi-brittle materials
ISHIDA et al. Topology optimization for maximizing linear buckling load based on level set method

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