带有冰纹效果的蜡染印花图案仿真方法

摘要

本发明涉及数字图像处理，为提供带有冰纹效果的蜡染印花图案仿真方法，使其在这种古老艺术的保护、记录与传承方面发挥积极作用，本发明采取的技术方案是，带有冰纹效果的蜡染印花图案仿真方法，包括如下步骤：1)起始开裂冰纹的确定：采用随机法确定起始开裂冰纹；2)冰纹扩展过程的模拟：基于断裂力学中根据格里菲斯断裂准则求解二维裂纹断裂的扩展有限元方法模拟冰纹的开裂过程，模拟出起始开裂冰纹扩展的过程，进而获得每一步扩展后冰纹的冰纹尖端端点信息，得到一系列的冰纹尖端端点后，再将其连成连续平滑的流线冰纹；3）冰纹的绘制：产生不同深浅颜色的蜡染冰纹。本发明主要应用于数字图像处理、蜡染场合。

说明书

技术领域

本发明涉及非物质文化遗产保护、特效制作、数字图像处理与数字艺术设计等领域，具 体讲，涉及带有冰纹效果的蜡染印花图案仿真方法。

技术背景

蜡染印花技术是一种古老而独特的防染织布印花技术。在蜡染印花图案的制作过程中， 由于涂有蜡的织布在染液中不断的翻卷浸染中，作为防染剂的蜡层自然开裂，染料随着裂缝 渗透到织布上而形成人工难以描绘的天然纹理，生成了独特的出于自然的冰纹艺术作品。即 使同样的印花图案在浸染之后，因蜡的开裂其带有偶然性和随机性，得到的冰纹也不会相同， 因开裂而形成的冰纹效果是其区别于其它印花技术的最大特征，可以说冰纹是蜡染印花图案 的灵魂，这种带有特殊纹理为蜡染印花织物增添了独特的韵味。本发明根据以上冰纹效果的 蜡染图案的特征，从冰纹产生的物理原理着手，设计了一种带有冰纹效果的蜡染印花图案仿 真方法。

目前，对于蜡染图案的冰纹效果的模拟，主要有以下方法：

一是运用Photoshop、Illustrator等平面设计软件来制作带有冰纹效果的蜡染印花图 案。Lawrence等[1]通过扫描真实的手工蜡染印花花布或是搜集印刷蜡染印花图案的设计稿 来收集带有冰纹效果的数字图像，然后使用Photoshop等平面设计软件直接从收集的图像中 提取冰纹，将其复制到新创作的蜡染印花图案上，进而产生新的带有冰纹的蜡染印花图案。 Michael[2]和Li[3]也通过使用Photoshop软件来设计带有冰纹效果的蜡染印花图案，但他 们主要采用纯手工的方式绘制各种冰纹，通过在软件中在已设计完成的数字蜡染印花图像中 添加冰纹图形最后合成蜡染冰纹印花图案。这种主要依靠手工方式在软件中直接设计图案的 方法虽然操作简单，但重复性的工作量比较大。

二是基于几何特性的模拟方法，比如距离变换算法和基于Voronoi图的方法。Wyvill等 [4]根据蜡染印花图案中的冰纹更容易在在涂蜡区域比较窄的地方发生开裂的特点，基于距离 变换算法[5]，将初始开裂点设定在距离边界最大处，通过计算涂蜡区域内各部分与已有冰纹 的距离信息，以及冰纹总是向离边界最近的方向生成扩展的原理，达到了模拟蜡染印花图案 中冰纹的目的。该方法可以方便快捷地模拟冰纹的大体形态，但是没有考虑冰纹的一些更为 细致的特征，使得生成的冰纹出现这线型，同时不够连续，在冰纹的仿真效果上还不够自然、 逼真。唐颖等[6]发现蜡染印花图案中的冰纹在一般情况下具有网络性和随机性等分布特点， 且Voronoi图的结构特征与蜡染印花图案中冰纹自然龟裂的情况非常符合，提出了一种基于 Voronoi图的蜡染印花图案冰纹效果的模拟，但是由于Voronoi图的结构比较规则，而自然 龟裂的方向具有随机型，这种方法对于冰纹的开裂方向不易控制；因此，Tang等[7]在距离 变换算法的基础上结合Voronoi图来模拟蜡染印花图案的冰纹效果，利用了两种算法各自的 优势，产生了连续的冰纹效果，同时冰纹的开裂方向接近真实蜡染印花中的冰纹分布情况。 但是由于采用Voronoi图构成的冰纹都是内部连续的冰纹，而真实的冰纹效果并不完全是连 续的，它内部存在一定的随机断裂性。

