CN115203847B - 一种基于mpm的各向异性相场断裂算法的仿真方法 - Google Patents
一种基于mpm的各向异性相场断裂算法的仿真方法 Download PDFInfo
- Publication number
- CN115203847B CN115203847B CN202210837074.5A CN202210837074A CN115203847B CN 115203847 B CN115203847 B CN 115203847B CN 202210837074 A CN202210837074 A CN 202210837074A CN 115203847 B CN115203847 B CN 115203847B
- Authority
- CN
- China
- Prior art keywords
- crack
- phase field
- anisotropic
- fracture
- stress
- 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
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 33
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 14
- 238000000034 method Methods 0.000 title claims description 42
- 230000008878 coupling Effects 0.000 claims abstract description 7
- 238000010168 coupling process Methods 0.000 claims abstract description 7
- 238000005859 coupling reaction Methods 0.000 claims abstract description 7
- 239000000463 material Substances 0.000 claims description 61
- 239000002245 particle Substances 0.000 claims description 21
- 239000011159 matrix material Substances 0.000 claims description 12
- 238000009792 diffusion process Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 4
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 3
- 238000007670 refining Methods 0.000 claims description 2
- 230000015556 catabolic process Effects 0.000 abstract description 9
- 238000006731 degradation reaction Methods 0.000 abstract description 9
- 239000013013 elastic material Substances 0.000 abstract description 2
- 230000000694 effects Effects 0.000 description 10
- 230000003044 adaptive effect Effects 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000005336 cracking Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 230000005489 elastic deformation Effects 0.000 description 3
- 230000002452 interceptive effect Effects 0.000 description 3
- 230000001902 propagating effect Effects 0.000 description 3
- 239000007787 solid Substances 0.000 description 3
- 239000006185 dispersion Substances 0.000 description 2
- 230000000977 initiatory effect Effects 0.000 description 2
- 230000011218 segmentation Effects 0.000 description 2
- 230000007306 turnover Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 101100016034 Nicotiana tabacum APIC gene Proteins 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000011438 discrete method Methods 0.000 description 1
- 238000002224 dissection Methods 0.000 description 1
- 238000001035 drying Methods 0.000 description 1
- 229920001971 elastomer Polymers 0.000 description 1
- 239000000806 elastomer Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 239000000835 fiber Substances 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 239000012634 fragment Substances 0.000 description 1
- 238000013467 fragmentation Methods 0.000 description 1
- 238000006062 fragmentation reaction Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000003362 replicative effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000010008 shearing Methods 0.000 description 1
- 239000002002 slurry Substances 0.000 description 1
- 230000000153 supplemental effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- 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/25—Design optimisation, verification or simulation using particle-based methods
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Geometry (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于MPM的各向异性相场断裂算法的仿真方法，涉及计算机图形学、基于物理仿真的领域，为了解决现有的断裂算法中，裂纹尖端较模糊，裂纹表面不准确，且损伤区域的裂纹随着裂纹演化会越来越宽的问题，本发明一种基于MPM的各向异性相场断裂算法，通过细化断裂区域的网格，使其具有较高的分辨率，能够保留裂纹的细节，并且在一定程度上衰减裂纹的厚度；设计了正交各向异性超弹性本构模型和各向异性相场断裂模型，发展了弹性材料各向异性能量部分和裂纹相场的耦合。该算法的正交各向异性能量容易通过退化函数与断裂耦合，自适应的调整网格的分辨率，能够增强裂纹模拟的细节，更准确的确定裂纹表面。
Description
技术领域
本发明属于计算机图形学、基于物理的仿真领域，涉及一种基于MPM的各向异性相场断裂算法的仿真方法。
背景技术
随着计算机图形学的发展，许多真实世界中的材料和自然现象都能进行逼真的再现，断裂的模拟作为计算机图形学领域最热门的研究课题之一，已经广泛应用在虚拟手术、影视特效、游戏等多个领域。然而由于断裂模拟通常需要处理大规模拓扑结构的变化、复杂的碰撞检测以及在断裂过程中产生的数值不稳定问题，此外还需要跟踪裂纹的演化、分支、合并等，因此，如何快速、真实地模拟断裂效果的多样性和细节仍是计算机图形学的难题。
目前大多数方法基于Griffith[Griffith A A.VI.The phenomena of ruptureand flow in solids[J].Philosophical transactions of the royal society oflondon.Series A,containing papers of a mathematical or physical character,1921,221(582-593):163-198.]的断裂力学上进行模拟，最早出现的基于物理的模拟断裂的方法是质点弹簧模型，文献[Hirota K,Tanoue Y,Kaneko T.Generation of crackpatterns with a physical model[J].The visual computer,1998,3(14):126-137.]、[Hirota K,Tanoue Y,Kaneko T.Simulation of three-dimensional cracks[J].TheVisual Computer,2000,16(7):371-378.]提出使用质点弹簧模型表示物体的应力和断裂，当应力超过材质最大应力时弹簧断开，来模拟泥浆干燥时表面开裂以及物体内部裂纹的扩展。文献[Aoki K,Dong N H,Kaneko T,et al.Physically based simulation of crackson drying 3D solids[C]//Proceedings Computer Graphics International,2004.IEEE,2004:357-364.]结合水分模型模拟了物体在不同湿度下表面及内部裂纹扩展效果。基于质点弹簧模型断裂效果的模拟受模型构型的限制，在突然去掉某些弹簧时会出现明显的伪影。因此，出现了基于网格的方法，通过改变网格的拓扑结构来追踪断裂过程。文献[O'brien J F,Hodgins J K.Graphical modeling and animation of brittlefracture[C]//Proceedings of the 26th annual conference on Computer graphicsand interactive techniques.1999:137-146.]首次使用FEM模拟了各向同性脆性断裂，通过分析物体形变时产生的应力，分割与断裂面相交的单元、修改相邻的单元划分裂纹周围的网格。[O'brien J F,Bargteil A W,Hodgins J K.Graphical modeling and animationof ductile fracture[C]//Proceedings of the 29th annual conference on Computergraphics and interactive techniques.2002:291-294.]扩展了1999年的工作，通过将应变加性分解为塑性部分和弹性部分，成功模拟了各向同性韧性断裂。有限元将物体离散为连续的网格单元，不连续的裂纹通过网格节点的分离及网格的分割来表示，由此，许多研究者针对网格重建和细化方法开展了大量工作。文献[Wicke M,Ritchie D,Klingner B M,etal.Dynamic local remeshing for elastoplastic simulation[J].ACM Transactionson graphics(TOG),2010,29(4):1-11.]提出动态网格算法，能够局部细化和粗化网格，以尽可能少的改变四面体网格，从而提高模拟精度。文献[Busaryev O,Dey T K,WangH.Adaptive fracture simulation of multi-layered thin plates[J].ACMTransactions on Graphics(TOG),2013,32(4):1-6.]提出了一种分层模型来模拟薄壳的断裂，并且利用基于Delaunay的三角剖分动态重建网格，增强了裂纹的细节。文献[PfaffT,Narain R,De Joya J M,et al.Adaptive tearing and cracking of thin sheets[J].ACM Transactions on Graphics(TOG),2014,33(4):1-9.]提出了自适应网格划分模拟了多种材质薄壳碎裂效果。文献[Chen Z,Yao M,Feng R,et al.Physics-inspired adaptivefracture refinement[J].ACM Transactions on Graphics(TOG),2014,33(4):1-7.]基于梯度下降流自适应地细化裂纹区域曲面网格。尽管如此，有限元方法依然具有较强的网格依赖性以及大量的网格分割过程，因此，文献[Molino N,Bao Z,Fedkiw R.A virtual nodealgorithm for changing mesh topology during simulation[J].ACM Transactions onGraphics(TOG),2004,23(3):385-392.]提出虚拟节点方法，通过复制体网格或曲面网格单元，并且为每个复制的单元分配部分材料，避免了网格的剖分和重建的过程。文献[Wang Y,Jiang C,Schroeder C A,et al.An Adaptive Virtual Node Algorithm with RobustMesh Cutting[C]//Symposium on Computer Animation.2014:77-85.]扩展了虚拟节点工作，使得切割能够通过单元的顶点，并且能够自适应地分配材料。文献[Koschier D,BenderJ,Thuerey N.Robust eXtended finite elements for complex cutting ofdeformables[J].ACM Transactions on Graphics(TOG),2017,36(4):1-13.]提出了一种基于XFEM方法和全隐式时间积分的切割算法，该算法不需要复杂的网格拓扑结构的变化，避免了网格分辨率对模拟效果的影响，且在大时间步长下也能产生稳定的切割效果。
相比于有限元的体积离散，边界元方法更关注材料表面的离散。文献[Hahn D,Wojtan C.High-resolution brittle fracture simulation with boundary elements[J].ACM Transactions on Graphics(TOG),2015,34(4):1-12.]首次使用边界元方法模拟了脆性断裂，使用基于应力和应力强度的准则分开处理裂纹萌生和扩展，并且在裂纹区域使用拉格朗日精细的网格、在弹性形变区域使用粗糙网格，以增强断裂面的细节。文献[ZhuY,Bridson R,Greif C.Simulating rigid body fracture with surface meshes[J].ACMTrans.Graph.,2015,34(4):150:1-150:11.]通过三角网格来表示模拟脆性材料的断裂过程，同时，通过改进后的网格表面跟踪技术来追踪材料裂纹的扩展，可以获得细尺度的裂缝效果。文献[Hahn D,Wojtan C.Fast approximations for boundary element basedbrittle fracture simulation[J].ACM Transactions on Graphics(TOG),2016,35(4):1-11.]扩展了Hahn D的工作，通过快速估计应力强度加速裂纹扩展，设计了边界元网格、隐式网格以及碰撞网格等不同分辨率的网格供用户自由选择，此外，应用Tikhonov正则化器，解决了Neumann边值问题。
MPM方法在形变体和流体模拟中表现出了较强的优越性，MPM将材料离散为粒子，通过在欧拉和拉格朗日空间上分别移动粒子和计算，能够自动处理材料拓扑结构的变化。文献[Wang S,Ding M,Gast T F,et al.Simulation and visualization of ductilefracture with the material point method[J].Proceedings of the ACM on ComputerGraphics and Interactive Techniques,2019,2(2):1-20.]提出一种改进的预评分策略，用四面体单元代替传统MPM的粒子，显式构建断裂表面。
除Griffith断裂力学外，研究者也开始针对连续损伤力学(CDM)展开探索，CDM通过将尖锐裂纹近似为扩散裂纹带从而将裂纹视为连续体，引入一些损伤变量来表示裂纹的演化，其中两种常用的方法PD和PFF也相继应用到图形学中。文献[Chen W,Zhu F,Zhao J,et al.Peridynamics-Based Fracture Animation for Elastoplastic Solids[C]//Computer Graphics Forum.2018,37(1):112-124.]利用基于近场动力学的积分形式代替基于线性断裂力学的偏微分形式对裂纹进行建模，避免了在不连续裂纹处求导。该方法通过积分形式获得材料点周围其他材料点对其产生的作用，能够容易地模拟裂纹的产生和任意方向的扩展。文献[Wolper J,Fang Y,Li M,et al.CD-MPM:continuum damage materialpoint methods for dynamic fracture animation[J].ACM Transactions on Graphics(TOG),2019,38(4):1-15.]提出将CDM与MPM结合来模拟各向同性韧性断裂。通过用连续相场变量表示材料的损伤状态，将裂纹扩展转化为能量变分最小化问题，此外，通过逐渐降低应力对弹性模型进行退化，自然地将相场与弹性能量进行了耦合，能够模拟多种复杂的裂纹效果。文献[Wolper J,Chen Y,Li M,etal.AnisoMPM:Animating anisotropic damagemechanics:Supplemental document[J].ACM Trans.Graph,2020,39(4):37:1–37:16.]扩展了CD-MPM的工作，基于非局部CDM几何方法模拟了各向异性动态断裂和韧性断裂，通过结构张量来构建各向异性相场，并且提出硬约束使嵌入方向不可扩展，从而能够模拟坚硬的材料，在极端变形下具有鲁棒性。文献[Cao W,Lyu L,Ren X,et al.Fracture PatternsDesign for Anisotropic Models with the Material Point Method[C]//ComputerGraphics Forum.2020,39(7):93-104.]提出帧场的概念来控制各向异性的方向，丰富了裂纹的灵活性，此外，重新设计了各向异性弹性和各向异性相场，并对各向异性相场进行了退化。文献[Fan L,Chitalu F M,Komura T.Simulating brittle fracture with materialpoints[J].ACM Transactions on Graphics(TOG),2022.]基于Voronoi细分从损伤材料中提取裂纹表面，对材料碎片体积进行建模，模拟了脆性断裂。
虽然基于MPM的弹性体相场断裂能够自动处理裂纹的产生、扩展及复杂的裂纹分支等，但仍存在下列问题：
(1)裂纹表面没有显示建模，因此难以很好的表现裂纹细节问题；
(2)由于裂纹被具有一定宽度的裂纹平滑表示，因此裂纹尖端较模糊，裂纹表面不准确，且损伤区域的裂纹随着裂纹演化会越来越宽。
发明内容
为了克服上述基于MPM相场断裂模拟存在的问题，本发明提出了一种基于MPM的各向异性相场断裂算法的仿真方法，其目的在于增强裂纹细节，使得模拟的裂纹更接近准确的裂纹表面，首先通过细化即将断裂区域的网格，使其具有较高的分辨率，能够保留裂纹的细节，并且在一定程度上衰减裂纹的厚度；其次，设计了正交各向异性超弹性本构模型，该能量扩展了2019年[Kim T,De Goes F,Iben H.Anisotropic elasticity for inversion-safety and element rehabilitation[J].ACM Transactions on Graphics(TOG),2019,38(4):1-15.]横向各向异性超弹性本构模型，该能量具有翻转安全特性；最后，设计了基于各向异性裂纹密度的各向异性相场断裂模型，且发展了弹性材料各向异性能量部分和裂纹的耦合。为了实现其目的，本发明所采用的技术方案是：
一种基于MPM的各向异性相场断裂算法的仿真方法，其包括以下步骤：
步骤1、初始化参数；
步骤2、当相场超过设定的阈值时，通过窄带对周围网格进行细化；
步骤3、运算相场的梯度和拉普拉斯算子及海森矩阵，并更新相场；
步骤4、将材料点的质量、动量传输到网格，并运算材料点的应力；
步骤5、在网格上更新网格速度，并设置边界条件。
本方法的进一步改进在于：所述步骤1中参数包括材料方向、杨氏模量、泊松比、粒子位置、速度、形变梯度、仿射速度场、体积、密度、损伤、体积比、欧拉网格大小、网格速度、质量、拉普拉斯算子、海森矩阵、网格损伤、时间步长、边界范围和断裂应力阈值。
本方法的进一步改进在于：所述步骤2中将相场值超过0.5的粒子区域，通过插值函数传输到精细网格，其中，插值函数模型为
本方法的进一步改进在于：所述步骤3中基于各向异性裂纹表面密度函数相场模型，实现裂纹扩展呈各向异性，则裂纹表面应用以下函数来表示：
其中，为单位体积裂纹表面密度泛函，通过引入对称二阶结构张量A构建各向异性裂纹表面密度从而将各向异性添加到损伤模型中，/>为裂纹相场的梯度，b1为横向同向裂纹方向，b2为正交各向异性裂纹方向，βi控制材料裂纹方向的强度，l0为正则化扩散裂纹的长度比例，d为相场，F为材料形变梯度。
本方法的进一步改进在于：所述步骤4中根据正交各向异性弹性能量密度与相场耦合更新材料点的应力，正交各向异性弹性能量密度用以下函数来表示：
Ψ(F)＝Ψiso(F)+Ψortho(F)
其中Ψiso(F)和Ψortho(F)分别表示为各项同性部分能量密度和各向异性部分能量密度。
其中F为材料的形变梯度，μ为剪切模量、k为体积模量，C＝FTF为右柯西应力张量，IC＝tr(C)表示C的迹，J＝detF表示F的行列式。I5为各向异性不变量，A＝aaT，a为纤维方向方向量，||·||2表示欧几里得范数。I4＝aTSa，S为形变梯度F极分解F＝RS中的S矩阵。
正交各向异性弹性能量密度退化为：
其中g(d)＝(1-d)2(1-r)+r(d∈[0,1],g(d)∈[0,1])为退化函数。
则应力的计算表示为：
由于采用了上述技术方案，本发明取得的技术进步是：
(1)设计的正交各向异性能量具有翻转安全、体积保持的特性，且容易通过退化函数与断裂耦合；不仅分析断裂对各向同性弹性的退化，也发展了对各向异性弹性部分的退化；
(2)通过窄带，自适应的调整裂纹区域和超弹性区域网格的分辨率，能够增强裂纹模拟的细节，且使用精细的裂纹能在一定程度上减轻裂纹厚度，更准确的确定裂纹表面。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是带有窄带的欧拉网格；
图2是仿真步骤流程图；
图3是扩散裂纹与l0的关系图；
图4为各向同性弹性和横向、纵向各向异性相场；
图5为各向同性弹性和各向同性相场；
图6为各向异性弹性和各向异性相场。
具体实施方式
为实现发明目的，结合技术方案通过具体实施例和附图，对本发明做进一步说明。物质点法是一种混合拉格朗日和欧拉离散方法，将无网格拉格朗日粒子和背景欧拉网格结合，为方便计算，以下变量中下标为p表示拉格朗日粒子的变量，下标为i表示欧拉网格的变量。
仿真步骤包括如下内容，仿真步骤流程图附图3所示：
(1)初始化拉格朗日粒子和欧拉网格相关参数。
(2)当相场d＞0.5时，使用窄带细化周围网格并标记，采用窄带放射粒子法，当相场变量超过阈值时，细化该粒子及其周围粒子的欧拉网格，没有超过阈值的粒子及弹性区域的粒子依旧使用粗糙网格进行模拟，如图1所示。
(3)计算相场的梯度和拉普拉斯算子及海森矩阵，并更新相场。
受Wolper J等人AnisoMPM和Wei Cao等人的启发，引入损伤变量d∈[0,1]来表示裂纹相场，其中d＝1表示完全损伤的材料，d＝0表示没有损伤的材料。利用扩散裂纹来近似尖锐裂纹，本方案采用空间正则化裂纹表面来进行扩散裂纹建模，裂纹表面密度函数如式(1)
其中l0为正则化扩散裂纹的长度尺度参数，为单位体积裂纹表面密度泛函，通过引入对称二阶结构张量A构建各向异性裂纹表面密度从而将各向异性添加到损伤模型中，/>及A的表示如式(2)(3)尖锐裂纹及扩散裂纹随l0的变化如图2所示。
其中▽d为裂纹相场的梯度，b1为横向同向裂纹方向，b2为正交各向异性裂纹方向，βi控制材料裂纹方向的强度，l0为正则化扩散裂纹的长度比例，d为相场，F为材料形变梯度。
相场的演化与相场、相场的梯度以及局部裂纹驱动力相关，局部裂纹驱动力取决于整个材料。根据扩散裂纹拓扑的最小化原则，局部损伤演化方程如式(4)，
其中，η为控制裂纹演化速度的常数，为局部裂纹驱动力，Dc为几何裂纹抵抗力。相场的演化是裂纹驱动力和几何抵抗力之间的平衡，当驱动力大于抵抗力时，会加剧损伤，反之，损伤没有变化。Dc为裂纹表面密度的变分导数，如公式(5)。
其中表示拉普拉斯算子，:表示张量的缩并，/>表示梯度，/>为相场的海森矩阵。
裂纹驱动力与裂纹驱动状态函数的历史有关，因此，将定义为模拟期间每个时刻裂纹驱动力/>的最大值。材料失效有多种准则，如基于材料有阈值应变、无阈值应变、有阈值应力、无阈值应力、主应力准则失效等，本方案采用的是有阈值应力失效准则，/>和/>的定义如式(6)(7)。
其中<>的计算形式为<x>:＝(x+|x|)/2，ζ控制裂纹驱动力的斜率，σc为临界应力阈值，该因子决定了材料断裂前能承受的最大应力。
因此，最终的局部裂纹相场演化公式为
为进一步得到下一个状态的相场，设当前状态为下一步为n+1，两步之间的时间间隔为Δt，当前状态相场为式(9)，
为避免d>1，因此，最得到终的下一个状态的相场更新为
为了方便计算相场d的拉普拉斯算子和海森矩阵，将相场传输到网格上进行计算。首先，通过二次B-spline样条Ni(xp)插值将粒子的相场传输到网格，则
其中，wip＝Ni(xp)，为粒子和节点之间的插值的权重。
利用中心差分计算d的拉普拉斯算子和海森矩阵，并将网格上相场的拉普拉斯算子和海森矩阵传输回粒子。
d_laplacep＝∑iwipd_laplacei (15)
d_hessianp＝∑iwipd_hessiani (16)
根据公式(12)更新粒子的相场。
(4)将材料点的质量、动量传输到网格，并计算材料点的应力，具体步骤如下：
a.采用二次B-Spline插值样条将材料的质量传输到网格节点。
b.发展了弹性各向异性部分的退化。
在相场裂纹与超弹性耦合中，采用单调退化函数g(d)来降低超弹性形变中拉伸、剪切的贡献，而压缩过程中没有退化，从而使材料沿裂纹区域分离。本方案采用常用的二次退化函数耦合弹性能量密度和相场裂纹，如公式(21)。
g(d)＝(1-d)2(1-r)+r (d∈[0,1],g(d)∈[0,1]) (21)
其中r＜＜1为残差应力，以确保即使在完全损伤的区域，也有较小的应力，避免形变梯度可能无限增长。
c.计算第一皮奥拉-基尔霍夫应力。
将模拟的材料看作是域Ω上的连续介质，定义材料空间Ω0和形变空间Ωt的超弹性体，形变映射为x＝φ(X,t)∈Ωt，其中x和X分别为世界坐标和材料坐标，形变梯度其行列式J＝detF表示无限小材料在Ωt和Ω0中的体积比，当J＞1表示材料被拉伸，J<1表示材料被压缩。超弹性体的本构模型由能量密度函数Ψ(F)表示，该能量密度测量了材料点X附近无限小区域dv内的每单位未变形体积应变能。
本发明采用2019年[Wolper J,Fang Y,Li M,et al.CD-MPM:continuum damagematerial point methods for dynamic fracture animation[J].ACM Transactions onGraphics(TOG),2019,38(4):1-15.]的超弹性各向同性本构模型，如式(22)。
其中，Ψiso(F)为各向同性能量密度函数，μ为剪切模量、k为体积模量，C＝FTF为右柯西应力张量，IC＝tr(C)表示C的迹。
Kim T等人提出了一种翻转安全的大变形横向各向同性超弹性本构模型ΨAA，如式(23)。
其中I5为各向异性不变量，A＝ααT，||·||2表示欧几里得范数。
为了解决模拟中奇点(即I5＝0)带来的不稳定状态，引入低阶各向异性不变量I4，并且通过符号函数S(I4)保留奇点处I4符号信息。其定义形式如式(24)(25)。
I4＝αTSα (24)
其中，S为形变梯度F极分解F＝RS中的S矩阵。
在本发明中将Kim T横向各向异性超弹性扩展到正交各向异性超弹性，见式(26)(27)。
其中ai、αi(i＝1,2,3)分别为材料方向强度和对应的材料方向，I4＝αi TSαi，S存在于F的极分解中F＝RS。
本发明最终使用的各向异性超弹性本构模型如式(28)。
Ψ(F)＝Ψiso(F)+Ψortho(F) (28)
因此，各向同性弹性形变部分与相场裂纹的耦合为式(29)
此外，本发明发展了各向异性弹性形变部分与相场裂纹的耦合，如果形变材料在材料方向ai上拉伸，则I5≥1，否则I5＜1，因此，各向异性弹性形变的退化为式(30)
第一皮奥拉-基尔霍夫应力如式(31)(32)
d.通过affine速度场计算动量并传输到网格节点。
其中为APIC速度梯度。
(5)更新网格速度并设置边界条件。
(6)网格速度传输回粒子，并更新粒子速度、位置、形变梯度及Cp。
本发明采用了多个方法进行实验，这里展示部分实验结果。图4为各向同性弹性和横向、纵向各向异性相场，图5为各向同性弹性和各向同性相场，图6为各向异性弹性和各向异性相场。
以上所述的实施例仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明装置权利要求书确定的保护范围内。
Claims (3)
1.一种基于MPM的各向异性相场断裂算法的仿真方法，其特征在于，其包括以下步骤：
步骤1、初始化参数，参数包括材料方向、杨氏模量、泊松比、粒子位置、速度、形变梯度、仿射速度场、体积、质量、密度、损伤、体积比、欧拉网格大小、网格速度、拉普拉斯算子、海森矩阵、网格损伤、时间步长、边界范围和断裂应力阈值；
步骤2、当相场超过设定的阈值时，通过窄带对周围网格进行细化；
步骤3、运算相场的梯度和拉普拉斯算子及海森矩阵，并更新相场；
步骤4、将材料点的质量、动量传输到网格，并运算材料点的应力；
所述步骤4中根据正交各向异性弹性能量密度与相场耦合更新材料点的应力，正交各向异性弹性能量密度用以下函数来表示：
Ψ(F)＝Ψiso(F)+Ψortho(F)
其中Ψiso(F)和Ψortho(F)分别表示为各项同性部分能量密度和各向异性部分能量密度；
其中F为材料的形变梯度，μ为剪切模量、k为体积模量，C＝FTF为右柯西应力张量，IC＝tr(C)表示C的迹，J＝detF表示F的行列式；I5为各向异性不变量，I4为低阶各向异性不变量，ai、αi(i＝1,2,3)分别为材料方向强度和对应的材料方向，I4＝αi TSαi，S为形变梯度F极分解F＝RS中的S矩阵，S(I4)为符号函数；
正交各向异性弹性能量密度退化为：
其中g(d)＝(1-d)2(1-r)+r(d∈[0,1],g(d)∈[0,1])为退化函数，r＜＜1为残差应力；
则应力的计算表示为：
步骤5、在网格上更新网格速度，并设置边界条件。
2.根据权利要求1所述的一种基于MPM的各向异性相场断裂算法的仿真方法，其特征在于，所述步骤2中将相场值超过0.5的粒子区域，通过插值函数传输到精细网格，其中，插值函数模型为
3.根据权利要求1所述的一种基于MPM的各向异性相场断裂算法的仿真方法，其特征在于，所述步骤3中基于各向异性裂纹表面密度函数相场模型，实现裂纹扩展呈各向异性，则裂纹表面应用以下函数来表示：
其中，为单位体积裂纹表面密度泛函，通过引入对称二阶结构张量A构建各向异性裂纹表面密度从而将各向异性添加到损伤模型中，/>为裂纹相场的梯度，b1为横向同向裂纹方向，b2为正交各向异性裂纹方向，βi控制材料裂纹方向的强度，l0为正则化扩散裂纹的长度比例，d为相场，F为材料形变梯度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210837074.5A CN115203847B (zh) | 2022-07-15 | 2022-07-15 | 一种基于mpm的各向异性相场断裂算法的仿真方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210837074.5A CN115203847B (zh) | 2022-07-15 | 2022-07-15 | 一种基于mpm的各向异性相场断裂算法的仿真方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115203847A CN115203847A (zh) | 2022-10-18 |
| CN115203847B true CN115203847B (zh) | 2024-05-28 |
Family
ID=83582355
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210837074.5A Active CN115203847B (zh) | 2022-07-15 | 2022-07-15 | 一种基于mpm的各向异性相场断裂算法的仿真方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115203847B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115618691B (zh) * | 2022-11-10 | 2024-01-26 | 东南大学 | 基于纤维增强复合材料各向异性损伤破裂的相场分析方法 |
| CN115714024B (zh) * | 2022-11-22 | 2023-11-21 | 东南大学 | 组织液-纤维环流固耦合的椎间盘软组织损伤演变预测方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104574472A (zh) * | 2014-12-31 | 2015-04-29 | 北京大学 | 基于嵌入网格的固体碎裂模拟和动画方法 |
| CN113360992A (zh) * | 2021-06-29 | 2021-09-07 | 大连理工大学 | 岩土结构大变形断裂分析的相场物质点方法 |
| CN114186456A (zh) * | 2021-12-02 | 2022-03-15 | 大连理工大学 | 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2018213607A1 (en) * | 2017-05-17 | 2018-11-22 | The Regents Of The University Of California | Computerized rendering of objects having anisotropic elastoplasticity for codimensional frictional contact |
-
2022
- 2022-07-15 CN CN202210837074.5A patent/CN115203847B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104574472A (zh) * | 2014-12-31 | 2015-04-29 | 北京大学 | 基于嵌入网格的固体碎裂模拟和动画方法 |
| CN113360992A (zh) * | 2021-06-29 | 2021-09-07 | 大连理工大学 | 岩土结构大变形断裂分析的相场物质点方法 |
| CN114186456A (zh) * | 2021-12-02 | 2022-03-15 | 大连理工大学 | 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法 |
Non-Patent Citations (4)
| Title |
|---|
| AnisoMPM：Animating Anisotropic Damage Mechanics Supplemental Document;Wolper J.et al;;ACM Trans. Graph;39(4);1-17 * |
| 基于增强的物质点法的非均质弹性材料仿真方法研究;赵静;唐勇;李胜;汪国平;;计算机学报(第11期);180-193 * |
| 基于物质点法的延性断裂建模与并行仿真;赵子朋;中国优秀博士论文全文数据库;1-81 * |
| 形变体仿真中材质本构模型的应用与设计综述;赵静;唐勇;李胜;刘学慧;汪国平;;软件学报(第09期);280-301 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115203847A (zh) | 2022-10-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115203847B (zh) | 一种基于mpm的各向异性相场断裂算法的仿真方法 | |
| Wolper et al. | CD-MPM: continuum damage material point methods for dynamic fracture animation | |
| Macklin et al. | Unified particle physics for real-time applications | |
| Xiao et al. | A novel CNN-based Poisson solver for fluid simulation | |
| Bergou et al. | Tracks: toward directable thin shells | |
| Teran et al. | Robust quasistatic finite elements and flesh simulation | |
| Chelloug et al. | Real Objects Understanding Using 3D Haptic Virtual Reality for E-Learning Education. | |
| Muguercia et al. | Fracture modeling in computer graphics | |
| Wu et al. | Interactive High-Resolution Boundary Surfaces for Deformable Bodies with Changing Topology. | |
| Liu et al. | Dehydration of core/shell fruits | |
| Kan et al. | An immersed MMALE material point method for FSI problems with structure fracturing | |
| Mandal et al. | Interactive physics-based virtual sculpting with haptic feedback | |
| Mercier-Aubin et al. | Adaptive rigidification of elastic solids | |
| Petit et al. | Tracking fractures of deformable objects in real-time with an RGB-D sensor | |
| Nesme et al. | Animating shapes at arbitrary resolution with non-uniform stiffness | |
| Tu et al. | MPM‐driven dynamic desiccation cracking and curling in unsaturated soils | |
| Li et al. | D‐Cloth: Skinning‐based Cloth Dynamic Prediction with a Three‐stage Network | |
| Li et al. | Novel adaptive SPH with geometric subdivision for brittle fracture animation of anisotropic materials | |
| Huang et al. | Geometrically based potential energy for simulating deformable objects | |
| Achar et al. | A Comparative Study of Garment Draping Techniques | |
| Strauss et al. | A semi-automatic surface reconstruction framework based on T-surfaces and isosurface extraction methods | |
| Basloom | A Survey On Physical Methods For Deformation Modeling | |
| Patkar et al. | A new sharp-crease bending element for folding and wrinkling surfaces and volumes | |
| Zhang et al. | Integrating Mesh and Meshfree Methods for Physics-Based Fracture and Debris Cloud Simulation. | |
| Nguyen et al. | A method for improving the accuracy of PODI-RBF solutions for the indentation of an elastic body by a rigid indenter |
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 |