(19)国家知识产权局 (12)发明 专利 (10)授权公告 号 (45)授权公告日 (21)申请 号 202111546154.7 (22)申请日 2021.12.16 (65)同一申请的已公布的文献号 申请公布号 CN 114330053 A (43)申请公布日 2022.04.12 (73)专利权人 北京应用物理与计算数 学研究所 地址 100083 北京市海淀区花园路6号 专利权人 北京理工大 学 (72)发明人 刘瑜　黄厚兵　宋海峰　盛杰　 武帅　杨超　王越超　韩国民　 孙博　施小明　 (74)专利代理 机构 北京荟英捷创知识产权代理 事务所(普通 合伙) 11726 专利代理师 张阳(51)Int.Cl. G06F 30/23(2020.01) G06F 17/13(2006.01) C25F 3/00(2006.01) G06F 119/14(2020.01) (56)对比文件 CN 111899819 A,2020.1 1.06 CN 113555072 A,2021.10.26 CN 113704997 A,2021.1 1.26 CN 112487568 A,2021.0 3.12 审查员 李孜孜 (54)发明名称 模拟金属表面氢化点腐蚀形貌演化的方法 及装置 (57)摘要 本发明公开了一种模拟金属表面氢化点腐 蚀形貌演化的方法及装置， 采用各向异性的弹性 常数建立了耦合各向异性弹性能的氢化点腐蚀 相场模型， 并利用插值函数的方法构造等效弹性 常数来表征体系的非均匀性， 适用于各向异性的 金属材料， 可模拟金属材料因弹性各向异性生长 导致的形貌差异， 更好地展示了氢化点腐蚀的生 长形貌特征； 在点腐蚀相场模型中， 采用有限元 方法求解非周期边界条件的力学平衡方程， 可以 避免傅里叶变换求解周期性弹性场引起非周期 性物理参数周期化的问题。 权利要求书4页 说明书10页 附图2页 CN 114330053 B 2022.12.02 CN 114330053 B 1.一种模拟金属表面氢化 点腐蚀形貌演化的方法， 其特 征在于， 包括以下步骤： 获取初始 物理参数； 所述初始 物理参数包括界面能、 动力学系数； 确定对应时间步下的相场和浓度场的分布参 量； 根据所述分布参 量计算所述对应时间步下的化学能驱动力； 基于有限元方法求解所述对应时间步下的力学平衡方程， 得到弹性应变和应力分量； 所述力学平衡方程包含预设弹性常数和预设边界条件， 所述预设弹性常数包括各向异性 弹 性常数， 所述预设边界条件 包括非周期边界条件； 根据所述弹性应 变和应力分量计算所述对应时间步下的弹性应 变能驱动力； 根据所述弹性应变能驱动力与所述化学能驱动力， 计算得到所述对应时间步下的总驱 动力； 将所述对应时间步下的总驱动力 输入相场演化和浓度扩散方程中进行演化求解， 以及 增加相应的时间总步长； 若判断所述时间总步长满足输出 条件， 则输出 所述对应时间步下的演化结果； 若判断所述 时间总步长不满足所述输出条件， 则将所述对应时间步下的演化结果作为 新的相场和浓度场的初始分布参数， 以进行新 一轮的演化计算， 直至满足所述输出 条件； 所述基于有限元方法求解所述对应时间步下的力学平衡方程， 得到弹性应变和应力分 量， 包括： 初始化全局矩阵， 所述全局矩阵包括刚度矩阵、 力矩阵； 获取金属基体及氢化物的弹性常数， 以及根据插值函数计算得到等效弹性常数； 所述 弹性常数为各向同性弹性常数或各向异性弹性常数； 将所述弹性常数、 所述 等效弹性常数转换成全局刚度矩阵； 应用初始设置的本征应 变或者应力矩阵， 应用全局载荷向量， 应用预设边界条件； 求解位移场分量， 以及根据所述 位移场分量计算得到弹性应 变和应力分量； 所述相场演化和浓度扩散方程， 包括： 其中， F是系统的总自由能， L是相界面迁移率的动力学系数， M是氢浓度的扩散迁移率， D是扩散系数， ΔGch和ΔGel分别是化学能驱动力和弹性应 变能驱动力； 其中， 总自由能F的计算公式如下： F＝Fbulk+Fint+Fel＝ ∫[fbulk(c, η )+fint+fel]dV 其中， Fbulk是化学自 由能， Fint是界面梯度能， Fel是弹性应变能， fbulk(c, η )是体自由能密 度， fint是与弥散界面有关的梯度能密度， fel是弹性应 变能密度， η是相场变量， c是 氢浓度；权　利　要　求　书 1/4 页 2 CN 114330053 B 2在KKS相场模型假设下， fbulk(c, η )写为 fbulk(c, η )＝h( η )fβ(cβ)+[1‑h( η )]fα(cα)+wg( η ) c＝h( η )cβ+[1‑h( η )]cα 其中， cα和cβ分别表示金属基体和氢化物中的氢的摩尔分数； h( η )是单调变化的插值函 数， h( η )＝3 η2‑2 η3； fα(cα)和fβ(cβ)分别是金属基体和氢化物的自由能密度， w是双势阱g( η ) ＝ η2(1‑η )2的高度； 其中， κη是与相场变量 η有关的梯度能系数； 其中， 和 分别是应力分量和弹性应 变。 2.根据权利要求1所述的方法， 其特征在于， 所述预设边界条件还包括周期边界条件； 所述非周期边界条件 包括固定、 载荷、 自由边界条件。 3.根据权利要求1所述的方法， 其特 征在于， 所述 等效弹性常数的计算公式如下： 其中， Cijkl( η )是等效弹性常数， 其为相场变量η的函数， 和 分别是金属基体 和氢化物的弹性常数， h( η )是 单调变化的插值 函数， ΔCijkl是 和 之差。 4.根据权利要求1 ‑3任一项所述的方法， 其特征在于， 基于微弹性理论计算弹性应变能 在金属表面氢化 点腐蚀形貌演化中的物理效应， 如下： 其中， 和 分别是弹性应变和应力分量， εij(r)是位移场变量u的分量ui和uj表示 的总应变， 与晶格本征应变 相关， δij是克罗内克 符号， h( η )是单调变化的插值函权　利　要　求　书 2/4 页 3 CN 114330053 B 3