三是基于物理的模拟方法。O'Brien与Hodgins[8]采用基于四面体网格的有限元模型对 脆性材料进行建模，根据断裂力学的断裂准则判断该网格所对应的材料是否破碎，并根据断 裂力学理论计算破碎裂纹的开裂大小和扩展方向，最终得到了真实感较强的像玻璃、墙壁、 瓷碗等材料断裂模拟效果。后来他们[9]在以上有限元模型中加入对塑性材料产生形变的因 素，模拟实现了带有形变效果的塑性材料的破碎断裂现象。Muller和Teschner[10]在上述方 法基础上把将有限元网格模型统一为若干个大小相等的六面体网格，每个立方体网格表面又 包含了一系列的三角网格，根据三角网格划分破碎区域，进而更新立方体网格模型，最后实 现了像建筑物、车子等三维物理的破碎断裂现象。Parker等[11]同样延续了O'Brien的研究 成果，对有限元计算模型进行了简化处理，提出一种新的物理模拟引擎能够比较高效地实现 游戏场景中的脆性物体破碎特效。Bao等[12]同样基于有限元分析的方法，提出将模拟对象 离散成刚体四面体网格模型，在确定破碎裂纹之后，采用虚节点拓扑重组生成新的面片，进 而形成新的四面体网格结构，这种方法不用担心时间步长的约束，且避免了在求解非线性方 程时的时间开销，能够高效地模拟脆性物体的破碎效果。Su等[13]在Bao等人[12]的研究基 础上充分考虑了物理破碎时的动量和能量守恒，提出了碰撞中心点的预处理机制，提高了模 拟的实时性。基于物理的方法对于各种裂纹的模拟真实感较高，但是复杂的网格划分和数学 方程的计算过程带来了巨大的开销。

[1]Lawrence G M.Digital Printing and Traditional Surface Design Techniques.North Carolina  State University,2002.

[2]Michael C K,Jaffrey R,Marlborough N H,et.al.Method of producing a batik type image on  cloth,1995.

[3]Li D,Yan Y N.Design of dig‐art imitating batik pattern.Advanced materials research,2011, 214(2):260‐264

[4]Wyvill B,Overveld K V,Carpendale S.Rendering cracks in Batik.Proceedings of the 3rd  international symposium on Non‐photorealistic animation and rendering,2004.61‐69.

[5]Jain A K.Fundamentals of Digital Image Processing.Prentice‐Hall Englewood Cliffs,1989,3: 342‐430.

[6]唐颖，房宽峻，付少海等，蜡染冰纹的计算机仿真研究，纺织学报，2010，31（3）：128‐132.

[7]Tang Ying,Fang,Kuanjun,Fu,Shaohai et al.An improved algorithm for simulating  wax‐printing patterns.Textile Research Journal,2011,81(14):1510‐1512,

[8]O'Brien J,Hodgins J.Graphical modeling and animation of brittle fracture.Proceedings of  the26th Annual Conference on Computer Graphics and Interactive Techniques,1999. 137‐146.

[8]O'Brien J,Bargteil A,Hodgins J.Graphical modeling and animation of ductile fracture.ACM  Transactions on Graphics,2002,21(3):291‐294.

