CN107025339B - 一种位错对铁电材料畴结构影响机理的分析方法和系统 - Google Patents

Abstract

本发明实施例公开了一种位错对铁电材料畴结构影响机理的分析方法和系统，该方法包括：步骤S1，建立包含压电效应的第一压电效应计算模型；步骤S2，基于位错在铁电材料中引起的应力应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型；步骤S3，基于包含位错的第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；步骤S4，基于包含位错的第二压电效应计算模型和挠曲电效应计算模型，生成极化矢量图和/或极化云图。本发明实施例与传统的仅考虑压电效应的位错相场模型相比，考虑了挠曲电效应，获得了与实验结果更吻合的模拟结果。

Description

一种位错对铁电材料畴结构影响机理的分析方法和系统

技术领域

本发明属于铁电材料模拟技术领域，特别涉及一种位错对铁电材料畴结构影响机理的分析方法和系统。

背景技术

传统半导体存储器主要有两大体系：易失性存储器和非易失性铁电存储器。

非易失性铁电存储器(FRAM)是由美国Ramtron公司生产的非易失性铁电介质读写存储器。其核心是铁电晶体材料，铁电晶体材料使得铁电存储产品同时拥有随机存储器(RAM)和非易失性存储器的特性。非易失性铁电存储器因其众多优良的特性有着广泛的应用潜力，被誉为最有潜力的存储器。然而，非易失性铁电存储器存在一系列亟待解决的失效问题(例如，铁电材料的退化问题)，这些失效问题都与铁电材料的极化(即畴结构)有关。铁电材料中存在的缺陷(例如，位错)通常会是极化翻转过程中新畴结构的潜在成核点，因而深入研究缺陷对极化翻转的影响，对提高铁电性能和防止铁电存储器的退化来说是迫切需要的。

位错是最常见的一种缺陷，它会在铁电材料制备过程中不可避免地产生。由于实验无法直观测量到位错对极化(即畴结构)的具体影响作用，因而常通过理论模型来进行内在机理的研究。

为了研究位错对铁电材料极化(即畴结构)的具体影响作用，Alpay等人通过基于Landau-Devonshire形式的热力学模型研究位错处的极化梯度对铁电材料铁电性能的影响，发现铁电材料中的位错会形成影响电畴翻转的死层，该死层也会减弱铁电材料其他铁电性能；Zheng等人通过相场理论研究铁薄膜/基底界面处位错对铁电材料电学性能的影响，发现界面位错会在高温下极大地减小铁电材料的平均极化，从而数倍地减小有效的居里温度，并对局部介电性能产生非常大的影响。

虽然以上研究中所采用模型均考虑了极化与位错应变场的耦合，即压电效应，它们在一定程度能揭示位错与铁电材料极化的相互作用。但是由于位错附近存在巨大的应变梯度，这个巨大的应变梯度会与极化通过挠曲电耦合效应发生耦合，最终产生挠曲电极化，从而显著改变附近极化场的分布。而现有的位错理论模型只考虑了压电效应，均未考虑挠曲电效应，挠曲电效应的忽略会导致位错在铁电材料极化(即畴结构)演化及失效中所起作用被严重低估。

因此，有必要提出一种同时考虑压电效应和挠曲电效应的理论模型，以观察铁电材料内不同角度的刃型位错对铁电畴结构的影响，使得模拟结果更准确。

发明内容

本发明的目的是提供一种位错对铁电材料畴结构影响机理的分析方法和系统，本发明针对现有位错相场模型存在的问题，提供一种简单、准确的模拟方法，通过分析不同柏氏矢量方向的位错对铁电材料畴结构变化的影响，揭示了现有技术中较难直接测量到的位错影响机理。

为实现上述目的，本发明的一个方面提供了一种位错对铁电材料畴结构影响机理的分析方法，所述方法包括：

步骤S1，建立包含压电效应的第一压电效应计算模型；

步骤S2，基于位错在铁电材料中引起的应力应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型；

步骤S3，基于包含位错的第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；

步骤S4，基于包含位错的第二压电效应计算模型和所述挠曲电效应计算模型，生成极化矢量图和/或极化云图。

根据本发明的另一方面，提供了一种位错对铁电材料畴结构影响机理的分析系统，包括：第一压电模型建立单元、第二压电模型建立单元、挠曲电模型建立单元和极化图生成单元；

第一压电模型建立单元，用于建立包含压电效应的第一压电效应计算模型；

第二压电模型建立单元，用于基于位错在铁电材料中引起的应力场、应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型；

挠曲电模型建立单元，用于基于包含位错的第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；

极化图生成单元，用于基于包含位错的第二压电效应计算模型和所述挠曲电效应计算模型，生成极化矢量图和/或极化云图。

如上所述，本发明的模拟方法简单、高效，通过模拟系统分析了不同柏氏矢量方向的位错对铁电材料畴结构变化的影响，揭示了现有技术中较难直接测量到的位错影响机理；与传统的仅考虑压电效应的位错相场模型相比，该方法加入了挠曲电效应的作用，获得了与实验结果更吻合的模拟结果。

附图说明

图1是本发明实施例一提供的位错对铁电材料畴结构影响机理的分析方法的流程示意图；

图2是根据本发明实施例二提供的建立包含铁电材料的相场几何模型的流程示意图；

图3是根据本发明实施例三提供的步骤S1的流程示意图；

图4是本发明实施例二提供的在COMSOL Multiphysics软件中建立的相场几何模型的结构示意图；

图5是本发明实施例二提供的在COMSOL Multiphysics软件中模拟的铁电材料在初始时刻的极化示意图；

图6是本发明实施例四提供的步骤S2的流程示意图；

图7是本发明实施例提供的包含

方向刃型位错的相场几何模型的示意图；

图8是本发明实施例提供的实验观察到的

