CN111914445A - 基于fft法确定非均质岩石裂纹扩展路径的计算方法 - Google Patents

Abstract

本发明公开了一种基于FFT法确定非均质岩石裂纹扩展路径的计算方法，与原有的基于有限元框架的相场方法相比，基于FFT均匀化方法的相场断裂计算方法可以克服有限元框架的不足，减小尺寸效应对计算的影响，并且在保证计算精度的同时，提高了数值计算效率，降低了计算成本；从微观角度出发，通过均匀化理论建立起微观尺度的裂纹开展和宏观尺度的材料破坏之间的联系，能够准确地反映材料断裂的机理，为预测非均质材料断裂提供了可靠依据，此外该基于FFT法确定非均质岩石裂纹扩展路径的计算方法易于程序化，便于通过编程语言加以实现，通过计算机进行断裂预测，大大减少了人工计算量。

Description

基于FFT法确定非均质岩石裂纹扩展路径的计算方法

技术领域

本发明涉及计算力学数值模拟方法与工程断裂力学技术领域，尤其涉及一种基于FFT法确定非均质岩石裂纹扩展路径的计算方法。

背景技术

非均质岩石材料的断裂机理及预测一直是实际工程问题中的研究热点。非均质材料 中存在初始孔隙、夹杂及微裂纹等缺陷，将导致材料受到外部荷载作用时，其微观、宏观裂纹的开展和演化过程变得十分复杂。而传统的线弹性断裂力学不能描述材料的微观结构，无法预测非均质材料中的裂纹开展过程，所以提出一种能够很好地模拟裂纹的形 成和发展的新方法尤为重要。

近年来，在计算力学领域，相场方法开始应用于模拟断裂。相场法通过引入一个耗散的相场变量d∈[0,1]，其中0代表完整的材料，1代表完全的断裂状态即裂纹。与传 统方法相比，相场方法克服了经典断裂力学理论框架需要复杂的裂纹发展判定法则的不 足，通过相场变量显式地追踪裂缝扩展，且不需要额外追踪裂缝的几何形状，在计算裂 缝三维扩展和裂缝分叉等裂缝复杂扩展路径时具有独特优势。然而，目前的相场方法在 数值实现上多是基于传统的有限元框架，数值结果对网格依赖性较大，在对模型中复杂 的微观结构进行空间离散时，需要精细的网格划分，导致计算效率降低，计算成本增高。

与传统有限元框架不同，FFT均匀化方法通过在傅立叶空间中离散 Lippmann-Schwinger周期性方程来求解非均质材料的局部和全局响应。该方法可直接利 用数字图像技术(如SEM图像等)描述材料的微观结构特征，并满足周期性边界条件 来减小尺寸效应。与有限元方法相比，FFT均匀化方法无需对复杂结构进行网格划分和 组装整体刚度矩阵，在保证计算精度的同时，提高了数值计算效率。因此，将FFT均匀 化方法框架引入到相场方法模拟中，可有效提高数值计算效率，降低计算成本。

综上，目前针对多相复合材料和多孔介质材料等非均质材料的相场模拟方法，计算 效率较低，且并未从微观角度出发，未能通过均匀化理论建立起微观尺度的裂纹开展和宏观尺度的材料破坏之间的联系，不能完全准确地反映材料断裂的机理。FFT均匀化方 法框架虽可以克服有限元框架的不足，但并未应用到三维准脆性非均质材料的相场断裂 方法模拟中。因此，目前缺乏一种基于FFT均匀化方法的三维相场断裂科学计算方法。

发明内容

针对以上问题，本发明提出一种基于FFT法确定非均质岩石裂纹扩展路径的计算方 法，以解决现有技术中的相场方法由于未能与FFT均匀化理论相结合，因而不能描述微观尺度的裂纹开展和宏观尺度的材料破坏之间的联系，计算效率低且不能准确地反映材料断裂机理的技术问题。

为实现本发明的目的，提供一种基于FFT法确定非均质岩石裂纹扩展路径的计算方 法，包括如下步骤：

S20，在非均质岩石材料的特征单元体计算网格上施加宏观应变增量，根据参数输入值，初始化数值计算参数，根据初始应变场计算相场初始相场历史变量；