[10]Muller M,Teschner M.Physically‐based simulation of objects represented by surface  meshes.Proceedings of International Computer Graphics,2004.26‐33

[11]Parker E G,O’Brien J F.Real‐time deformation and fracture in a game environment. Proceedings of the2009ACM SIGGRAPH/Eurographics Symposium on Computer  Animation,2009.165‐175.

[12]Bao Z,Hong J M,Teran J,et al.Fracturing rigid materials.IEEE Transactions on  Visualization and Computer Graphics,2007.370‐378.

[13]Su J,Schroeder C,Fedkiw R.Energy stability and fracture for frame rate rigid body  simulations.Proceedings of the2009ACM SIGGRAPH/Eurographics Symposium on  Computer Animation,2009.155‐164.

发明内容

为克服现有技术的不足，利用计算机模拟技术在很大程度上再现带有冰纹效果的蜡染印 花图案，设计一种带有冰纹效果的蜡染印花图案仿真方法，使其在这种古老艺术的保护、记 录与传承方面发挥积极作用。为此，本发明采取的技术方案是，带有冰纹效果的蜡染印花图 案仿真方法，包括如下步骤：

1)起始开裂冰纹的确定：采用随机法确定起始开裂冰纹，即记录起始冰纹尖端的位置信息， 同时为了得到一些有意义的带有冰纹的蜡染印花图案，主要考虑分布比较规律的冰纹效果；

2)冰纹扩展过程的模拟：基于断裂力学中根据格里菲斯断裂准则求解二维裂纹断裂的扩展 有限元方法模拟冰纹的开裂过程，模拟出起始开裂冰纹扩展的过程，进而获得每一步扩展后 冰纹的冰纹尖端端点信息，得到一系列的冰纹尖端端点后，再将其连成连续平滑的流线冰纹；

3）冰纹的绘制：根据每步扩散后的冰纹尖端端点信息确定冰纹上各点的颜色值，进而产 生不同深浅颜色的蜡染冰纹。

起始开裂冰纹的确定具体为：

首先，将彩色图像转化成二值图像。在该二值图像中，黑色代表背景，白色代表目标区 域；接下来，采用距离变换算法对处理后的二值图像进行处理，得到相应的距离变换图像， 一幅二值图像的距离变换图代表白色目标区域中的点距离背景区域的最小距离图，图像中每 个像素点的距离值L(p)的计算如下：

L(p)＝min(|p-ν|),p∈M,v∈B    （1）

其中，M表示整个图像区域，B代表黑色背景区域，经过图像的距离变换后，把待处理 的二值图像变换为一幅灰度图像，其中每个像素点的灰度值表示该像素点到达背景区域的距 离，规定亮度越大的区域，其所在位置距离非涂蜡区域越远；反之，若亮度越小，其所在位 置距离非涂蜡区域越近。

基于扩展有限元的冰纹扩展过程模拟具体为：

将扩展有限元的方法应用到冰纹产生过程的模拟，即：通过水平集方法确定可扩展区域内 各节点相对于已扩展冰纹的位置，进而计算各节点的自由度，确定冰纹尖端的富集区域，更 新该区域的蜡质材料的刚度矩阵和应力场的变化，由此得到冰纹尖端的位移变化，提取冰纹 尖端的扩展方向和程度，获得此刻的冰纹扩展趋势；根据格里菲斯断裂准则，判断冰纹是否 继续扩展，并更新冰纹的尖端信息，根据每次冰纹发生扩展后的冰纹尖端端点信息的记录， 得到一系列有先后扩展顺序的尖端端点，采用平滑线产生方法将其连接成一条完成的冰纹；

对于起始冰纹的扩展过程，采用断裂力学中基于扩展有限元求解二维裂纹问题的方法进 行仿真；