方向刃型位错附近的极化示意图；

图9是本发明实施例提供的考虑压电效应后模拟出的位错附近的极化示意图；

图10是本发明实施例五提供的步骤S3的流程示意图；

图11是本发明实施例五提供的考虑挠曲电效应后模拟出的位错附近极化示意图；

图12是与图11对应的位错附近的极化分量P₂的云图；

图13是本发明实施例提供的包含[100]方向刃型位错的相场几何模型的示意图；

图14是本发明实施例提供的考虑挠曲电效应后模拟出的[100]方向位错附近极化示意图；

图15是本发明实施例提供的与图15对应的[100]方向位错附近极化分量P₂的云图；

图16是本发明实施例提供的

方向刃型位错的相场几何模型的示意图；

图17是本发明实施例提供的考虑挠曲电效应后模拟出的

方向刃型位错附近极化示意图；

图18是与图18对应的

方向刃型位错附近极化分量P₂的云图；

图19是本发明实施例提供的a畴中含有位错的相场几何模型的示意图；

图20是本发明实施例提供的实验观测到的在a畴中位错的存在导致a畴中产生了一个c畴的示意图；

图21是本发明实施例模拟出的极化示意图；

图22是本发明实施例七提供的位错对铁电材料畴结构影响机理的分析系统的结构示意图；

图23是本发明实施例八提供的相场几何模型建立单元5的结构示意图；

图24是本发明实施例九提供的第一压电模型建立单元1的结构示意图；

图25是本发明实施例十提供的第二压电模型建立单元2的结构示意图；

图26是本发明实施例十一提供的挠曲电模型建立单元3的结构示意图；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

需要说明的是：

铁电畴(简称电畴)是指铁电体中偶极子有序排列、自发极化方向一致的区域。根据薄膜中电畴的数量，可分为单畴和多畴两种。按照极化方向和铁电材料平面位向可将铁电畴分为a畴和c畴两大类。极化轴垂直于铁电材料平面的畴称c畴，极化轴平行于铁电材料的畴称a畴。

在材料科学中，位错指晶体材料的一种内部微观缺陷，即原子的局部不规则排列(晶体学缺陷)。从几何角度看，位错属于一种线缺陷，可视为晶体中已滑移部分与未滑移部分的分界线，其存在对材料的物理性能，尤其是力学性能，具有极大的影响。理想位错主要有两种形式:刃型位错(edgedislocations)和螺旋位错(screw dislocations)，混合位错(mixed dislocations)兼有刃型位错和螺旋位错两者的特征。

柏氏矢量是描述位错性质的一个重要物理量。柏氏矢量的确定方法：

1.人为假定位错线方向；一般是从纸背向纸面或由上向下为位错线正向；

2.用右手螺旋法则来确定柏格斯回路的旋转方向；使位错线的正向与右螺旋的正向一致；

3.将含有位错的实际晶体和理想的完整晶体相比较；在实际晶体中作柏氏回路，在完整晶体中按相同的路线和步法作回路，实际晶体中的回路必然不重合，路线终点指向起点的矢量，即“柏氏矢量”。极化(polarization)，指具有铁电性的铁电材料电偶极子的正负电荷中心不重合而产生的电偶极矩。

居里温度(T_C)，是指铁电体从铁电相转变成顺电相的相变温度。低于居里点温度时该物质成为铁电体。当温度高于居里点温度时，该物质成为顺电体。

压电效应：某些电介质在沿一定方向上受到外力的作用而变形时，其内部会产生极化现象，同时在它的两个相对表面上出现正负相反的电荷。当外力去掉后，它又会恢复到不带电的状态，这种现象称为正压电效应。当作用力的方向改变时，电荷的极性也随之改变。相反，当在电介质的极化方向上施加电场，这些电介质也会发生变形，电场去掉后，电介质的变形随之消失，这种现象称为逆压电效应。通常将正压电效应和逆压电效应都简称为压电效应。

应力，是指物体由于外因(受力、湿度、温度场变化等)而变形时，在物体内各部分之间产生相互作用的内力，以抵抗这种外因的作用，并试图使物体从变形后的位置恢复到变形前的位置，这种物体内部的作用力为内力，也就是我们通常所说的应力。

影响畴结构变化的主要因素：

(1)静电能：对于一个绝缘的晶体，它的极化引起面电荷，从而增加了静电能，如果晶体中出现若干畴就会使静电能降低。静电能的极小值是与晶体的形状有关，所以，晶体中的畴结构就依赖于晶体的形状。此外，晶体的电导率以及包围晶体的介质的电导率也影响着静电能，从而影响晶体内的电畴结构。

(2)应力场的能量：铁电体中的自发极化总是同自发形变相伴随，多畴的出现可以降低应力场的能量。如果晶体中有着不完整性，这种不完整性一方面可以增加晶体的电导率，从而降低极化场能量；另一方面，它们的存在使晶体内部区域处于受力状态，从而附加了应力场的能量，因此，畴的结构趋于使静电场能和应力场能皆尽可能减少的趋势。

再者，由于畴结构，畴同畴之间出现分界面(畴壁)。因此，附加了界面能，亦即畴壁能。还可概括以下几部分：

(3)偶极矩相互作用能：在畴壁上极化的方向是随处而异的，因而存在着偶极矩的相互作用能。

(4)各项异性能：在铁电体中，某些轴向的极化对应于最小能量，但畴的存在，往往使极化不能取得最有利的方向，因而增加了各项异性能，较之铁磁的情形，铁电体中各项异性能的单位值要大得多。

(5)弹性能：在铁电相时，晶体总是具有压电性的，在畴壁上，压电形变偏离于自发极化，因而出现弹性能。

实施例一

图1是本发明实施例一提供的位错对铁电材料畴结构影响机理的分析方法的流程示意图。