S30，将初始相场历史变量代入相场控制方程，将相场控制方程进行快速傅里叶变换，推导出傅里叶空间中固定点相场迭代格式；基于斯特芬森加速固定点迭代方法，迭 代计算相场变量，直至满足迭代收敛条件，结束迭代并进行傅里叶逆变换，计算真实空 间中初始相场变量；将初始相场变量代入材料本构关系获得计算网格材料点的应力和应 变初始值，根据应力和应变初始值确定应力场；

S40，根据应力场进行应力平衡条件收敛性判断，如果满足收敛条件，进入步骤S50中的宏观应变增量的施加；否则，在FFT计算框架中，采用交错求解方法，先根据应变 初始值，计算Lippmann-Schwinger周期性方程，更新网格材料点应变场；再固定应变场， 计算网格材料点相场历史变量，依据步骤S30迭代计算本时步网格材料点的相场变量及 应力，直至满足迭代应力平衡收敛条件，结束迭代并更新本时步的网格材料点的相场变 量及应力；

S50，根据各网格材料点的应力，确定所有网格材料点的体积平均，计算出特征单元体的宏观应力值，完成本时步的计算，返回执行步骤S20，并以本时步的宏观应力值 作为步骤S20的参数输入值，直至加载到指定宏观应变值，以获得整个加载过程的相场 d的演化云图，确定裂纹的扩展路径。

在一个实施例中，在步骤S20之前，还包括：

S10，基于相场模型的相场变量，根据材料类别，确定材料劣化函数及相场方法计算控制方程。

在一个实施例中，基于相场模型的相场变量，根据材料类别，确定材料劣化函数及相场方法计算控制方程包括：

基于相场模型的相场变量，根据材料类别，确定材料劣化函数，采用材料劣化函数对弹性应变能进行劣化，获得材料的总势能，使总势能包括弹性应变能和断裂能，以依 据能量最小变分原理确定相场方法计算控制方程。

具体地，所述材料劣化函数包括：

其中，d表示相场变量，函数

参数A的值大于等于1，B 的值大于0。

在一个实施例中，所述参数输入值包括上一时步网格材料点的应变ε(t_n,x)和相场变 量d(t_n,x)，以及施加的宏观应变增量ΔE(t_n+1)和时间增量Δt_n+1；

需要初始化的数值计算参数包括：宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以及初 始相场变量d^0,0(t_n+1,x)，并计算初始相场历史变量

在一个实施例中，所述宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以及初始相场变量 d^0,0(t_n+1,x)计算方式分别为：

E(t_n+1)＝E(t_n)+ΔE(t_n+1)，

ε⁰(t_n+1,x)＝ε(t_n,x)+ΔE(t_n+1)，

d^0,0(t_n+1,x)＝d(t_n,x)，

式中，t_n表示上一时间步，t_n+1＝t_n+Δt_n+1表示本时间步，Δt_n+1为时间增量，x表示位置，E(t_n)表示上一时步的宏观应变，E(t_n+1)表示本时步的宏观应变，ΔE(t_n+1)表示本 时步内的宏观应变变化量，ε(t_n,x)表示材料点x在t_n时刻的应变，d(t_n,x)表示材料点x 在t_n时刻的相场变量值；

所述初始相场历史变量包括：

式中，ψ为弹性自由能，

及

分别为各向同性弹性自由能的正、负分量，k是 一微小正值，λ、μ为拉梅常数，εⁱ、n_i分别为应变张量ε的主应变和主应变方向，符 号

N表示维度。

在一个实施例中，所述相场变量的迭代求解公式为：

式中，φ为相场变量d的函数，上标0代表初始FFT应力平衡迭代循环步，上标i 代表第i步斯特芬森加速固定点迭代；

网格材料点的应力计算公式为：

式中，I为单位二阶张量。

在一个实施例中，所述宏观应力计算公式为：

式中，Ω表示计算空间的体积，t_n+1表示本时间步。

上述基于FFT法确定非均质岩石裂纹扩展路径的计算方法，与原有的基于有限元框 架的相场方法相比，基于FFT均匀化方法的相场断裂计算方法可以克服有限元框架的不足，减小尺寸效应对计算的影响，并且在保证计算精度的同时，提高了数值计算效率， 降低了计算成本；从微观角度出发，通过均匀化理论建立起微观尺度的裂纹开展和宏观 尺度的材料破坏之间的联系，能够准确地反映材料断裂的机理，为预测非均质材料断裂 提供了可靠依据，此外该基于FFT法确定非均质岩石裂纹扩展路径的计算方法易于程序 化，便于通过编程语言加以实现，通过计算机进行断裂预测，大大减少了人工计算量。