采用水平集方法描述不连续面的几何位置及其扩展过程，使得不连续面的描述独立于有 限元网格，从而在不连续面扩展时可以不进行复杂的网格重构。

冰纹的绘制具体步骤为：

根据不同时间扩展出来的冰纹段采用渐进的颜色进行绘制，根据每次扩展后的冰纹尖端 端点信息，判断出冰纹的产生时间，进而采用如下公式计算在任意时刻每一位置冰纹的颜色 值：

c(x,t)＝c(x₀,t)·(1-d(x,x₀))，（2）

其中，c(x,t)表示在t时刻x位置处的冰纹的颜色值，x₀表示冰纹起始开裂点的位置， d(x,x₀)代表x与x₀间的距离值，其取值范围为(0,1)，d(x,x₀)＝0表示在t时刻x位置位于起 始冰纹开裂位置；类似地，d(x,x₀)＝1表示在t时刻x位置位于此时冰纹的尖端位置；c(x₀,t)表 示在t时刻冰纹起始开裂点的颜色值，其颜色随时间t值的增大而加深。

扩展有限元法的核心内容是将裂纹等不连续结构嵌入至单元内部，在普通有限元位移模 式的基础上，增加一类描述非连续性现象的附加函数，使得在进行裂纹模拟数值计算时，无 需重新划分网格就可以顺利进行,其位移的一般形式：

$u\left(x\right)=\underset{i\in N}{\Sigma }{N}_i\left(x\right)u_i+\underset{i\in N_{\Gamma }}{\Sigma }{N}_i\left(x\right)H\left(x\right)a_i+\underset{i\in N_{\Lambda }}{\Sigma }{N}_i\left(x\right)\underset{a=1}{\overset{4}{\Sigma }}{\phi }_a\left(x\right)b_i^a,---\left(3\right)$

公式（2）右边第一项为标准有限元的部分，第二项反映被裂纹截断区域的位移间断特 性，第三项反映裂纹尖端区域的部分，其中，Γ表示裂纹表面，Λ表示裂纹尖端，N是模拟 区域的节点集合；N_Γ表示属于裂纹表面节点集合，N_Λ表示属于裂纹尖端的节点集合；N_i为 单元节点i的有限元形函数；u_i为节点i的自由度；H(x)是Heaviside函数,当所考察点在 裂纹的上方时H(x)取1，而在裂纹的下方时H(x)取-1；φ_a(x)是附加非连续函数；a_i和为 节点i的附加自由度，节点i的真实位移是由两部分组成，一部分为标准有限元情况下的位 移，另一部分裂纹表面或裂纹尖端的移动位移；

将公式（2）代入积分方程，可以得到扩展有限元方法的离散形式，即

Kd＝f            （4）

其中，d是节点位移；K是刚度矩阵，由每个网格单元刚度组成；f是外力向量，由每 个网格单元的外力集成得到,e代表单元向量，i表示节点i，即

$K=\left(\begin{array}{ccc}{k}_{\mathrm{ij}}^{\mathrm{uu}}& k_{\mathrm{ij}}^{\mathrm{ua}}& k_{\mathrm{ij}}^{\mathrm{ub}}\\ k_{\mathrm{ij}}^{\mathrm{au}}& k_{\mathrm{ij}}^{\mathrm{aa}}& k_{\mathrm{ij}}^{\mathrm{ab}}\\ k_{\mathrm{ij}}^{\mathrm{bu}}& k_{\mathrm{ij}}^{\mathrm{ba}}& k_{\mathrm{ij}}^{\mathrm{bb}}\end{array}\right),$

$k_{\mathrm{ij}}^{\mathrm{rs}}={\int }_{{\Omega }^k}{\left(B_i^r\right)}^T{\mathrm{DB}}_j^s\mathrm{d\Omega },r,s=u,a,b---\left(5\right)$

