一种有效计算固体电解质断裂韧性的并行多尺度方法

1.本发明涉及一种多尺度计算方法，特别是计算考虑高温时力-电化学耦合场作用下原子到连续介质的模拟过程，计算材料的断裂韧性。


背景技术：

2.现有的原子到连续介质的计算方法可以计算纳米级别的均匀加载的断裂问题[1]以及热的扩散问题[2]，关于毫米级别的模型的断裂问题以及复杂边界条件的断裂问题都还没有合适的计算模型。
[0003]
[1]zimmerman,jonathan a.,et al."calculation of stress in atomistic simulation."modellingandsimulation in materials science andengineering 12.4(2004):s319.
[0004]
[2]wagner,gregory j.,et al."an atomistic-to-continuum coupling method for heat transfer in solids."computer methods inappliedmechanics and engineering 197.41-42(2008):3351-3365.


技术实现要素：

[0005]
本发明的目的是为了解决目前毫米级别的模型的断裂问题以及复杂边界条件的断裂问题还没有合适的计算模型，提供一种原子到连续介质的多尺度计算方法计算高温下受力-电化学耦合场作用的电解质的断裂韧性。
[0006]
为实现上述目的，本发明采取的技术方案如下：
[0007]
一种有效计算固体电解质断裂韧性的并行多尺度计算方法，所述方法为：
[0008]
(1)使用abaqus计算毫米级别的模型；
[0009]
(2)使用等效应力强度因子将毫米级别的模型转化为纳米级别的模型；
[0010]
(3)多尺度计算纳米级别模型的核心区域，边界条件由abaqus计算结果提供。步骤(3)得到使用多尺度计算得到的核心区域的应力应变图，根据这个会判断加载到多少时达到临界状态，再将这时的加载带入公式得到断裂韧性。
[0011]
进一步地，步骤一中，如图1所示模型为平面应变模型，采用平面八节点四边形单元，单轴拉伸，使用abaqus有限元软件计算得到。
[0012]
进一步地，所述步骤(2)具体为：
[0013]
等效应力强度因子计算公式为
[0014][0015]
其中，ki为等效应力强度因子，f为与加载和几何形状有关的参数，s为单轴拉伸加载的载荷大小，a是裂纹长度的一半；
[0016]
对于上述模型，f的计算公式为：
[0017][0018]
α＝a/b
ꢀꢀ
(3)
[0019]
α为a与b的比值，无实际意义，b为裂纹方向模型长度的一半；
[0020]
在加载时，通过保持等效应力强度因子不变，将微米模型转化为微观纳米模型。
[0021]
进一步地，所述步骤(3)具体为：
[0022]
(1)通过局部化函数来建立原子量和连续介质量之间的关系如下，
[0023]
ρ
*
(x,t)＝∑m
α
δ[x-x
α
(t)]
ꢀꢀ
(4)
[0024]
其中，ρ
*
(x,t)为原子量与连续介质量之间的关系，m
α
为原子α的质量，δ为dirac算子，x为连续介质节点，x
α
(t)为α原子在t时刻的坐标；
[0025]
(2)通过加权余量法的形式，建立起原子质量和连续介质节点量的关系，从而用原子质量表达节点量为，
[0026][0027]
其中，ρi是节点i处的连续介质量，ω是计算区域，||.....||2表示平方范数，ρ
*
是对应节点处的原子量，i是连续介质网格节点i，ni是节点i处的有限元形函数，v是积分区域，j分别代表连续介质网格节点j，nj是节点j处的有限元形函数，ρj是节点j处的连续介质量，m
ij
是中括号对应处公式，即x
α
为原子α的坐标，n
iα
为原子α坐标处的形函数值，是分别对j和α求和，为m
ij
的逆矩阵，
[0028]
设一质点系在理想的(线性或非线性)约束以及主动力的作用下从某一时刻出发，则对于符合约束的可能加速度，建立约束函数，
[0029][0030]
其中，mi是质点质量，ri是质点位移，是位移的两阶导数，fi是原子受力；
[0031]
利用高斯最小约束原理和拉格朗日乘子法，有
[0032][0033]
其中，是原子真实受力，是由势函数计算得到的原子受力，λi是拉格朗日乘子，是约束方程，约束方程的广义形式为
[0034][0035]
其中，为原子α的原子量，和分别为原子α的坐标和速度，ai是节点i处的连续介质节点量；
[0036]
采用高斯最小约束法，对方程(7)取极值，
[0037][0038][0039]
其中，φ是势函数；
[0040]
整理方程(10)，原子力表示成由标准分子动力学势函数和连续介质对代表性体积内原子影响的两部分：
[0041][0042]
力耦合约束表示代表体积内的总原子动量与对应的节点动量一致：
[0043][0044]
其中，是节点j处动量的导数；
[0045]
建立了原子量和连续介质量之间的双向信息传递过程，考虑整个系统的控制方程，把式(5)的系统动量分解成连续介质区域和原子区域两部分，再同时方程两边对时间求导：
[0046][0047]
ω
fe
是有限元区域、n
jα
是α原子处的形函数值，ρ是连续介质量，是连续介质节点速度；
[0048]
以上建立了多尺度模型计算的算法。
[0049]
本发明相对于现有技术的有益效果为：本发明能研究具有中心裂纹的gdc(gadolinium doped ceria，gdc，如gdc前面有系数则表示掺杂的钆离子占阳离子的比例，例如10gdc表示阳离子中钆离子数量占10％)电解质在不同温度下的单轴拉伸断裂行为。基于断裂力学理论，通过等效应力强度因子的方法将gdc电解质的宏观结构转化为微观中间过渡模型，并通过多尺度方法给出了分析gdc宏观结构断裂韧性的详细计算过程。此外，该多尺度方法研究了力-电化学耦合场下gdc宏观结构的断裂韧性。
附图说明
[0050]
图1为10gdc和20gdc非化学计量氧空位计算模型图；
[0051]
图2为本发明的方法流程图；
[0052]
图3为10gdc 800℃时裂尖附近非化学计量氧空位浓度分布图；
[0053]
图4为10gdc 900℃时裂尖附近非化学计量氧空位浓度分布图；
[0054]
图5为20gdc 800℃时裂尖附近非化学计量氧空位浓度分布图；
[0055]
图6为20gdc 900℃时裂尖附近非化学计量氧空位浓度分布图；
[0056]
图7为使用多尺度方法计算加载为0mpa温度为800℃时10gdc的应力应变图；
[0057]
图8为使用多尺度方法计算加载为0.4mpa温度为800℃时10gdc的应力应变图；
[0058]
图9为使用多尺度方法计算加载为0mpa温度为900℃时10gdc的应力应变图；
[0059]
图10为使用多尺度方法计算加载为0.4mpa温度为900℃时10gdc的应力应变图；
[0060]
图11为使用多尺度方法计算加载为0mpa温度为800℃时20gdc的应力应变图；
[0061]
图12为使用多尺度方法计算加载为0.4mpa温度为800℃时20gdc的应力应变图；
[0062]
图13为使用多尺度方法计算加载为0mpa温度为900℃时20gdc的应力应变图；
[0063]
图14为使用多尺度方法计算加载为0.4mpa温度为900℃时20gdc的应力应变图。
具体实施方式
[0064]
为使本发明的目的、技术方案及优点更加清楚明白，以下结合具体实施方式，对本发明进行进一步的详细说明。应当理解的是，此处所描述的具体实施方式仅用以解释本发明，并不限定本发明的保护范围。
[0065]
原子到连续介质的多尺度计算方法是一种将原子级别的模拟与连续介质级别的模拟相结合的方法，用于研究物质的性质和行为。在这种方法中，使用不同的计算技术和模型来描述原子级别和连续介质级别的行为，以便在不同尺度上获得准确和全面的信息。原子级别的模拟通常使用分子动力学(molecular dynamics，md)或量子力学方法，可以模拟原子之间的相互作用、原子的运动和化学反应等微观过程。这种模拟能够提供原子级别的细节信息，但由于计算资源和时间限制，其适用范围通常受到限制。连续介质级别的模拟使用连续介质力学(continuum mechanics)或有限元方法，将物质视为连续的，基于宏观物理规律建立数学模型。这种模拟可以描述大尺度的物质行为，如流体流动、固体变形等，但无法提供原子级别的细节。多尺度计算方法将这两种级别的模拟相结合，通过在原子级别和连续介质级别之间建立耦合关系，实现信息的交流和传递。这可以通过嵌入原子级别模拟结果到连续介质模型中，或者通过使用粗粒化模型将连续介质级别的信息传递到原子级别模拟中来实现。
[0066]
本发明的原理是：通过推导电解质中非化学计量氧空位浓度的分布，得到裂尖处非化学计量氧空位浓度。通过等效应力强度因子的方法，将毫米级别的模型在保持应力强度因子不变的前提下转化为纳米级别的模型，对纳米级别模型的核心区域采用多尺度计算方法，其他区域采用有限元软件abaqus进行计算，给多尺度计算区域的边界节点每个节点施加位移，位移大小由有限元软件abaqus计算得到。
[0067]
实施例1：
[0068]
计算高温下10gdc电解质在一侧低氧分压时的断裂韧性。
[0069]
如图1所示，这是计算10gdc模型中非化学计量氧空位浓度分布的模型图。平面应变模型长宽均为1毫米，中心裂纹长度为0.04毫米，模型左侧的低氧分压浓度为logp(o2)＝18。如图2所示，此图是使用多尺度方法计算断裂韧性的流程图。其中图2(a)是毫米级模型，图2(a)中模型与图1中模型相同。图2(b)是纳米级模型，通过保持应力强度因子不变的方法将图2(a)转化为图2(b)，图2(c)是裂纹附近网格划分，多尺度的边界条件由abaqus计算图
(b)中模型得到。图2(d)是多尺度计算示意图。使用图1中模型计算得到裂尖处的非化学计量氧空位浓度分布后，再采用图2中多尺度模型计算裂纹附近区域，最后计算可得到断裂韧性。
[0070]
(1)使用abaqus模型中各计算节点的位移。
[0071]
(2)通过节点位移计算节点的非化学计量氧空位浓度。
[0072]
(3)采用多尺度模型计算裂纹附近的区域，得到模型的断裂韧性。
[0073]
对图1模型进行加载，载荷范围是0-0.4mpa，载荷增幅是0.05mpa。电化学边界条件是氧分压为logp(o2)＝18，温度条件为为800℃和900℃时。计算结果如图3和4所示，图3是800℃时10gdc裂纹尖端的非化学计量氧空位浓度分布的计算结果，图4是900℃时10gdc裂纹尖端的非化学计量氧空位浓度分布的计算结果。选取加载为0mpa和0.4mpa，温度为800℃和900℃，使用多尺度模型计算此时裂尖核心区域的应力应变曲线，计算流程图如图2所示。计算结果如图7～10所示，图7～10分别表示加载为0mpa温度为800℃，加载为0.4mpa温度为800℃,加载为0mpa温度为900℃,加载为0.4mpa温度为900℃时多尺度计算得到的应力应变图，各图还给出了模拟初始阶段，应力最大值时，模拟结束时原子的构型图。读取应力最大值时，单轴拉伸加载的大小，将加载大小带入公式计算得到断裂韧性。10gdc在800℃和900℃的断裂韧性(k
ic
表示断裂韧性)如下表所示：
[0074][0075]
实施例2：
[0076]
计算高温下20gdc电解质在一侧低氧分压时的断裂韧性。
[0077]
如图1所示，这是计算20gdc模型中非化学计量氧空位浓度分布的模型图。平面应变模型长宽均为1毫米，中心裂纹长度为0.04毫米，模型左侧的低氧分压浓度为logp(o2)＝18。如图2所示，此图是使用多尺度方法计算断裂韧性的流程图。其中图2(a)是毫米级模型，图2(a)中模型与图1中模型相同。图2(b)是纳米级模型，通过保持应力强度因子不变的方法将图2(a)转化为图2(b)，图2(c)是裂纹附近网格划分，多尺度的边界条件由abaqus计算图(b)中模型得到。图2(d)是多尺度计算示意图。使用图1中模型计算得到裂尖处的非化学计量氧空位浓度分布后，再采用图2中多尺度模型计算裂纹附近区域，最后计算可得到断裂韧性。
[0078]
(1)使用abaqus模型中各计算节点的位移。
[0079]
(2)通过节点位移计算节点的非化学计量氧空位浓度。
[0080]
(3)采用多尺度模型计算裂纹附近的区域，得到模型的断裂韧性。
[0081]
对图1模型进行记载，载荷范围是0-0.4mpa，载荷增幅是0.05mpa。电化学边界条件是氧分压为logp(o2)＝18，温度条件为为800℃和900℃时。计算结果如图5和6所示，图5是800℃时20gdc裂纹尖端的非化学计量氧空位浓度分布的计算结果，图6是900℃时20gdc裂纹尖端的非化学计量氧空位浓度分布的计算结果。选取加载为0mpa和0.4mpa，温度为800℃
和900℃，使用多尺度模型计算此时裂尖核心区域的应力应变曲线，计算流程图如图2所示。计算结果如图11～14所示，图11～14分别表示加载为0mpa温度为800℃，加载为0.4mpa温度为800℃,加载为0mpa温度为900℃,加载为0.4mpa温度为900℃时多尺度计算得到的应力应变图，各图还给出了模拟初始阶段，应力最大值时，模拟结束时原子的构型图。读取应力最大值时，单轴拉伸加载的大小，将加载大小带入公式计算得到断裂韧性。20gdc在800℃和900℃的断裂韧性(k
ic
表示断裂韧性)如下表所示：
[0082]

