一种临氢厚壁圆柱壳弹性应力应变的预测方法

技术领域

本发明涉及临氢承载结构的弹性响应预测领域，具体是基于氢脆的HELP理论预测厚壁圆柱壳的弹性应力应变。

背景技术

圆柱壳是氢系统中储氢和输氢过程中的关键承压元件，确定厚壁圆柱壳在氢环境下的弹性响应是设计储氢罐、加氢反应器等临氢设备时必须考虑的问题。由于氢损伤本构关系的复杂性，现有的预测材料在氢环境下的力学响应的方法对研究人员的素质要求很高，同时需要大量编程工作，难以在工程上推广。也没有出现专门用于预测临氢圆柱壳应力应变的方法。因此提出一种简单的预测临氢圆柱壳应力应变的方法具有工程意义。本发明将氢脆的氢促进塑性局部化理论(HELP)应用于厚壁圆柱壳，提出了一种用大气环境圆柱弹性应力应变预测临氢圆柱壳的弹性应力应变的方法。

发明内容

本发明的目的在于针对现有技术的不足，提供一种临氢厚壁圆柱壳弹性应力应变预测方法，该方法基于经典的大气环境中的圆柱壳弹性应力应变来预测临氢圆柱壳的弹性应力应变。

本发明的目的是通过以下技术方案来实现的：一种临氢厚壁圆柱壳弹性应力应变的预测方法，该方法包括以下步骤：

步骤1：按照以下公式计算厚壁圆柱壳的周向应力σ_θ及径向应力σ_r：

${\sigma }_{\theta }=\frac{a^{2}p}{b^2-a^2}(1+\frac{b^2}{r^2})$

${\sigma }_r=\frac{a^{2}p}{b^2-a^2}(1-\frac{b^2}{r^2})$

其中a为圆柱壳内半径，b为圆柱壳外半径，p为内压载荷，r为需要求取应力应变的位置。

步骤2：求解以下方程组，计算氢浓度c：

${\sigma }_{kk}=\frac{2(1+\nu )a^{2}p}{b^2-a^2}-\frac{E}{3}\frac{V_H}{V_M}(c-c_0)$

$c=(\frac{{\mathrm{\beta K}}_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}+\frac{{\mathrm{\alpha N}}_T}{N_L}\frac{K_T{\theta }_L}{1-{\theta }_L+K_T{\theta }_L})$

$K_L=\exp \left(\frac{{\sigma }_{kk}{V}_H}{3RT}\right)$

$K_T=\exp \left(\frac{W_B}{RT}\right)$

${\theta }_L=\frac{K_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}$

${\theta }_L^0=c_0/\beta$

其中，E为材料的弹性模量，ν为材料的泊松比，V_H表示氢在母材中的偏摩尔体积，V_M为母材的摩尔体积，α为每单位晶格的氢陷阱数，β每单位晶格的晶格间隙结合位点数，N_L＝N_A/V_M为每单位体积的金属原子数量，N_A为阿伏伽德罗常数，N_T表示每单位体积的陷阱个数。R为理想气体常数，T为开氏温标，c₀为无载荷时圆柱壳内的氢浓度，W_B为材料的氢陷阱结合能。

步骤3：计算轴向应力σ_z

${\sigma }_z=\nu ({\sigma }_{\theta }+{\sigma }_r)-\frac{E}{3}\frac{V_H}{V_M}(c-c_0)$

步骤4：计算周向应变ε_θ及径向应变ε_r

${ϵ}_{\theta }=\frac{1}{E}{\sigma }_{\theta }-\frac{\nu }{E}({\sigma }_r+{\sigma }_z)+\frac{1}{3}\frac{V_H}{V_M}(c-c_0)$

${ϵ}_r=\frac{1}{E}{\sigma }_r-\frac{\nu }{E}({\sigma }_{\theta }+{\sigma }_z)+\frac{1}{3}\frac{V_H}{V_M}(c-c_0).$

进一步地，所述步骤2采用牛顿法求解氢浓度c，具体包括以下子步骤：

步骤201：设置氢浓度初值c＝c₀；

步骤202：根据当前氢浓度c，按照下式计算各个参数：

${\sigma }_{kk}=\frac{2(1+\nu )a^{2}p}{b^2-a^2}-\frac{E}{3}\frac{V_H}{V_M}(c-c_0)$

$K_L=\exp \left(\frac{{\sigma }_{kk}{V}_H}{3RT}\right)$

$K_T=\exp \left(\frac{W_B}{RT}\right)$

${\theta }_L=\frac{K_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}$

${\theta }_L^0=c_0/\beta$

步骤203：计算以下函数值g：

$g=c-(\frac{{\mathrm{\beta K}}_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}+\frac{{\mathrm{\alpha N}}_T}{N_L}\frac{K_T{\theta }_L}{1-{\theta }_L+K_T{\theta }_L})$

步骤204：判断g值大小，若|g|≥ε_err，则执行步骤205至步骤206，否则结束计算，得到氢浓度c，其中ε_err为收敛容差，可取ε_err＝10^-6；

步骤205：按照下式计算氢浓度c₁

$c_1=c-\frac{g}{g^{\prime }}$

其中