附图说明

图1是一个实施例的基于FFT法确定非均质岩石裂纹扩展路径的计算方法流程图；

图2是一个实施例中参数A取值对材料劣化函数的影响示意图；

图3是一个实施例中参数B取值对材料劣化函数的影响示意图；

图4是一个实施例中不同粘性参数下包含一个中心球状空隙的单胞宏观轴向应力- 应变关系示意图；

图5为一个实施例中单轴拉伸试验中单胞的裂纹演化规律示意图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请， 并不用于限定本申请。

在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式 地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。

参考图1所示，图1为一个实施例的基于FFT法确定非均质岩石裂纹扩展路径的计算方法流程图，包括如下步骤：

S20，在非均质岩石材料的特征单元体计算网格上施加宏观应变增量，根据参数输入值，初始化数值计算参数，根据初始应变场计算相场初始相场历史变量。

S30，将初始相场历史变量代入相场控制方程，将相场控制方程进行快速傅里叶变换，推导出傅里叶空间中固定点相场迭代格式；基于斯特芬森加速固定点迭代方法，迭 代计算相场变量，直至满足迭代收敛条件，结束迭代并进行傅里叶逆变换，计算真实空 间中初始相场变量；将初始相场变量代入材料本构关系获得计算网格材料点的应力和应 变初始值，根据应力和应变初始值确定应力场。

S40，根据应力场进行应力平衡条件收敛性判断，如果满足收敛条件，进入步骤S50中的宏观应变增量的施加；否则，在FFT计算框架中，采用交错求解方法，先根据应变 初始值，计算Lippmann-Schwinger周期性方程，更新网格材料点应变场；再固定应变场， 计算网格材料点相场历史变量，依据步骤S30迭代计算本时步网格材料点的相场变量及 应力，直至满足迭代应力平衡收敛条件，结束迭代并更新本时步的网格材料点的相场变 量及应力。

S50，根据各网格材料点的应力，确定所有网格材料点的体积平均，计算出特征单元体的宏观应力值，完成本时步的计算，返回执行步骤S20，并以本时步的宏观应力值 作为步骤S20的参数输入值，直至加载到指定宏观应变值，以获得整个加载过程的相场 d的演化云图，确定裂纹的扩展路径。

上述基于FFT法确定非均质岩石裂纹扩展路径的计算方法，与原有的基于有限元框 架的相场方法相比，基于FFT均匀化方法的相场断裂计算方法可以克服有限元框架的不足，减小尺寸效应对计算的影响，并且在保证计算精度的同时，提高了数值计算效率， 降低了计算成本；从微观角度出发，通过均匀化理论建立起微观尺度的裂纹开展和宏观 尺度的材料破坏之间的联系，能够准确地反映材料断裂的机理，为预测非均质材料断裂 提供了可靠依据，此外该基于FFT法确定非均质岩石裂纹扩展路径的计算方法易于程序 化，便于通过编程语言加以实现，通过计算机进行断裂预测，大大减少了人工计算量。

在一个实施例中，在步骤S20之前，还包括：

S10，基于相场模型的相场变量，根据材料类别，确定材料劣化函数及相场方法计算控制方程。

在一个实施例中，基于相场模型的相场变量，根据材料类别，确定材料劣化函数及相场方法计算控制方程包括：

基于相场模型的相场变量，根据材料类别，确定材料劣化函数，采用材料劣化函数对弹性应变能进行劣化，获得材料的总势能，使总势能包括弹性应变能和断裂能，以依 据能量最小变分原理确定相场方法计算控制方程。

本实施例所选择的双参数材料劣化函数不同于传统的相场的二次形劣化函数，该劣 化函数主要对弹性应变能进行劣化，以此可获得材料的总势能，其中总势能由弹性应变能和断裂能组成，依据能量最小变分原理，可确定相场控制方程。

具体地，所述材料劣化函数包括：

其中，d表示相场变量，函数

参数A的值大于等于1，B 的值大于0。