如图1所示，本发明实施例一提供的位错对铁电材料畴结构影响机理的分析方法，该方法包括：

步骤S1，建立包含压电效应的第一压电效应计算模型。本步骤中，首先，计算相场几何模型的第一自由能密度，根据第一自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程，进一步，将位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场的弱形式，设置铁电材料的材料参数，再将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何计算域的初始值和边界条件，然后对几何计算域划分网格，得到第一压电效应计算模型。

该步骤的具体执行过程可参见下图3所示。

步骤S2，基于位错在铁电材料中引起的应力场、应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型。

本步骤中，在第一压电效应计算模型的基础上，考虑位错在铁电材料中引起的应力场和应变场，建立包含位错的第二压电效应计算模型。

步骤S3，基于包含位错的第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型。

本步骤中，在包含位错的第二压电效应计算模型的基础上，考虑挠曲电效应能量密度方程，建立挠曲电效应计算模型。

步骤S4，基于包含位错的第二压电效应计算模型和挠曲电效应计算模型，生成极化矢量图和/或极化云图。

本步骤中，根据包含位错的第二压电效应计算模型，生成包含位错的极化矢量图；根据挠曲电效应计算模型生成包含挠曲电效应的极化矢量图；进一步，生成包含位错的极化云图。

本发明实施例通过建立包含压电效应的第一压电效应计算模型，根据位错在铁电材料中的引起的应力场、应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型；根据包含位错的第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；进一步，根据包含位错的第二压电效应计算模型和挠曲电效应计算模型，生成极化矢量图和/或极化云图，通过本发明实施例的技术方案能够直观测量到位错影响机理。与传统的仅考虑压电效应的位错模型相比，该方法加入了挠曲电效应的作用，获得了与实验结果更吻合的模拟结果。

实施例二

图2是根据本发明实施例二提供的建立包含铁电材料的相场几何模型的流程示意图。

如图2所示，在本发明实施例一的基础上，在步骤S1之前，还包括：建立包含铁电材料的相场几何模型，其包括：

步骤S01，获取位错在铁电材料中的第一坐标位置及柏氏矢量方向。

本步骤中，获取某个位错(铁电材料中可能存在多个位错)在铁电材料中的第一坐标位置及该位错的柏氏矢量方向。第一坐标位置为位错在铁电材料中的坐标位置，柏氏矢量为位错在铁电材料中的方向，柏氏矢量方向的确定，可以通过现有技术中确定柏氏矢量方向的方法，此处不再赘述。

在本发明实施例中，位错的第一坐标位置及柏氏矢量方向的获取方式包括但不限于：通过显微镜可以观察到位错在铁电材料中的第一坐标位置及位错的柏氏矢量方向，也可以通过其他方式获取位错在铁电材料中的第一坐标位置及位错的柏氏矢量方向。

步骤S02，建立坐标系。

本步骤中，建立坐标系，例如，可以以铁电材料下表面的中间位置为坐标原点建立坐标系，横轴为x轴，纵轴为y轴，并规定向右和向上为两轴(x轴和y轴)的正方向。

这里，步骤S01和步骤S02的执行顺序不唯一，可以先后执行也可以同时执行。

步骤S03，基于位错在铁电材料中的第一坐标位置，确定位错在坐标系中的第二坐标位置。

本步骤中，第二坐标位置为位错在坐标系中的坐标位置。假设，位错在铁电材料中的第一坐标位置为(x₁，x₂)，根据该第一坐标位置(x₁，x₂)，确定位错在该坐标系中的第二坐标位置。

步骤S04，基于位错在铁电材料中的柏氏矢量方向，确定位错在该坐标系中的几何方向。

本步骤中，根据位错在铁电材料中的柏氏矢量方向，确定位错在坐标系中的几何方向。

本发明实施例中，柏氏矢量为位错在铁电材料中的方向，柏氏矢量方向的确定，可以通过现有技术中确定柏氏矢量方向的方法，此处不再赘述。几何方向为位错在坐标系中与横轴(x轴)正方向的夹角。

这里，步骤S03和步骤S04的执行顺序不唯一，可以先后执行也可以同时执行。

步骤S05，基于坐标系、第二坐标位置和几何方向，建立包含铁电材料的相场几何模型。

本步骤中，在上述坐标系的基础上，根据位错的第二坐标位置和几何方向，创建包含铁电材料的相场几何模型，该相场几何模型用于模拟位错对铁电材料极化的影响机理。

如图4所示，图4是本发明实施例二提供的在COMSOL Multiphysics软件中建立的相场几何模型的结构示意图，为了方便模拟本发明的实验结果，本发明在COMSOLMultiphysics软件中建立相场几何模型(该相场几何模型是初始形态)，该相场几何模型的尺寸为80*30，该相场几何模型对应的实际尺寸为80nm*30nm，宽度W为80nm，高度H为30nm。如图5所示，图5是本发明实施例二提供的在COMSOL Multiphysics软件中模拟的铁电材料在初始时刻的极化示意图。

本发明实施例中，相场几何模型的形状不限于上述大小的长方形，可以是其他大小的长方形或正方形，只要是能够实现本发明对相场几何模型的要求即可。

本发明实施例通过获取位错在铁电材料中的第一坐标位置和柏氏矢量方向，确定位错在坐标系中的第二坐标位置和几何方向，在坐标系的基础上，基于位错的第二坐标位置和几何方向，建立包含铁电材料的相场几何模型。

实施例三

图3是根据本发明实施例三提供的步骤S1的流程示意图。

如图3所示，在本发明实施例一的基础上，在本发明实施例三提供的步骤S1进一步包括：

步骤S11，对铁电材料的体自由能密度、弹性应变能密度、极化与应变的耦合能量密度、梯度能密度和静电能密度求和，得到相场几何模型的第一自由能密度h。