$f_i^e={\{f_i^u,f_i^v,f_i^{b1},f_i^{b2},f_i^{b3},f_i^{b4}\}}^T{,f}_i^u={\int }_{{\Omega }^k}{N}_i\mathrm{bd\Omega }+{\int }_{\Gamma }{N}_i\mathrm{Pd\Gamma },$

$f_i^v={\int }_{{\Omega }^k}{N}_i\mathrm{Hbd\Omega }+{\int }_{\Gamma }{N}_i\mathrm{HPd\Gamma },f_i^{b_m}={\int }_{{\Omega }^k}{N}_i{\phi }_a\mathrm{bd\Omega }+{\int }_{\Gamma }{N}_i{\phi }_a\mathrm{Pd\Gamma },m=\mathrm{1,2,3,4},---\left(6\right)$

其中，D为弹性本构矩阵，P为边界应力矢量，Ω^k为裂纹经过的有限元单元，和分别为形函数的变形矩阵，分别为：

$B_i^u=\left(\begin{array}{cc}{N}_{i,x}& 0\\ 0& N_{i,y}\\ N_{i,y}& N_{i,x}\end{array}\right),B_i^a=\left(\begin{array}{cc}{\left(N_{i}H\right)}_{,x}& 0\\ 0& {\left(N_{i}H\right)}_{,y}\\ {\left(N_{i}H\right)}_{,y}& {\left(N_{i}H\right)}_{,x}\end{array}\right),$

$B_i^b=\left(\begin{array}{cccc}{B}_i^{b1}& B_i^{b2}& B_i^{b3}& B_i^{b4}\end{array}\right),B_i^{\mathrm{ba}}=\left(\begin{array}{cc}{\left(N_i{\phi }_a\right)}_{,x}& 0\\ 0& {\left(N_i{\phi }_a\right)}_{,y}\\ {\left(N_i{\phi }_a\right)}_{,y}& {\left(N_i{\phi }_a\right)}_{,x}\end{array}\right),a=\mathrm{1,2,3,4};---\left(7\right)$

N_i,x、N_i,y分别代表N_i的x偏导数及N_i的y偏导数，(N_iH)_,x，(N_iH)_,y分别代表N_iH的x 偏导数及N_iH的y偏导数，H是Heaviside函数.(N_iφ_a)_,x、(N_iφ_a)_,y分别表示N_iφ_a的的x、 y偏导数，a=1，2，3，4；分别表示和四个函数。

水平集方法是一种跟踪界面移动的数值技术,它将界面的变化表示成比界面高一维的水 平集曲线,在扩展有限元方法中可用来确定裂纹界面的位置和跟踪其裂纹的生长，水平集函 数常取下列符号距离函数表示,即

$\phi (x,t)=\underset{x_{\Gamma }\in \Gamma }{\pm \min }\left|\right|x-x_{\Gamma }\left|\right|,---\left(8\right)$

如果x位于Γ(t)所定义的裂纹上方,那么公式（6）右边项前面的符号取正,如果x位于 Γ(t)所定义的裂纹上方,那么右边项前面的符号取负。

本发明具备下列技术效果：

本发明能够用于生成较为逼真的带有冰纹效果的蜡染印花图案，可以作为艺术图案创作、 特效制作的设计工具，对于这种古老艺术形式的保护与传承也将起到积极作用。本发明支持 用户自由设定不同的参数、鼠标及键盘操作来对系统进行操作。本发明运行环境简便，成本 低，支持Windows XP、Windows vista、Windows7环境，2G以上内存、1G以上显存配置 的微机。

附图说明

图1展示了利用本发明实现带有冰纹效果的蜡染印花图案效果。

图2冰纹的轴平行可扩展矩形区域。

图3本发明流程图。

具体实施方式

本发明提出了一种带有冰纹效果的蜡染印花图案仿真方法，主要包括以下步骤：