技术特征：
1.一种有效计算固体电解质断裂韧性的并行多尺度计算方法，其特征在于：所述方法为：(1)使用abaqus计算毫米级别的模型；(2)使用等效应力强度因子将毫米级别的模型转化为纳米级别的模型；(3)多尺度计算纳米级别模型的核心区域，边界条件由abaqus计算结果提供。2.根据权利要求1所述的一种有效计算固体电解质断裂韧性的并行多尺度计算方法，其特征在于：步骤一中，采用平面八节点四边形单元，单轴拉伸，使用abaqus有限元软件计算得到。3.根据权利要求1所述的一种有效计算固体电解质断裂韧性的并行多尺度计算方法，其特征在于：所述步骤(2)具体为：等效应力强度因子计算公式为其中，k
i
为等效应力强度因子，f为与加载和几何形状有关的参数，s为单轴拉伸加载的载荷大小，a是裂纹长度的一半；对于上述模型，f的计算公式为：α＝a/b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)α为a与b的比值，无实际意义，b为裂纹方向模型长度的一半；在加载时，通过保持等效应力强度因子不变，将微米模型转化为微观纳米模型。4.根据权利要求1所述的一种有效计算固体电解质断裂韧性的并行多尺度计算方法，其特征在于：所述步骤(3)具体为：(1)通过局部化函数来建立原子量和连续介质量之间的关系如下，ρ
*
(x,t)＝∑m
α
δ[x-x
α
(t)]
ꢀꢀꢀꢀꢀꢀ
(4)其中，ρ
*
(x,t)为原子量与连续介质量之间的关系，m
α
为原子α的质量，δ为dirac算子，x为连续介质节点，x
α
(t)为α原子在t时刻的坐标；(2)通过加权余量法的形式，建立起原子质量和连续介质节点量的关系，从而用原子质量表达节点量为，其中，ρ
i
是节点i处的连续介质量，ω是计算区域，||.....||2表示平方范数，ρ
*
是对应节点处的原子量，i是连续介质网格节点i，n
i
是节点i处的有限元形函数，v是积分区域，j分别代表连续介质网格节点j，n
j
是节点j处的有限元形函数，ρ
j
是节点j处的连续介质量，m
ij
是中括号对应处公式，即x
α
为原子α的坐标，n
iα
为原子α坐标处的形函数值，是分别对j和α求和，为m
ij
的逆矩阵，设一质点系在理想的(线性或非线性)约束以及主动力的作用下从某一时刻出发，则对于符合约束的可能加速度，建立约束函数，其中，m
i
是质点质量，r
i
是质点位移，是位移的两阶导数，f
i
是原子受力；利用高斯最小约束原理和拉格朗日乘子法，有其中，是原子真实受力，是由势函数计算得到的原子受力，λ
i
是拉格朗日乘子，是约束方程，约束方程的广义形式为其中，为原子α的原子量，和分别为原子α的坐标和速度，a
i
是节点i处的连续介质节点量；采用高斯最小约束法，对方程(7)取极值，采用高斯最小约束法，对方程(7)取极值，其中，φ是势函数；整理方程(10)，原子力表示成由标准分子动力学势函数和连续介质对代表性体积内原子影响的两部分：力耦合约束表示代表体积内的总原子动量与对应的节点动量一致：其中，是节点j处动量的导数；建立了原子量和连续介质量之间的双向信息传递过程，考虑整个系统的控制方程，把式(5)的系统动量分解成连续介质区域和原子区域两部分，再同时方程两边对时间求导：ω
fe
是有限元区域、n
jα
是α原子处的形函数值，ρ是连续介质量，是连续介质节点速度；
以上建立了多尺度模型计算的算法。

技术总结
一种有效计算固体电解质断裂韧性的并行多尺度计算方法，所述方法为：使用ABAQUS计算毫米级别的模型；使用等效应力强度因子将毫米级别的模型转化为纳米级别的模型；多尺度计算纳米级别模型的核心区域，边界条件由ABAQUS计算结果提供。本发明能研究具有中心裂纹的GDC电解质在不同温度下的单轴拉伸断裂行为。基于断裂力学理论，通过等效应力强度因子的方法将GDC电解质的宏观结构转化为微观中间过渡模型，并通过多尺度方法给出了分析GDC宏观结构断裂韧性的详细计算过程。此外，该多尺度方法研究了力-电化学耦合场下GDC宏观结构的断裂韧性。韧性。韧性。


技术研发人员：杨志强 黄润泽
受保护的技术使用者：哈尔滨工业大学
技术研发日：2023.05.31
技术公布日：2023/9/14