具体地，根据下式计算第一自由能密度h：

h＝f_landau+f_strain+f_coup+f_grad+f_elec(1.1)

其中，f_landau表示体自由能密度，f_strain表示弹性应变能密度，f_coup表示极化与应变的耦合能量密度，f_grad表示梯度能密度，f_elec表示静电能密度。

已知，体自由能密度f_landau表示为：

其中，α₁表示介电刚度系数，介电刚度系数α₁可以由居里-外斯定律α₁＝(T-T_C)/2ε₀C₀得到，其中，T表示温度，T_C表示居里温度，ε₀表示真空介电常数，C₀表示居里外斯常数；α_ij，β_ijkl，γ_ijklmn表示唯象朗道-德文希尔系数，根据不同的材料取不同的值；P_i、P_j、P_k、P_l、P_m、P_n表示不同方向上的极化分量，例如，P_i表示极化在i方向上的分量，P_j、P_k、P_l、P_m、P_n依次类推。

已知，弹性应变能密度f_strain表示为：

其中，c_ijkl表示弹性系数，ε_ij、ε_kl表示不同方向的应变分量。

已知，极化与应耦变的合能量密度f_coup表示为：

其中，q_ijkl表示电致伸缩系数。

已知，梯度能密度f_grad表示为：

其中，g_ijkl表示梯度系数，P_i,j、P_k,l表示不同极化分量在不同方向上的极化梯度，x_j表示坐标。

已知，静电能密度f_elec表示为：

其中，ε表示介电常数，ε₀表示真空介电常数，E_i表示电场。

步骤S12，基于第一自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程。

具体的，根据式子(1.1)，分别对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程：

根据式子(1.1)，对位移场求偏导，得到位移场的本构方程为：

其中，σ_ij表示薄膜的应力。

根据式子(1.1)，对电场求偏导，得到电场的本构方程为：

其中，D_i表示电位移。

根据式子(1.1)，对极化场求偏导，得到极化场的本构方程为：

其中，η_i表示局部有效电场力，Λ_ij表示高阶局部电场力。

综上，通过式子(1.8)-(2.2)可知，位移场、电场、极化场完全耦合。

这里举例说明一下：

对于式子(2.1)，当i＝j＝1时，以此类推，可以求出相场几何模型下的σ₁₁、σ₁₂、σ₂₂、D₁、D₂、η₁、η₂、Λ₁₁、Λ₁₂、Λ₂₁、Λ₂₂的具体形式。根据剪应力互等定律，可知σ₁₂＝σ₂₁。

下面结合COMSOL软件说明本构方程写入COMSOL软件的具体过程：

将上面求解出的本构方程(1.8)-(2.2)写入COMSOL软件中：在ModelBulider→Component 1→Definitions，右键单击Definitions，选择Variables，在右侧输入各个变量。左边第一列Name，表示变量名称，如σ₁₁在这里写为sigma11，变量名只要不重复，不是x，y，z，t等软件内默认字符即可随意定义。第二列Expression即为变量对应的具体表达式，如sigma11的表达式在Expression下写为c11*e11+c12*e22-q11*p1^2-q12*p2^2，以此类推。第三列Unit表示变量的量纲，本模型采用的是无量纲，这一列可以不写。第四列为Description，即为对变量的描述，也可以不写，这一列的目的就是为了防止变量过多情况下混淆，sigma11的Description可以记为Stressxx。

此外还要写入其他变量，因为在这个模型中，因变量是下面提到的位移场中的位移u、v，电场中的电势V，极化场中的极化P₁、P₂。需要添加的变量有：应变

电场E_i＝-V_,i，应变梯度

需要注意的是，例如在写时，左边第一列Name可记为e11，右边Expression列对应的要写成ux，即表示u对坐标x求偏导；

可分别写为uxx，即表示u对x求二阶偏导数。还需要注意的是乘号不能少，比如写e12y，其Expression就是0.5*(uyy+vxy)，而不能写成0.5(uyy+vxy)。

步骤S13，将位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在计算域内积分，得到三个物理场的弱形式。

已知，位移场的力平衡方程为：

σ_ik,i+b_k＝0(i,j,k,l＝1,2)(2.3)

其中，σ_ik表示薄膜的应力，b_k表示铁电材料的体力(一个物体所受到的外力可分为两种：一种是分布在物体表面的作用力，称为面力，例如一个物体对另一个物体的压力，水压力等。另一种是分布在物体体积内的力，称为体力，如重力，磁场力，惯性力等)，此处不考虑铁电材料的体力，故体力b_k＝0。

由麦克斯韦方程组可知，电场的平衡方程为：

D_i,i-ξ＝0 (2.4)

在电场的平衡方程(2.5)中，不考虑铁电材料中的载流子，因此体电荷密度ξ＝0。故电场的平衡方程也可简写为：

D_i,i＝0 (2.5)

其中，D_i,i表示电位移的梯度(需要说明的是：电介质中，电位移是由外电场强度和电极化强度两部分叠加构成。电荷分布在物体内部时，单位体积内的电量称为体电荷密度。)

已知，对于极化场的方程：根据Landau-Devonshire热力学唯象理论，总自由能F可以展开为：

h＝f_landau+f_strain+f_coup+f_grad+f_elec (1.1)

如前所述，h表示第一自由能密度，f_landau表示体自由能密度，f_strain表示弹性应变能密度，f_coup表示极化与应变的耦合能量密度，f_grad表示梯度能密度，f_elec表示静电能密度，V表示体积。

对于铁电材料而言，极化翻转过程中，自发极化矢量P是随着时间而改变的，故极化场从非平衡态到平衡态的动态过程，可以用与时间有关的Ginzburg-Landau方程描述为：

r＝(x₁+x₂+x₃) (2.8)