1)起始开裂冰纹的确定。与图案无关联的冰纹因其开裂位置不受蜡染图案的影响，本发明 采用随机法确定起始开裂冰纹，即记录起始冰纹尖端的位置信息。为了得到一些有意义的带 有冰纹的蜡染印花图案，主要考虑了分布比较规律的冰纹效果。

2)冰纹扩展过程的模拟。基于断裂力学中根据格里菲斯断裂准则求解二维裂纹断裂的扩展 有限元方法模拟冰纹的开裂过程，模拟出起始开裂冰纹扩展的过程，进而获得每一步扩展后 冰纹的冰纹尖端端点信息。得到一系列的冰纹尖端端点后，再将其连成连续平滑的流线冰纹。

3）冰纹的绘制。由于冰纹开裂的先后导致各个冰纹段的染色时间也有所不同，根据每步 扩散后的冰纹尖端端点信息确定冰纹上各点的颜色值，进而产生不同深浅颜色的蜡染冰纹。

本发明的具体技术方案如下：

1）起始开裂冰纹的确定

本发明中待处理的输入图像为未带冰纹效果的蜡染印花图案，这类图案大多是轮廓边界 比较清晰的彩色图像。首先，将彩色图像转化成二值图像。在该二值图像中，黑色代表背景， 白色代表目标区域。接下来，采用距离变换算法对处理后的二值图像进行处理，得到相应的 距离变换图像。一幅二值图像的距离变换图代表白色目标区域中的点距离背景区域的最小距 离图。图像中每个像素点的距离值L(p)的计算如下：

L(p)＝min(|p-ν|),p∈M,v∈B    （1）

其中，M表示整个图像区域，B代表黑色背景区域。经过图像的距离变换后，把待处理 的二值图像变换为一幅灰度图像，其中每个像素点的灰度值表示该像素点到达背景区域的距 离。规定亮度越大的区域，其所在位置距离非涂蜡区域越远；反之，若亮度越小，其所在位 置距离非涂蜡区域越近。

2)基于扩展有限元的冰纹扩展过程的模拟

本发明将扩展有限元的方法应用到冰纹产生过程的模拟，即：通过水平集方法确定可扩展 区域内各节点相对于已扩展冰纹的位置，进而计算各节点的自由度，确定冰纹尖端的富集区 域，更新该区域的蜡质材料的刚度矩阵和应力场的变化，由此得到冰纹尖端的位移变化，提 取冰纹尖端的扩展方向和程度，获得此刻的冰纹扩展趋势；根据格里菲斯断裂准则，判断冰 纹是否继续扩展，并更新冰纹的尖端信息。根据每次冰纹发生扩展后的冰纹尖端端点信息的 记录，我们可以得到一系列有先后扩展顺序的尖端端点，采用平滑线产生方法将其连接成一 条完成的冰纹。

对于起始冰纹的扩展过程，本发明采用断裂力学中基于扩展有限元求解二维裂纹问题的 方法进行仿真。扩展有限元法的核心内容是将裂纹等不连续结构嵌入至单元内部，在普通有 限元位移模式的基础上，增加一类描述非连续性现象的附加函数，使得在进行裂纹模拟数值 计算时，无需重新划分网格就可以顺利进行,其位移的一般形式：

$u\left(x\right)=\underset{i\in N}{\Sigma }{N}_i\left(x\right)u_i+\underset{i\in N_{\Gamma }}{\Sigma }{N}_i\left(x\right)H\left(x\right)a_i+\underset{i\in N_{\Lambda }}{\Sigma }{N}_i\left(x\right)\underset{a=1}{\overset{4}{\Sigma }}{\phi }_a\left(x\right)b_i^a,---\left(3\right)$