上述材料劣化函数的参数A和B同时控制着材料劣化函数的演化，劣化函数还具有g(1)＝0、g(0)＝1以及g′(1)＝0的特性。

在一个实施例中，所述参数输入值包括上一时步网格材料点的应变ε(t_n,x)和相场变 量d(t_n,x)，以及施加的宏观应变增量ΔE(t_n+1)和时间增量Δt_n+1；

需要初始化的数值计算参数包括：宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以及初 始相场变量d^0,0(t_n+1,x)，并计算初始相场历史变量

在一个实施例中，所述宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以及初始相场变量 d^0,0(t_n+1,x)计算方式分别为：

E(t_n+1)＝E(t_n)+ΔE(t_n+1)，

ε⁰(t_n+1,x)＝ε(t_n,x)+ΔE(t_n+1)，

d^0,0(t_n+1,x)＝d(t_n,x)，

式中，t_n表示上一时间步，t_n+1＝t_n+Δt_n+1表示本时间步，Δt_n+1为时间增量，x表示位置，E(t_n)表示上一时步的宏观应变，E(t_n+1)表示本时步的宏观应变，ΔE(t_n+1)表示本 时步内的宏观应变变化量，ε(t_n,x)表示材料点x在t_n时刻的应变，d(t_n,x)表示材料点x 在t_n时刻的相场变量值；

所述初始相场历史变量包括：

式中，ψ为弹性自由能，

及

分别为各向同性弹性自由能的正、负分量，k是 一微小正值，λ、μ为拉梅常数，εⁱ、n_i分别为应变张量ε的主应变和主应变方向，符 号

N表示维度，符号

表示卷积。

在一个实施例中，所述相场变量的迭代求解公式为：

式中，φ为相场变量d的函数，上标0代表初始FFT应力平衡迭代循环步，上标i 代表第i步斯特芬森加速固定点迭代；

网格材料点的应力计算公式为：

式中，I为单位二阶张量。

本实施例基于所述斯特芬森加速固定点迭代方法，首先在FFT框架下，求解相场变量d^0,i+1(t_n+1,x)，当计算满足迭代收敛条件后，结束迭代并更新初始相场变量d⁰(t_n+1,x)及计算材料点的初始应力σ⁰(t_n+1,x)。

进一步地，根据上述应力计算公式获得的应力场，进行应力平衡条件收敛性判断，如果满足收敛条件，进入下一步宏观应变增量的施加；否则，在FFT计算框架中，采用 交错求解方法，先根据上一迭代步的应变场，计算Lippmann-Schwinger周期性方程，更 新材料点应变场；再固定应变场，计算材料点相场历史变量

迭代计算本时步 材料点的相场变量d^j+1,i+1(t_n+1,x)及应力，直至满足迭代应力平衡收敛条件，结束迭代并 更新本时步的材料点的相场变量d^j+1(t_n+1,x)及应力σ^j+1(t_n+1,x)。

在一个示例中，相场变量的迭代求解公式也可以为：

式中，上标“j+1”代表第j+1步FFT应力平衡迭代，上标“i+1”代表第i+1步斯 特芬森加速固定点迭代。

相应地迭代步微观应力计算公式为：

式中各变量定义与前文所述相同。

在一个实施例中，所述宏观应力计算公式为：

式中，Ω表示计算空间的体积，t_n+1表示本时间步。

在一个实施例中，以定量计算某多孔材料断裂为例，对上述基于FFT法确定非均质岩石裂纹扩展路径的计算方法进行说明：

单胞孔隙率约为5％，空间离散为127×127×127分别率，弹性模量E_s的值为100GPa， 泊松比ν的值为0.15。

在步骤1中，不同于传统的相场所选择的二次形劣化函数，本实施例所选择的双参数材料劣化函数形式为：

式中，d表示相场变量，函数

参数A的值为2，B的值 为5，参数A取值对材料劣化函数的影响及参数B取值对材料劣化函数的影响分别如图 2、图3所示；

基于传统的相场模型，本实施例的相场控制方程形式为：

式中，g_c是材料临界表面能释放率，其值为2.0×10^-3kN/mm，l是裂纹开展控制参数，其值为0.015，η是材料粘性参数，其取值对一个包含中心球状空隙的单胞宏观轴向 应力-应变关系的影响如图4所示，