其中，P_i表示极化，r为空间矢量，t表示时间，L表示动力学常数，F表示总自由能，δF/δP_i(r,t)表示热力学的驱动力，x₁、x₂为面内水平和薄膜厚度方向，x₃表示垂直于x₁、x₂的方向。P＝(P₁,P₂,P₃)为极化矢量，用三个方向的极化场来描述畴的结构。

下面计算三个物理场的弱形式：

其中，位移场的平衡方程为：σ_ik,i+b_k＝0，位移场的试函数为δu_x，则位移场的弱形式为：

电场的平衡方程为：D_i,i-ξ＝0，电场的试函数为

则电场的弱形式为

极化场的平衡方程为：

极化场的试函数为δP_i，则极化场的弱形式为：

步骤S14，设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何域的初始值和边界条件；

步骤S15，对几何计算域划分网格，得到第一压电效应计算模型。

下面结合COMSOL软件说明如何添加铁电材料的材料参数：

不同的铁电材料对应不同的材料参数，预先将每种铁电材料对应的材料参数进行关联。因为本模型中并没有在Material库中选择具体的铁电材料，而是通过自定义铁电材料PbTiO₃的某些材料参数来完成铁电材料的选择，因此需要在Global Definitions下的Parameters中输入材料参数。输入框分为四列，依次为Name、Expression、Value和Description。将(1h)-(1l)中提到的材料参数无量纲化处理后输入。

需要注意的是，第一列参数Name可以随意设定，只要不是前面提到的u、v、V、P₁、P₂及x、y、t等软件系统内默认具有特殊意义的符号均可。第二列Expression中输入参数的数值，第三列则自动显示数值，第三列不可编辑。第四列为对该参数的描述，仅仅为方便回顾，可以不用填写。

下面结合COMSOL软件说明如何添加位移场、电场和极化场的弱形式和因变量：

在Model下添加位移场、电场和极化场的弱形式，然后位移场定义两个因变量u和v，电场定义1个变量V，极化场定义两个变量P₁，P₂。

这里，为本模型中用到的无量纲化处理方法：采用室温下PbTiO₃的自发极化P_s＝0.757C/m²、梯度能系数g₁₁₀＝1.73×10^-10m⁴N/C²和朗道体自由能系数α₀＝-1.725×10⁸m²N/C²作为基本量进行无量纲化。

以下列出了本发明实施例所使用的所有物理量的量纲和无量纲方程式：

其中，上标星号*所表示的是对应的无量纲符号。例如PbTiO₃的朗道系数α₁₁＝-6.6489e7m⁶N/C⁴，则其对应的无量纲系数α₁*₁＝α₁₁*P₀ ²/|α₀|＝-0.2209。

有量纲和无量纲形式的各场变量间相互关系如下：

x^*＝x/(1×10^-9m)

P^*＝P/(0.757Cm^-2)

E^*＝E/(1.3058×10⁸C^-1N)

(▽ε)^*＝▽ε/(1×10⁹m^-1)

需要说明的是：模拟中的无量纲长度“1”对应真实空间中的1nm。无量纲极化值“1”对于真实空间中的0.757Cm²。

无量纲处理后的参数如下表1所示，注意表1中为唯象朗道-德文希尔系数α_ij的具体展开形式的无量纲值，

为唯象朗道-德文希尔系数β_ijkl的具体展开形式的无量纲值，

为唯象朗道-德文希尔系数γ_ijklmn的具体展开形式的无量纲值，

为弹性系数c_ijkl具体展开形式的Voigt简写记法的无量纲值，

为电致伸缩系数q_ijkl具体展开形式的Voigt简写记法的无量纲值，

为梯度系数g_ijkl具体展开形式的Voigt简写记法的无量纲值。

表1 Pb(Zr_0.1Ti_0,9)O₃的无量纲系数

下面结合COMSOL软件说明如何设置几何计算域的边界条件和初始值：

位移场设置为左右位移边界条件为周期性连续边界条件，上下表面为应力自由，即为自由边界条件。位移场初始值u₀＝v₀＝0。

电场设置左右电势边界为周期性连续边界条件，上下表面的电势设置为V＝0，即为短路边界条件。电势的初始值V₀＝0。

极化场设置左右极化边界为周期性连续边界条件，上下表面的极化边界条件考虑有较大的尺寸效应，即设置为-dP₂/dx₂都为0，极化场初始值为P₁＝0，P₂＝1。

本步骤中，基于位移场、电场和极化场的弱形式，对几何计算域划分网格，每个网格单元是1*1的正方形网格，相场几何模型进行80*30的网格划分，对应着真实空间的80nm*30nm，计算每个网格的电势，得到整个几何计算域的电势。

本发明实施例通过计算相场几何模型的第一自由能密度，基于第一自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程；进一步，将位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在计算域内积分，得到三个物理场的弱形式，设置材料参数后再将位移场、电场和极化场的本构方程分别代入至三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何计算域的初始值和边界条件，然后对几何计算域划分网格，得到第一压电效应计算模型。

实施例四

图6是本发明实施例四提供的步骤S2的流程示意图。

如图6所示，本发明在实施一的基础上，本发明实施例四提供的步骤S2进一步包括：

步骤S21，计算位错在铁电材料中引起的应力场。

根据下式计算位错引起的应力

其中，G表示剪切模量，b表示柏氏矢量的模，ν表示泊松比。

为了模拟不同柏氏矢量方向的位错，引进极坐标，将上述式子(4.4)改为：

最终可得位错在薄膜内产生的位移如下：

步骤S22，计算位错在铁电材料中引起的应变场。

已知，位错引起的应变

的本构方程为：

其中，E表示弹性模量，ν表示泊松比，已知，Θ＝σ_x+σ_y+σ_z＝σ_kk，当i＝j时，δ_ij取1，当i≠j时，δ_ij取0，δ_ij称为“kronecker delta”，具有对称性，δ_ij＝δ_ji。