$g^{\prime }=1+[\beta +\alpha \frac{N_T}{N_L}\frac{K_T}{{(1-{\theta }_L+K_T{\theta }_L)}^2}]\frac{{\theta }_L^0(1-{\theta }_L^0)}{{(1-{\theta }_L^0+K_L{\theta }_L^0)}^2}\frac{K_L{V}_H}{3RT}\frac{E}{3}\frac{V_H}{V_M}$

步骤206：令c＝c₁，返回步骤202。

本发明具有以下优点：采用求解非线性方程组的方法，根据易得的圆柱壳在大气环境中的弹性应力应变直接预测氢环境下圆柱壳的弹性应力应变，不需要编写复杂的氢损伤材料本构模型的有限元程序，应用门槛较低。

附图说明

图1为本发明的实施对象简图；

图2为本发明实例计算得到的应力在不同半径上的分布；

图3为本发明实例计算得到的应变在不同半径上的分布。

具体实施方式

以下以图1和表1所示的实例为实施对象，对本发明作进一步说明。

图1所示的实例是一个两端约束的厚壁圆柱壳，其内半径和外半径分别为a和b，承受着恒定内压载荷p，在无应力状态下浓度为c₀的氢在圆柱壳内均匀分布。本发明可以预测处于氢环境中的圆柱壳在内压p作用下的应力和应变。

表1实例用到的材料参数和几何参数

本发明方法的实现过程如下：

步骤1：选取一个半径r(a≤r≤b)，例如取r＝a，将表1中的内半径a，外半径b以及内压载荷p带入以下公式计算厚壁圆柱壳的周向应力σ_θ及径向应力σ_r：

${\sigma }_{\theta }=\frac{a^{2}p}{b^2-a^2}(1+\frac{b^2}{r^2})$

${\sigma }_r=\frac{a^{2}p}{b^2-a^2}(1-\frac{b^2}{r^2})$

步骤2：求解以下方程组，计算氢浓度c：

${\sigma }_{kk}=\frac{2(1+\nu )a^{2}p}{b^2-a^2}-\frac{E}{3}\frac{V_H}{V_M}(c-c_0)$

$c=(\frac{{\mathrm{\beta K}}_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}+\frac{{\mathrm{\alpha N}}_T}{N_L}\frac{K_T{\theta }_L}{1-{\theta }_L+K_T{\theta }_L})$

$K_L=\exp \left(\frac{{\sigma }_{kk}{V}_H}{3RT}\right)$

$K_T=\exp \left(\frac{W_B}{RT}\right)$

${\theta }_L=\frac{K_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}$

${\theta }_L^0=c_0/\beta$

其中各参数均列于表1中，步骤2具体包括以下子步骤：

步骤201：设置氢浓度初值c＝c₀；

步骤202：根据当前氢浓度c带入下式计算各个参数：

${\sigma }_{kk}=\frac{2(1+\nu )a^{2}p}{b^2-a^2}-\frac{E}{3}\frac{V_H}{V_M}(c-c_0)$

$K_L=\exp \left(\frac{{\sigma }_{kk}{V}_H}{3RT}\right)$

$K_T=\exp \left(\frac{W_B}{RT}\right)$

${\theta }_L=\frac{K_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}$

${\theta }_L^0=c_0/\beta$

步骤203：将步骤202得到的参数带入以下函数值，计算函数值g：

$g=c-(\frac{{\mathrm{\beta K}}_L{\theta }_L^0}{1-{\theta }_L^0+K_L{\theta }_L^0}+\frac{{\mathrm{\alpha N}}_T}{N_L}\frac{K_T{\theta }_L}{1-{\theta }_L+K_T{\theta }_L})$

步骤204：判断g值大小，若|g|≥ε_err，则执行步骤205至步骤206，否则结束计算，得到氢浓度c，其中收敛容差取ε_err＝10^-6；

步骤205：按照下式计算氢浓度c₁

$c_1=c-\frac{g}{g^{\prime }}$

其中

$g^{\prime }=1+[\beta +\alpha \frac{N_T}{N_L}\frac{K_T}{{(1-{\theta }_L+K_T{\theta }_L)}^2}]\frac{{\theta }_L^0(1-{\theta }_L^0)}{{(1-{\theta }_L^0+K_L{\theta }_L^0)}^2}\frac{K_L{V}_H}{3RT}\frac{E}{3}\frac{V_H}{V_M}$

所需各参数均由步骤202计算得到；

步骤206：令c＝c₁，返回步骤202。

步骤3：将步骤2得到的氢浓度c带入下式计算轴向应力σ_z

${\sigma }_z=\nu ({\sigma }_{\theta }+{\sigma }_r)-\frac{E}{3}\frac{V_H}{V_M}(c-c_0)$

步骤4：将步骤1和步骤3得到的周向应力σ_θ、径向应力σ_r和轴向应力σ_z带入下式计算周向应变ε_θ及径向应变ε_r

${ϵ}_{\theta }=\frac{1}{E}{\sigma }_{\theta }-\frac{\nu }{E}({\sigma }_r+{\sigma }_z)+\frac{1}{3}\frac{V_H}{V_M}(c-c_0)$

${ϵ}_r=\frac{1}{E}{\sigma }_r-\frac{\nu }{E}({\sigma }_{\theta }+{\sigma }_z)+\frac{1}{3}\frac{V_H}{V_M}(c-c_0).$

采用不同的r值重复步骤1至步骤4可以计算得到整个截面上的应力和应变分布，结果如图2、图3所示，另外可以看出氢浓度在圆柱壳内的分布是均匀的。