是相场历史变量。

在步骤2中，输入值(参数输入值)包括上一时步的微观应变ε(t_n,x)和相场变量 d(t_n,x)，以及宏观应变变化量ΔE(t_n+1)，其值为

时间增量Δt_n+1，其值 为25s；需要初始化的参数包括本时间步的宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以 及初始相场变量d^0,0(t_n+1,x)，并计算初始相场历史变量

参数计算形式分别为：

E(t_n+1)＝E(t_n)+ΔE(t_n+1)

ε⁰(t_n+1,x)＝ε(t_n,x)+ΔE(t_n+1)

d^0,0(t_n+1,x)＝d(t_n,x)

式中，t_n表示上一时间步，t_n+1＝t_n+Δt_n+1表示本时间步，Δt_n+1为时间增量，x表示位置，E(t_n)表示上一时步的宏观应变，E(t_n+1)表示本时步的宏观应变，ΔE(t_n+1)表示本 时步内的宏观应变变化量，ε(t_n,x)表示材料点x在t_n时刻的应变，d(t_n,x)表示材料点x 在t_n时刻的相场变量值；

初始相场历史变量的计算方式为：

其中，

式中，ψ为弹性自由能，

及

分别为各向同性弹性自由能的正、负分量，k是 一微小正值，λ、μ为拉梅常数，εⁱ、n_i分别为应变张量ε的主应变和主应变方向，符 号

N表示维度。

在步骤3中，基于所述交错求解方法，首先在FFT框架下，运用牛顿迭代法求解相场变量d^0,i+1(t_n+1,x)，当计算满足迭代收敛条件后，结束迭代并更新初始相场变量 d⁰(t_n+1,x)及计算初始微观应力σ⁰(t_n+1,x)；

相场变量的迭代求解公式为：

式中，φ为相场变量d的函数，变量的第一个上标“0”表示初始外循环迭代步，第二个上标“i”表示第i步内循环迭代。

初始微观应力计算公式为：

式中，I为单位二阶张量。

在步骤4中，基于所述交错求解方法，在FFT框架下，首先求解当前步的相场历史变量

再运用牛顿迭代法求解相场变量d^j+1,i+1(t_n+1,x)，当计算满足迭代收敛条件后，结束迭代并更新初始相场变量d^j+1(t_n+1,x)及计算初始微观应力σ^j+1(t_n+1,x)；

相场变量的迭代求解公式为：

式中，φ为相场变量d的函数，变量的第一个上标“j+1”表示第“j+1”步外循 环迭代，第二个上标“i”表示第i步内循环迭代。

初始微观应力计算公式为：

在步骤5中，按照

计算本时步的宏观应力，式中Ω表示计算空间的体积。

当η取值为10⁵Pa·s时，在单轴拉伸试验中，单胞的裂纹演化规律如图5所示。

本实施例所提供的基于FFT法确定非均质岩石裂纹扩展路径的计算方法与原有的基于有限元框架的相场方法相比，基于FFT均匀化方法的相场断裂计算方法，克服了有 限元框架的不足，在保证计算精度的同时，提高了数值计算效率，降低了计算成本；本 发明方法从微观角度出发，通过均匀化理论建立起微观尺度的裂纹开展和宏观尺度的材 料破坏之间的联系，能够准确地反映材料断裂的机理，为预测非均质材料断裂提供了可 靠依据。其易于程序化，并行化，便于通过编程语言加以实现，通过计算机进行断裂预 测，大大减少了人工计算量。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中 的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

需要说明的是，本申请实施例所涉及的术语“第一\第二\第三”仅仅是区别类似的对象，不代表针对对象的特定排序，可以理解地，“第一\第二\第三”在允许的情况下可 以互换特定的顺序或先后次序。应该理解“第一\第二\第三”区分的对象在适当情况下 可以互换，以使这里描述的本申请的实施例能够以除了在这里图示或描述的那些以外的 顺序实施。

本申请实施例的术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或模块的过程、方法、装置、产品或设备没有限定于已 列出的步骤或模块，而是可选地还包括没有列出的步骤或模块，或可选地还包括对于这 些过程、方法、产品或设备固有的其它步骤或模块。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不 能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的 保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (8)

1.一种基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，包括如下步骤：