根据下式计算位错引起的应变ε_ij ^d：

整理得到：

其中，θ表示柏氏矢量与x₁正方向的夹角，θ＝3π/2。

步骤S23，基于第一压电效应计算模型、位错的应力场和应变场，建立包含位错的第二压电效应计算模型。

本步骤中，基于第一压电效应计算模型、位错引起的应力场和位错引起的应变场，设置总时长和时间步长，以确保使第一压电效应计算模型达到稳定状态，选择求解器，计算包含位错的第二压电效应计算模型，图7是本发明实施例提供的包含

方向刃型位错的相场几何模型的示意图、图13是本发明实施例提供的包含[100]方向刃型位错的相场几何模型的示意图、图16是本发明实施例提供的

方向刃型位错的相场几何模型的示意图、图19是本发明实施例提供的a畴中含有位错的相场几何模型的示意图。图8是本发明实施例提供的实验观察到的

方向刃型位错附近的极化示意图、图20是本发明实施例提供的实验观测到的在a畴中位错的存在导致a畴中产生了一个c畴的示意图。如图9所示，图9是本发明实施例提供的考虑压电效应后模拟出的位错附近的极化示意图、图14是本发明实施例提供的考虑挠曲电效应后模拟出的[100]方向位错附近极化示意图、图17是本发明实施例提供的考虑挠曲电效应后模拟出的

方向刃型位错附近极化示意图。

本发明实施例通过将位错在铁电材料中引起的应力场、应变场与已经建立的第一压电效应计算模型结合起来，建立包含位错的第二压电效应计算模型。

实施例五

图10是本发明实施例五提供的步骤S3的流程示意图。

如图10所示，在本发明实施例一的基础上，本发明实施例五的步骤S3进一步包括：

步骤S31，对第一自由能密度和挠曲电能量密度求和，得到第二自由能密度f。

f＝f_landau+f_strain+f_coup+f_grad+f_elec+f_flexo (5.2)

其中，f₁₁、f₁₂、f₄₄分别为极化与纵向、横向、剪切应变梯度之间的耦合系数。

步骤S32，根据第二自由能密度对位移场、电场和极化场的变量求偏导，得到位移场、电场和极化场的本构方程；

在原本构方程(1.8)-(2.2)的基础上，考虑挠曲电效应导致的变化，考虑挠曲电效应之后，原位移场、电场和极化场的本构方程(1.8)-(2.2)变为：

其中，η′_i、D′_i分别表示广义电场和电位移，而t′_ijl、Λ′_ij分别与应变梯度

及极化梯度(

表示极化梯度)互为共轭量(即t′_ijl分别与应变梯度及和极化梯度互为共轭量，Λ′_ij分别与应变梯度及和极化梯度互为共轭量)。

S33，将含挠曲电效应的位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场含挠曲电效应的弱形式；

此处考虑挠曲电效应后，除位移场的弱形式发生变化外，电场和极化场的弱形式均与实施例三中保持一致。

考虑挠曲电效应后，位移场的平衡方程为：(σ_ik-t′_ijk,j)_,i+b_k＝0，位移场的试函数为δu_x，则位移场的弱形式为：

S34，设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式并添加至铁电材料的几何计算域。

S35，对几何计算域划分网格，得到第二压电效应计算模型。

新建一组参数，在原有参数的基础上加上挠曲电耦合系数，其中将f₁₁、f₁₂、f₄₄的无量纲值分别设置为-2，-20，-10。这一步的边界条件、初始值和实施例三相同。

如图11所示，图11是本发明实施例五提供的考虑挠曲电效应后模拟出的位错附近极化示意图。

本发明实施例通过计算包含挠曲电能量密度的第二自由能密度，根据第二自由能密度对位移场、电场和极化场的本构方程的变量求偏导，得到位移场、电场和极化场的本构方程，然后将含挠曲电效应的位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场含挠曲电效应的弱形式，设置材料参数后，再将位移场、电场和极化场的本构方程代入三个物理场的弱形式并添加至铁电材料的几何计算域，然后对几何计算域划分网格，得到第二压电效应计算模型。

实施例六

在本发明实施例一的基础上，本发明实施例六提供的步骤S4进一步包括：

基于包含位错的第二压电效应计算模型和挠曲电效应计算模型，生成包含位错的极化矢量图和/或挠曲电效应的极化矢量图和/或极化云图。

下面结合COMSOL软件说明包含位错的极化矢量图和挠曲电效应的极化矢量图的生成过程：

在桌面新建一个Excel文档，将上面导出的Polar_without_dislo.txt数据导入到Excel中去。将表头文件都删去，留下A-D四列，分别为x、y、VectorX(表示P₁大小)、VectorY(表示P₂大小)。在E列生成一列新数据，该数据数值是x+P_l，在F列生成y+P₂。然后，打开Origin作图软件，将Excel下的x、y、x+P_l、y+P₂都复制到Origin中，四列数据分别设定为X、Y、X、Y，然后点击Plot→specialized–vectorXYXY即可生成极化矢量图，最终导出图片，即图11所示。

同样，分别用相同方法处理另外两个极化矢量图数据，生成考虑单位错的极化矢量图和/或挠曲电的极化矢量图。

同样，用excel将P₂.txt中的表头删掉，将三列数据(x、y和Color)粘贴到Origin中去，三列分别命名为X、Y和Z，选中三列数据，然后点击Plot→Cntour→XYZ contour即可生成极化云图，最终导出极化云图，即12所示，图12是与图11对应的位错附近的极化分量P₂的云图。通过图11和图12可以看出，传统仅考虑压电效应模拟出的位错对极化的影响与实验不能够很好地吻合。相反，考虑了挠曲电效应之后，在位错附近形成了一个狭长的a畴，该a畴贯穿了薄膜的厚度方向，与实验观测吻合得很好。图15是本发明实施例提供的与图14对应的[100]方向位错附近极化分量P₂的云图，θ＝0。图18是与图17对应的