公式（2）右边第一项为标准有限元的部分，第二项反映被裂纹截断区域的位移间断特 性，第三项反映裂纹尖端区域的部分。其中，Γ表示裂纹表面，Λ表示裂纹尖端，N是模拟 区域的节点集合；N_Γ表示属于裂纹表面节点集合，N_Λ表示属于裂纹尖端的节点集合；N_i为 单元节点i的有限元形函数；u_i为节点i的自由度；H(x)是Heaviside函数,当所考察点在 裂纹的上方时H(x)取1，而在裂纹的下方时H(x)取-1；φ_a(x)是附加非连续函数；a_i和为 节点i的附加自由度。节点i的真实位移是由两部分组成，一部分为标准有限元情况下的位 移，另一部分裂纹表面或裂纹尖端的移动位移。

将公式（2）代入积分方程，可以得到扩展有限元方法的离散形式，即

Kd＝f    （4）

其中，d是节点位移；K是总体刚度矩阵，由单元刚度组成；f是外力向量，由单 元外力向量集成得到,e代表单元向量，i表示节点i，即

$K=\left(\begin{array}{ccc}{k}_{\mathrm{ij}}^{\mathrm{uu}}& k_{\mathrm{ij}}^{\mathrm{ua}}& k_{\mathrm{ij}}^{\mathrm{ub}}\\ k_{\mathrm{ij}}^{\mathrm{au}}& k_{\mathrm{ij}}^{\mathrm{aa}}& k_{\mathrm{ij}}^{\mathrm{ab}}\\ k_{\mathrm{ij}}^{\mathrm{bu}}& k_{\mathrm{ij}}^{\mathrm{ba}}& k_{\mathrm{ij}}^{\mathrm{bb}}\end{array}\right),$

$k_{\mathrm{ij}}^{\mathrm{rs}}={\int }_{{\Omega }^h}{\left(B_i^r\right)}^T{\mathrm{DB}}_j^s\mathrm{d\Omega },r,s=u,a,b---\left(5\right)$

$f_i^e={\{f_i^u,f_i^v,f_i^{b1},f_i^{b2},f_i^{b3},f_i^{b4}\}}^T,---\left(6\right)$

其中，$f_i^u={\int }_{{\Omega }^k}{N}_i\mathrm{bd\Omega }+{\int }_{\Gamma }{N}_i\mathrm{Pd\Gamma },f_i^v={\int }_{{\Omega }^k}{N}_i\mathrm{Hbd\Omega }+{\int }_{\Gamma }{N}_i\mathrm{GPd\Gamma },$

$f_i^{b_m}={\int }_{{\Omega }^k}{N}_i{\phi }_a\mathrm{bd\Omega }+{\int }_{\Gamma }{N}_i{\phi }_a\mathrm{Pd\Gamma },m=\mathrm{1,2,3,4},$

其中，D为弹性本构矩阵，P为边界应力矢量，Ω^k是裂纹经过的有限元单元，和分别为形函数的变形矩阵，分别为：

$B_i^u=\left(\begin{array}{cc}{N}_{i,x}& 0\\ 0& N_{i,y}\\ N_{i,y}& N_{i,x}\end{array}\right),B_i^a=\left(\begin{array}{cc}{\left(N_{i}H\right)}_{,x}& 0\\ 0& {\left(N_{i}H\right)}_{,y}\\ {\left(N_{i}H\right)}_{,y}& {\left(N_{i}H\right)}_{,x}\end{array}\right),$

$B_i^b=\left(\begin{array}{cccc}{B}_i^{b1}& B_i^{b2}& B_i^{b3}& B_i^{b4}\end{array}\right),B_i^{\mathrm{ba}}=\left(\begin{array}{cc}{\left(N_i{\phi }_a\right)}_{,x}& 0\\ 0& {\left(N_i{\phi }_a\right)}_{,y}\\ {\left(N_i{\phi }_a\right)}_{,y}& {\left(N_i{\phi }_a\right)}_{,x}\end{array}\right),a=\mathrm{1,2,3,4};---\left(7\right)$