S20，在非均质岩石材料的特征单元体计算网格上施加宏观应变增量，根据参数输入值，初始化数值计算参数，根据初始应变场计算相场初始相场历史变量；

S30，将初始相场历史变量代入相场控制方程，将相场控制方程进行快速傅里叶变换，推导出傅里叶空间中固定点相场迭代格式；基于斯特芬森加速固定点迭代方法，迭代计算相场变量，直至满足迭代收敛条件，结束迭代并进行傅里叶逆变换，计算真实空间中初始相场变量；将初始相场变量代入材料本构关系获得计算网格材料点的应力和应变初始值，根据应力和应变初始值确定应力场；

S40，根据应力场进行应力平衡条件收敛性判断，如果满足收敛条件，进入步骤S50中的宏观应变增量的施加；否则，在FFT计算框架中，采用交错求解方法，先根据应变初始值，计算Lippmann-Schwinger周期性方程，更新网格材料点应变场；再固定应变场，计算网格材料点相场历史变量，依据步骤S30迭代计算本时步网格材料点的相场变量及应力，直至满足迭代应力平衡收敛条件，结束迭代并更新本时步的网格材料点的相场变量及应力；

S50，根据各网格材料点的应力，确定所有网格材料点的体积平均，计算出特征单元体的宏观应力值，完成本时步的计算，返回执行步骤S20，并以本时步的宏观应力值作为步骤S20的参数输入值，直至加载到指定宏观应变值，以获得整个加载过程的相场d的演化云图，确定裂纹的扩展路径。

2.根据权利要求1所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，在步骤S20之前，还包括：

S10，基于相场模型的相场变量，根据材料类别，确定材料劣化函数及相场方法计算控制方程。

3.根据权利要求1所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，基于相场模型的相场变量，根据材料类别，确定材料劣化函数及相场方法计算控制方程包括：

基于相场模型的相场变量，根据材料类别，确定材料劣化函数，采用材料劣化函数对弹性应变能进行劣化，获得材料的总势能，使总势能包括弹性应变能和断裂能，以依据能量最小变分原理确定相场方法计算控制方程。

4.根据权利要求3所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，所述材料劣化函数包括：

其中，d表示相场变量，函数

参数A的值大于等于1，B的值大于0。

5.根据权利要求1所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，所述参数输入值包括上一时步网格材料点的应变ε(t_n,x)和相场变量d(t_n,x)，以及施加的宏观应变增量ΔE(t_n+1)和时间增量Δt_n+1；

需要初始化的数值计算参数包括：宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以及初始相场变量d^0,0(t_n+1,x)，并计算初始相场历史变量

6.根据权利要求1所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，所述宏观应变E(t_n+1)、初始微观应变ε⁰(t_n+1,x)以及初始相场变量d^0,0(t_n+1,x)计算方式分别为：

E(t_n+1)＝E(t_n)+ΔE(t_n+1)，

ε⁰(t_n+1,x)＝ε(t_n,x)+ΔE(t_n+1)，

d^0,0(t_n+1,x)＝d(t_n,x)，

式中，t_n表示上一时间步，t_n+1＝t_n+Δt_n+1表示本时间步，Δt_n+1为时间增量，x表示位置，E(t_n)表示上一时步的宏观应变，E(t_n+1)表示本时步的宏观应变，ΔE(t_n+1)表示本时步内的宏观应变变化量，ε(t_n,x)表示材料点x在t_n时刻的应变，d(t_n,x)表示材料点x在t_n时刻的相场变量值；

所述初始相场历史变量包括：

式中，ψ为弹性自由能，

及

分别为各向同性弹性自由能的正、负分量，k是一微小正值，λ、μ为拉梅常数，εⁱ、n_i分别为应变张量ε的主应变和主应变方向，符号

N表示维度。

7.根据权利要求1所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，所述相场变量的迭代求解公式为：

式中，φ为相场变量d的函数，上标0代表初始FFT应力平衡迭代循环步，上标i代表第i步斯特芬森加速固定点迭代；

网格材料点的应力计算公式为：

式中，I为单位二阶张量。

8.根据权利要求1所述的基于FFT法确定非均质岩石裂纹扩展路径的计算方法，其特征在于，所述宏观应力计算公式为：

式中，Ω表示计算空间的体积，t_n+1表示本时间步。

Images