方向刃型位错附近极化分量P₂的云图，角度θ＝7π/4。图21是本发明实施例模拟出的极化示意图，通过图21可以看出，与图20实验观察到的极化情况很吻合。

实施例七

图22是本发明实施例七提供的位错对铁电材料畴结构影响机理的分析系统的结构示意图。

如图22所示，本发明实施例七提供的一种位错对铁电材料畴结构影响机理的分析系统，所述系统包括：

第一压电模型建立单元1，用于建立包含压电效应的第一压电效应计算模型。

第二压电模型建立单元2与第一压电模型建立单元1连接，用于基于位错在铁电材料中引起的应力场、应变场以及所述第一压电效应计算模型，建立包含位错的第二压电效应计算模型。

挠曲电模型建立单元3与第二压电模型建立单元2连接，用于基于包含位错的所述第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；

极化图生成单元4分别与第二压电模型建立单元2和挠曲电模型建立单元3连接，用于基于包含位错的所述第二压电效应计算模型和所述挠曲电效应计算模型，生成极化矢量图和/或极化云图。

上述各模块的执行过程可参见前述方法部分相应的步骤，此处不再赘述。

图23是本发明实施例八提供的相场几何模型建立单元5的结构示意图。

如图23所示，在一实施方式中，所述相场几何模型建立单元5包括：

位错信息获取单元101，用于获取位错在铁电材料中的第一坐标位置及柏氏矢量方向；

坐标建立单元102，用于建立坐标系；

位错位置确定单元103与位错信息获取单元101连接，用于基于位错在铁电材料中的第一坐标位置，确定位错在坐标系中的第二坐标位置；

位错方向确定单元104与位错信息获取单元101连接，用于基于位错在铁电材料中的柏氏矢量方向，确定位错在坐标系中的几何方向；

相场模型建立单元105分别与坐标建立单元102、位错位置确定单元103和位错方向确定单元104连接，基于坐标系、第二坐标位置和几何方向，建立包含铁电材料的相场几何模型。

上述各模块的执行过程可参见前述方法部分相应的步骤，此处不再赘述。

图24是本发明实施例九提供的第一压电模型建立单元1的结构示意图。

如图24所示，在一实施方式中，所述第一压电模型建立单元1进一步包括：

第一自由能密度计算单元11，用于对铁电材料的体自由能密度、弹性应变能密度、极化与应变的耦合能量密度、梯度能密度和静电能密度求和，得到所述相场几何模型的第一自由能密度；

第一本构方程确定单元12与第一自由能密度计算单元11连接，用于基于第一自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程；

第一弱形式计算单元13与第一本构方程确定单元12连接，用于将位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场的弱形式；

第一参数设置单元14分别与第一本构方程确定单元12和第一弱形式计算单元13连接，用于设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式，并添加至铁电材料的几何计算域，设置几何计算域的初始值和边界条件；

第一网格划分单元15与第一参数设置单元14连接，用于对几何计算域划分网格，得到第一压电效应计算模型。

上述各模块的执行过程可参见前述方法部分相应的步骤，此处不再赘述。

图25是本发明实施例十提供的第二压电模型建立单元2的结构示意图。

如图25所示，在一实施方式中，所述第二压电模型建立单元2进一步包括：

应力计算单元21，用于计算位错在铁电材料中引起的应力场；

应变计算单元22，用于计算位错在铁电材料中引起的应变场；

第二压电模型生成单元23分别与应力计算单元21和应变计算单元22连接，用于基于所述第一压电效应计算模型、位错的应力场和应变场，建立包含位错的第二压电效应计算模型。

上述各模块的执行过程可参见前述方法部分相应的步骤，此处不再赘述。

图26是本发明实施例十一提供的挠曲电模型建立单元3的结构示意图。

如图26所示，在一实施方式中，所述挠曲电模型建立单元3进一步包括：

第二自由能密度计算单元31，用于对第一自由能密度和挠曲电能量密度求和，得到第二自由能密度；

第二本构方程确定单元32与第二自由能密度计算单元31连接，用于基于第二自由能密度对位移场、电场和极化场的变量求偏导，得到位移场、电场和极化场的本构方程；

第二弱形式计算单元33，将含挠曲电效应的位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场含挠曲电效应的弱形式；

第二参数设置单元34分别与第二本构方程确定单元32和第二弱形式计算单元33连接，设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何域的初始值和边界条件；

第二网格划分单元35，用于对几何计算域划分网格，得到第二压电效应计算模型。

上述各模块的执行过程可参见前述方法部分相应的步骤，此处不再赘述。

在一实施方式中，所述极化图生成单元4进一步包括：

极化矢量图生成单元41，用于基于第一压电效应计算模型和挠曲电效应计算模型，生成包含位错的极化矢量图和/或挠曲电效应的极化矢量图和/或极化云图。

上述各模块的执行过程可参见前述方法部分相应的步骤，此处不再赘述。

如上所述，详细介绍了本发明的位错对铁电材料畴结构影响机理的分析方法和系统，与现有的技术相比，本发明除了采纳前人工作优点之外，还考虑了传统压电效应和挠曲电效应，本发明可观察铁电材料内不同角度的刃型位错对铁电畴结构的影响，模拟结果更准确。

本发明实施例提供的方法简单、高效，本发明实施例分析了不同柏氏矢量方向的位错对铁电材料畴结构变化的影响，揭示了现有技术中较难直接测量到的位错影响机理；与传统的仅考虑压电效应的位错模型相比，该方法加入了挠曲电效应的作用，获得了与实验结果更吻合的模拟结果。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (10)