N_i,x、N_i,y分别代表N_i的x偏导数及N_i的y偏导数，(N_iH)_,x，(N_iH)_,y分别代表N_iH的 x偏导数及N_iH的y偏导数，H是Heaviside函数.(N_iφ_a)_,x、(N_iφ_a)_,y分别表示N_iφ_a的的 x、y偏导数，a=1，2，3，4；分别表示和四个函数。

另外，本发明采用水平集方法描述不连续面的几何位置及其扩展过程，使得不连续面的 描述独立于有限元网格，从而在不连续面扩展时可以不进行复杂的网格重构。水平集方法是 一种跟踪界面移动的数值技术,它将界面的变化表示成比界面高一维的水平集曲线,在扩展 有限元方法中可用来确定裂纹界面的位置和跟踪其裂纹的生长。水平集函数常取下列符号距 离函数表示,即

$\phi (x,t)=\underset{x_{\Gamma }\in \Gamma }{\pm \min }\left|\right|x-x_{\Gamma }\left|\right|,---\left(8\right)$

如果x位于Γ(t)所定义的裂纹上方,那么公式（6）右边项前面的符号取正,如果x位于 Γ(t)所定义的裂纹上方,那么右边项前面的符号取负。

3）冰纹的绘制

在真实冰纹的产生过程中，由于冰纹的产生是有时间先后的，染色的时间、渗入染料的 多少都有所不同。因此，冰纹的色彩大多呈现出深浅不同，且随着染色时间的持续，颜色不 断加深；为了模拟这一现象，本发明根据不同时间扩展出来的冰纹段采用渐进的颜色进行绘 制。根据每次扩展后的冰纹尖端端点信息，可以判断出冰纹的产生时间，进而采用如下公式 计算在任意时刻每一位置冰纹的颜色值。

c(x,t)＝c(x₀,t)·(1-l(x,x₀))，    （2）

其中，c(x,t)表示在t时刻x位置处的冰纹的颜色值，x₀表示冰纹起始开裂点的位置，l(x,x₀) 代表x与x₀间的距离值（其取值范围为(0,1)），l(x,x₀)＝0表示在t时刻x位置位于起始冰纹 开裂位置（类似地，l(x,x₀)＝1表示在t时刻x位置位于此时冰纹的尖端位置），c(x₀,t)表示在 t时刻冰纹起始开裂点的颜色值，其颜色随时间t值的增大而加深。

下面结合附图和具体实施方式，进一步详细说明本发明。

本发明涉及大量方程的求解，采用如下子区域划分的方式可提高计算效率。

以蜡染图案的距离变换图为基础进行估算子区域的划分——即以冰纹开裂点为中心的轴 平行矩形区域进行划分。由于冰纹大多倾向于垂直于图案轮廓边界，而冰纹与图案轮廓边界 的交点往往是涂蜡区域内各点距离整个轮廓边界最近的点，根据蜡染图案的距离变换图，可 以获得冰纹开裂起始点的基本扩展方向；如图所示2，P点为起始冰纹开裂点，通过比较与P 点相邻的各节点的与边界距离值得到取最小的点P₁，即

L(P₁)＝min(L(p′))，    （9）

其中，p′表示点P的上下左右四个像素点。

由于冰纹一般是沿方向相反的两个方向扩展的，则以P点为对称中心取P₁的对应点P₂， 可以得到其实开裂冰纹开裂段P₁P₂,则以P点为扩展其实点的冰纹P₁P₂的扩展区域为以P为原 点，P₁P₂为平行轴的矩形可扩展区域。通过实验可发现，矩形区域的边长r可由以下方程求 出：

r＝L(p)·t，    （10）

其中，根据多次实验，r为冰纹起始开裂点与图案轮廓边界距离值的1.5～2.5倍，即 t∈(1.5,2.5)时，可包含该处冰纹的可扩展区域。通过在确定的局部轴平行冰纹可扩展区域 内采用在断裂力学中基于扩展有限元求解二维裂纹问题的方法进行仿真，这样可以提高仿真 效率。