1.一种位错对铁电材料畴结构影响机理的分析方法，所述方法包括：

步骤S1，建立包含压电效应的第一压电效应计算模型；

步骤S2，基于位错在铁电材料中引起的应力场、应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型；

步骤S3，基于包含位错的所述第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；

步骤S4，基于包含位错的所述第二压电效应计算模型和所述挠曲电效应计算模型，生成极化矢量图和/或极化云图。

2.根据权利要求1所述的方法，其中，在所述建立包含压电效应的第一压电效应计算模型步骤之前，还包括：建立包含铁电材料的相场几何模型，其包括：

步骤S01，获取位错在铁电材料中的第一坐标位置及柏氏矢量方向；

步骤S02，建立坐标系；

步骤S03，基于位错在铁电材料中的所述第一坐标位置，确定位错在所述坐标系中的第二坐标位置；

步骤S04，基于位错在铁电材料中的所述柏氏矢量方向，确定位错在所述坐标系中的几何方向；

步骤S05，基于所述坐标系、第二坐标位置和几何方向，建立包含铁电材料的相场几何模型。

3.根据权利要求2所述的方法，其中，所述建立包含压电效应的第一压电效应计算模型的步骤包括：

步骤S11，对铁电材料的体自由能密度、弹性应变能密度、极化与应变的耦合能量密度、梯度能密度和静电能密度求和，得到所述相场几何模型的第一自由能密度；

步骤S12，基于第一自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程；

步骤S13，将位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场的弱形式；

步骤S14，设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入至三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何计算域的初始值和边界条件；

步骤S15，对几何计算域划分网格，得到第一压电效应计算模型。

4.根据权利要求1所述的方法，其中，所述建立包含位错的第二压电效应计算模型的步骤包括：

步骤S21，计算位错在铁电材料中引起的应力场；

步骤S22，计算位错在铁电材料中引起的应变场；

步骤S23，基于所述第一压电效应计算模型、位错的应力场和应变场，建立包含位错的第二压电效应计算模型。

5.根据权利要求3所述的方法，其中，所述建立挠曲电效应计算模型的步骤包括：

步骤S31，对第一自由能密度和挠曲电能量密度求和，得到第二自由能密度；

步骤S32，基于第二自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程；

步骤S33，将含挠曲电效应的位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场含挠曲电效应的弱形式；

步骤S34，设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何计算域的初始值和边界条件；

步骤S35，对几何计算域划分网格，得到第二压电效应计算模型。

6.一种位错对铁电材料畴结构影响机理的分析系统，所述系统包括：

第一压电模型建立单元(1)，用于建立包含压电效应的第一压电效应计算模型；

第二压电模型建立单元(2)，用于基于位错在铁电材料中引起的应力场、应变场和第一压电效应计算模型，建立包含位错的第二压电效应计算模型；

挠曲电模型建立单元(3)，用于基于包含位错的第二压电效应计算模型和挠曲电效应能量密度方程，建立挠曲电效应计算模型；

极化图生成单元(4)，用于基于包含位错的所述第二压电效应计算模型和所述挠曲电效应计算模型，生成极化矢量图和/或极化云图。

7.根据权利要求6所述的系统，其中，所述系统还包括：相场几何模型建立单元(5)，用于建立包含铁电材料的相场几何模型，其包括：

位错信息获取单元(101)，用于获取位错在铁电材料中的第一坐标位置及柏氏矢量方向；

坐标建立单元(102)，用于建立坐标系；

位错位置确定单元(103)，用于基于位错在铁电材料中的第一坐标位置，确定位错在坐标系中的第二坐标位置；

位错方向确定单元(104)，用于基于位错在铁电材料中的柏氏矢量方向，确定位错在坐标系中的几何方向；

相场模型建立单元(105)，基于坐标系、第二坐标位置和几何方向，建立包含铁电材料的相场几何模型。

8.根据权利要求7所述的系统，其中，所述第一压电模型建立单元(1)包括：

第一自由能密度计算单元(11)，用于对铁电材料的体自由能密度、弹性应变能密度、极化与应变的耦合能量密度、梯度能密度和静电能密度求和，得到所述相场几何模型的第一自由能密度；

第一本构方程确定单元(12)，用于基于第一自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程；

第一弱形式计算单元(13)，将位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在计算域内积分，得到三个物理场的弱形式；

第一参数设置单元(14)，用于设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入至三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何计算域的初始值和边界条件；

第一网格划分单元(15)，用于对几何计算域划分网格，得到第一压电效应计算模型。

9.根据权利要求6-8中任一项所述的系统，其中，所述第二压电模型建立单元(2)包括：

应力计算单元(21)，用于计算位错在铁电材料中引起的应力场；

应变计算单元(22)，用于计算位错在铁电材料中引起的应变场；

第二压电模型生成单元(23)，用于基于所述第一压电效应计算模型位错的应力场和应变场，建立包含位错的第二压电效应计算模型。

10.根据权利要求8所述的系统，其中，所述挠曲电模型建立单元(3)包括：

第二自由能密度计算单元(31)，用于对第一自由能密度和挠曲电能量密度求和，得到第二自由能密度；

第二本构方程确定单元(32)，基于第二自由能密度对位移场、电场和极化场求偏导，得到位移场、电场和极化场的本构方程；

第二弱形式计算单元(33)，将含挠曲电效应的位移场、电场和极化场各自的平衡方程分别与三个物理场的试函数相乘并在几何计算域内积分，得到三个物理场含挠曲电效应的弱形式；

第二参数设置单元(34)，设置铁电材料的材料参数，将位移场、电场和极化场的本构方程分别代入三个物理场的弱形式并添加至铁电材料的几何计算域，设置几何域的初始值和边界条件；

第二网格划分单元(35)，用于对几何计算域划分网格，得到第二压电效应计算模型。

Images