含损伤的推进剂蠕变型本构模型的构建与有限元应用方法

摘要

本发明公开了一种含损伤的推进剂蠕变型本构模型的构建与有限元应用方法，其中含损伤的推进剂蠕变型本构模型的构建方法包括：获取推进剂蠕变参数；预设推进剂有效承载面积相对变化量；根据获取的蠕变参数构建推进剂蠕变型本构模型；根据有效承载面积相对变化量确定损伤变量；根据确定的损伤变量确定损伤发展方程；根据确定的损伤发展方程和构建的推进剂蠕变型本构模型，构建考虑损伤的推进剂蠕变本构模型；其中含损伤的推进剂蠕变型本构模型的限元应用方法是通过分解、离散上述构建方法构建的考虑损伤的推进剂蠕变本构模型，推导出有限元材料子程序二次开发需要的一致切线刚度阵。本发明能够计算推进剂药柱结构在贮存阶段蠕变损伤力学行为。

权利要求

1.一种考虑损伤的推进剂蠕变本构模型的构建方法，其特征是，包括以下步骤：获取推进剂蠕变参数；

预设推进剂有效承载面积相对变化量；

根据获取的蠕变参数构建推进剂蠕变型本构模型；

根据有效承载面积相对变化量确定损伤变量；

根据确定的损伤变量确定损伤发展方程；

根据确定的损伤发展方程和构建的推进剂蠕变型本构模型，构建考虑损伤的推进剂蠕变本构模型；

所述推进剂蠕变型本构模型包括(1)‑(3)式，具体如下：其中，εij为机械应变，Jijkl()为蠕变张量函数，ξ为时间t的折算时间，ξ`为第一被积函数τ的折算时间，σkl为推进剂损伤前的机械应力，τ为第一被积函数，t为加载时间，T()为温度，t`为第二被积函数，aT为温度移位因子，C1为推进剂WLF方程第一参数，C2为推进剂WLF方程第二参数，T为当前温度，Tr为参考温度；

所述构建考虑损伤的推进剂蠕变本构模型包括通过(10)式更新机械应变，具体如下：其中， 为有效应力张量函数。

2.根据权利要求1所述的考虑损伤的推进剂蠕变本构模型的构建方法，其特征是，所述损伤变量包括(4)式，具体如下：其中，D为损伤变量，A为推进剂材料初始承载面积，为推进剂材料发生损伤后的有效承载面积。

3.根据权利要求2所述的考虑损伤的推进剂蠕变本构模型的构建方法，其特征是，所述损伤变量满足(5)式，具体如下：其中，σ为推进剂损伤前的机械应力，为推进剂损伤后的机械应力。

4.根据权利要求3所述的考虑损伤的推进剂蠕变本构模型的构建方法，其特征是，所述损伤发展方程包括(6)‑(9)式，具体如下：其中，D0为推进剂初始损伤，DΔV为扩展形成的新损伤，σth为损伤应力阈值，γ推进剂损伤第一参数，K推进剂损伤第二参数，β为推进剂损伤第三参数，h()为应力状态函数，x为应力状态函数的输入参数。

5.一种考虑损伤的推进剂蠕变本构模型的有限元应用方法，其特征是，包括以下步骤：分解权利要求1‑4任一项构建的考虑损伤的推进剂蠕变本构模型，获得偏应变张量和球应变张量；

分别离散偏应变张量和球应变张量，并根据离散后的偏应变张量和球应变张量推导考虑损伤的推进剂蠕变本构模型的增量形式；

根据考虑损伤的推进剂蠕变本构模型的增量形式推导有限元材料子程序二次开发需要的一致切线刚度阵；

对增量形式和一致切线刚度阵进行编程，得到应用程序，并使用有限元软件调用所述应用程序。

6.根据权利要求5所述的考虑损伤的推进剂蠕变本构模型的有限元应用方法，其特征是，所述分解考虑损伤的推进剂蠕变本构模型，获得偏应变张量和球应变张量包括(11)式，具体如下：其中， 为有效应力偏张量函数， 为有效应力球张量函数，eij()为偏应变函数，为球应变张量函数，J()为推进剂蠕变柔量函数，ν是推进剂泊松比；

所述推进剂蠕变柔量函数包括(12)式，具体如下：

其中，J0为初始蠕变柔量，NJ为蠕变柔量级数的项数，t为加载时间，Jn为第n项蠕变柔量， 为第n项蠕变柔量系数；

所述有效应力偏张量函数包括(13)式，具体如下：

其中，Sij()为偏应力张量函数；

所述有效应力球张量函数包括(14)式，具体如下：

其中，σkk()为球应力张量函数。

7.根据权利要求6所述的考虑损伤的推进剂蠕变本构模型的有限元应用方法，其特征是，所述考虑损伤的推进剂蠕变本构模型的增量形式包括(15)式，具体如下：其中，Δεij(tm+1)为tm+1时刻应变张量的增量，Δeij(tm+1)为tm+1时刻应变偏量的增量，为tm+1时刻球应变张量的增量，δij为第四克罗内克符号。

8.根据权利要求7所述的考虑损伤的推进剂蠕变本构模型的有限元应用方法，其特征是，所述一致切线刚度阵包括(16)‑(17)式，具体如下：其中，Cijkl()为切线刚度张量函数，Δσij(tm+1)为tm+1时刻应力增量，Δεkl(tm+1)为tm+1时刻应变张量的增量，C2222(tm+1)为tm+1时刻切线刚度张量的第一分量，C2233(tm+1)为tm+1时刻切线刚度张量的第二分量，C2323(tm+1)为tm+1时刻切线刚度张量的第三分量，C3333(tm+1)为tm+1时刻切线刚度张量的第四分量，C1133(tm+1)为tm+1时刻切线刚度张量的第五分量，C1313(tm+1)为tm+1时刻切线刚度张量的第六分量，C1111(tm+1)为tm+1时刻切线刚度张量的第七分量，C1122(tm+1)为tm+1时刻切线刚度张量的第八分量，C1212(tm+1)为tm+1时刻切线刚度张量的第九分量；

所述切线刚度张量函数tm+1时刻切线刚度张量的所有分量中，除了第一分量至第九分量，其他分量均为零。

说明书

含损伤的推进剂蠕变型本构模型的构建与有限元应用方法

技术领域

[0001] 本发明涉及一种含损伤的推进剂蠕变型本构模型的构建与有限元应用方法，属于含能材料本构模型技术领域。

背景技术

[0002] 固体火箭发动机是远程火箭弹的动力装置，固体推进剂是固体发动机的重要组成部分。远程火箭弹在长时间贮存过程中，会使得药柱结构发生蠕变和损伤效应。对于药柱结构的蠕变损伤力学行为必须采用含损伤的推进剂蠕变型粘弹性本构模型，目前商业有限元软件只包含简单的线粘弹性本构模型，不能精确地描述推进剂蠕变损伤效应。

[0003] 试验结果表明，药柱结构在贮存期间会产生较大的蠕变变形，同时会产生一定的损伤效应。固体推进剂的损伤力学行为与当前的应力应变状态有着直接的联系。目前存在的推进剂损伤本构模型大都是在推进剂松弛型本构模型基础上建立的，难以直接反映推进剂在贮存阶段的蠕变力学行为。为了实现药柱结构在贮存阶段的蠕变损伤行为的精细化建模，亟需构建推进剂含损伤蠕变型本构模型。

[0004] 现有的商业有限元软件中，粘弹性本构模型都是不考虑材料的损伤力学行为的，难以对实际的情况进行准确地分析。考虑到现有商业有限元软件可以提供二次开发接口，可以利用该接口实现蠕变损伤本构模型的有限元应用。

发明内容

[0005] 本发明的目的在于克服现有技术中的不足，提供一种含损伤的推进剂蠕变型本构模型的构建与有限元应用方法，能够计算推进剂药柱结构在贮存阶段蠕变损伤力学行为。

[0006] 为达到上述目的，本发明是采用下述技术方案实现的：

[0007] 一方面，本发明提供一种考虑损伤的推进剂蠕变本构模型的构建方法，包括以下步骤：

[0008] 获取推进剂蠕变参数；

[0009] 预设推进剂有效承载面积相对变化量；

[0010] 根据获取的蠕变参数构建推进剂蠕变型本构模型；

[0011] 根据有效承载面积相对变化量确定损伤变量；

[0012] 根据确定的损伤变量确定损伤发展方程；

[0013] 根据确定的损伤发展方程和构建的推进剂蠕变型本构模型，构建考虑损伤的推进剂蠕变本构模型。

[0014] 进一步地，所述推进剂蠕变型本构模型包括(1)‑(3)式，具体如下：

[0015]

[0016]

[0017]

[0018] 其中，εij为机械应变，Jijkl()为蠕变张量函数，ξ为时间t的折算时间，ξ`为第一被积函数τ的折算时间，σkl为推进剂损伤前的机械应力，τ为第一被积函数，t为加载时间，T()为温度，t`为第二被积函数，aT为温度移位因子，C1为推进剂WLF方程第一参数，C2为推进剂WLF方程第二参数，T为当前温度，Tr为参考温度。

[0019] 进一步地，所述根据获取的蠕变参数确定损伤变量包括通过(4)式，具体如下：

[0020]

[0021] 其中，D为损伤变量，A为推进剂材料初始承载面积， 为推进剂材料发生损伤后的有效承载面积。

[0022] 进一步地，所述损伤变量满足(5)式，具体如下：

[0023]

[0024] 其中，σ为推进剂损伤前的机械应力， 为推进剂损伤后的机械应力。

[0025] 进一步地，所述损伤发展方程包括(6)‑(9)式，具体如下：

[0026]

[0027]

[0028]

[0029]

[0030] 其中，D0为推进剂初始损伤，D△V为扩展形成的新损伤，σth为损伤应力阈值，γ推进剂损伤第一参数，K推进剂损伤第二参数，β为推进剂损伤第三参数，h()为应力状态函数，x为应力状态函数的输入参数。

[0031] 进一步地，所述构建考虑损伤的推进剂蠕变本构模型包括通过(10)式更新机械应变，具体如下：

[0032]

[0033] 其中， 为有效应力张量函数。

[0034] 另一方面，本发明提供一种考虑损伤的推进剂蠕变本构模型的有限元应用方法，包括以下步骤：

[0035] 分解上述构建方法构建的考虑损伤的推进剂蠕变本构模型，获得偏应变张量和球应变张量；

[0036] 分别离散偏应变张量和球应变张量，并根据离散后的偏应变张量和球应变张量推导考虑损伤的推进剂蠕变本构模型的增量形式；

[0037] 根据考虑损伤的推进剂蠕变本构模型的增量形式推导有限元材料子程序二次开发需要的一致切线刚度阵；

[0038] 对增量形式和一致切线刚度阵进行编程，得到应用程序，并使用有限元软件调用所述应用程序。

[0039] 进一步地，所述分解考虑损伤的推进剂蠕变本构模型，获得偏应变张量和球应变张量包括(11)式，具体如下：

[0040]

[0041] 其中， 为有效应力偏张量函数， 为有效应力球张量函数，eij()为偏应变函数， 为球应变张量函数，J()为推进剂蠕变柔量函数，ν是推进剂泊松比；

[0042] 所述推进剂蠕变柔量函数包括(12)式，具体如下：

[0043]

[0044] 其中，J0为初始蠕变柔量，NJ为蠕变柔量级数的项数，t为加载时间，Jn为第n项蠕变柔量， 为第n项蠕变柔量系数；

[0045] 所述有效应力偏张量函数包括(13)式，具体如下：

[0046]

[0047] 其中，Sij()为偏应力张量函数；

[0048] 所述有效应力球张量函数包括(14)式，具体如下：

[0049]

[0050] 其中，σkk()为球应力张量函数。

[0051] 进一步地，所述考虑损伤的推进剂蠕变本构模型的增量形式包括(15)式，具体如下：

[0052]

[0053] 其中，△εij(tm+1)为tm+1时刻应变张量的增量，△eij(tm+1)为tm+1时刻应变偏量的增量， 为tm+1时刻球应变张量的增量，δij为第四克罗内克符号。

[0054] 进一步地，所述一致切线刚度阵包括(16)‑(17)式，具体如下：

[0055]

[0056]

[0057] 其中，Cijkl()为切线刚度张量函数，△σij(tm+1)为tm+1时刻应力增量，△εkl(tm+1)为tm+1时刻应变张量的增量，C2222(tm+1)为tm+1时刻切线刚度张量的第一分量，C2233(tm+1)为tm+1时刻切线刚度张量的第二分量，C2323(tm+1)为tm+1时刻切线刚度张量的第三分量，C3333(tm+1)为tm+1时刻切线刚度张量的第四分量，C1133(tm+1)为tm+1时刻切线刚度张量的第五分量，C1313(tm+1)为tm+1时刻切线刚度张量的第六分量，C1111(tm+1)为tm+1时刻切线刚度张量的第七分量，C1122(tm+1)为tm+1时刻切线刚度张量的第八分量，C1212(tm+1)为tm+1时刻切线刚度张量的第九分量；

[0058] 所述切线刚度张量函数tm+1时刻切线刚度张量的所有分量中，除了第一分量至第九分量，其他分量均为零。

[0059] 与现有技术相比，本发明所达到的有益效果：

[0060] 本发明通过采用损伤发展方程客观表达考虑含损伤的推进剂粘弹性问题，将损伤发展方程融入构建的推进剂蠕变型本构模型，构建含损伤的推进剂蠕变型本构模型，不仅可以反映推进剂药柱结构的蠕变变形规律，还可以准确地计算推进剂药柱结构蠕变变形阶段的损伤效应。

[0061] 本发明通过分解、离散含损伤的推进剂蠕变型本构模型，推导考虑损伤的推进剂蠕变本构模型的增量形式和一致切线刚度阵，并将增量形式和一致切线刚度阵进行编程得到应用程序，使用有限元软件调用所述应用程序，能够为贮存阶段远程火箭弹药柱结构的精细化结构完整性分析提供理论基础和实现手段。

附图说明

[0062] 图1所示为本发明含损伤的推进剂蠕变型本构模型的构建方法的一种实施例流程图；

[0063] 图2所示为本发明含损伤的推进剂蠕变型本构模型的有限元应用方法的一种实施例流程图。

具体实施方式

[0064] 下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

[0065] 实施例1

[0066] 针对现有的推进剂本构模型不能精确预示蠕变损伤行为的问题，本实施例提供一种考虑损伤的推进剂蠕变本构模型的构建方法，参考图1，包括以下步骤：

[0067] 获取推进剂蠕变参数；

[0068] 预设推进剂有效承载面积相对变化量；

[0069] 根据获取的蠕变参数构建推进剂蠕变型本构模型；

[0070] 根据有效承载面积相对变化量确定损伤变量；

[0071] 根据确定的损伤变量确定损伤发展方程；

[0072] 根据确定的损伤发展方程和构建的推进剂蠕变型本构模型，构建考虑损伤的推进剂蠕变本构模型。

[0073] 本发明通过采用损伤发展方程客观表达考虑含损伤的推进剂粘弹性问题，将损伤发展方程融入构建的推进剂蠕变型本构模型，获得构建的含损伤的推进剂蠕变型本构模型，不仅可以反映推进剂药柱结构的蠕变变形规律，还可以准确地计算推进剂药柱结构蠕变变形阶段的损伤效应。

[0074] 实施例2

[0075] 在实施例1的基础上，本实施例介绍了：推进剂蠕变型本构模型构建方法、确定损伤变量的方法、损伤发展方程确定方法以及考虑损伤的推进剂蠕变本构模型。

[0076] (一)构建推进剂蠕变型本构模型

[0077] 根据获取的蠕变参数构建推进剂蠕变型本构模型，包括以下步骤：

[0078] S11通过(111)式获得推进剂蠕变张量：

[0079] Jijkl(t)＝[δkiδjl+v(δkiδjl‑δijδkl)]J(t)   (111)

[0080] 其中，Jijkl()为蠕变张量函数，δki为第一克罗内克符号，δjl为第二克罗内克符号，δkl为第三克罗内克符号，δij为第四克罗内克符号，J(t)为推进剂蠕变柔量函数，v为泊松比；

[0081] S12通过(112)式获得推进剂机械应变：

[0082]

[0083] 其中，εij为机械应变，σkl为推进剂损伤前的机械应力，τ为第一被积函数，t为加载时间。

[0084] S13通过(1)‑(3)获得推进剂蠕变型本构模型：

[0085] 依据时温等效原理，即高聚物的同一力学松弛现象可以在较高的温度、较短的时间或较低的作用频率观察到，也可以在较低的温度下、较长时间内观察到，将(112)式变形成(1)‑(3)，获得推进剂蠕变型本构模型，具体如下：

[0086]

[0087]

[0088]

[0089] 其中，Jijkl()为蠕变张量函数，ξ为时间t的折算时间，ξ`为第一被积函数τ的折算时间，T()为温度，t`为第二被积函数，aT为温度移位因子，C1为推进剂WLF方程第一参数，C2为推进剂WLF方程第二参数，T为当前温度，Tr为参考温度。

[0090] (二)确定损伤变量

[0091] 对于单轴拉伸载荷下的杆结构，结构在没有损伤时的初始承载面积为A，在推进剂材料发生内部损伤，伴随着微裂纹的出现和扩展，推进剂材料发生损伤后的有效承载面积为 则有效承载面积变化为AD，具体如下：

[0092]

[0093] 基于损伤各向同性假设条件，则可以定义单一损伤变量D来描述材料的损伤状态。为便于描述，本发明采用损伤前后的有效承载面积的相对变化量为损伤变量D，通过(4)式确定损伤变量D：

[0094]

[0095] (三)确定损伤发展方程

[0096] 首先，根据(5)式推导推进剂损伤后的机械应力，具体如下：

[0097]

[0098] 其中，σ为推进剂损伤前的机械应力， 为推进剂损伤后的机械应力。

[0099] 然后，根据周建平提出的微裂纹动力学方法，推进剂损伤发展方程包括(6)‑(9)式，具体如下：

[0100]

[0101]

[0102]

[0103]

[0104] 其中，D0为推进剂初始损伤，D△V为扩展形成的新损伤，σth为损伤应力阈值，γ推进剂损伤第一参数，K推进剂损伤第二参数，β为推进剂损伤第三参数，h()为应力状态函数，x为应力状态函数的输入参数。

[0105] 值得注意的是，在拉压异性的状态下，(7)式中的σ应该采用推进剂当量应力σ*代*替，通过(311)式获得当量应力σ；

[0106]

[0107] 其中，σij为损伤前推进剂应力张量，σkk为损伤前推进剂体积应力，hD为微裂纹闭合参数，ν是推进剂泊松比，<·>为麦克劳林算子；

[0108] 麦克劳林算子满足(312)式，具体如下：

[0109]

[0110] (四)考虑损伤的推进剂蠕变本构模型

[0111] 构建考虑损伤的推进剂蠕变本构模型包括通过(10)式更新机械应变，具体如下：

[0112]

[0113] 其中， 为有效应力张量函数。

[0114] 实施例3

[0115] 本实施例提供一种考虑损伤的推进剂蠕变本构模型的有限元应用方法，参考图2，包括以下步骤：

[0116] 分解实施例1或2构建的考虑损伤的推进剂蠕变本构模型，获得偏应变张量和球应变张量；

[0117] 分别离散偏应变张量和球应变张量，并根据离散后的偏应变张量和球应变张量推导考虑损伤的推进剂蠕变本构模型的增量形式；

[0118] 根据考虑损伤的推进剂蠕变本构模型的增量形式推导有限元材料子程序二次开发需要的一致切线刚度阵；

[0119] 对增量形式和一致切线刚度阵进行编程，得到应用程序，并使用有限元软件调用所述应用程序。

[0120] 本发明通过分解含损伤的推进剂蠕变型本构模型获得本构方程，编程获得的增量形式和一致切线刚度阵，获得应用程序，并使用有限元软件调用所述应用程序，能够为贮存阶段远程火箭弹药柱结构的精细化结构完整性分析提供理论基础和实现手段。

[0121] 实施例4

[0122] 在实施例3的基础上，本实施例详细介绍了分解考虑损伤的推进剂蠕变本构模型的方法、推导考虑损伤的推进剂蠕变本构模型的增量形式的方法以及推导一致切线刚度阵的方法。

[0123] (一)分解考虑损伤的推进剂蠕变本构模型

[0124] 分解实施例1或2构建的考虑损伤的推进剂蠕变本构模型，获得偏应变张量和球应变张量，包括以下步骤：

[0125] S411通过(411)式将考虑损伤的推进剂蠕变本构模型的应变张量分解成偏应变张量函数和球应变张量函数，具体如下：

[0126]

[0127] 其中，εij()为应变张量函数，Sij()为偏应力张量函数，σij()为球应力张量函数，eij()为偏应变张量函数，Y1()为推进剂第一蠕变函数，Y2()为推进剂第二蠕变函数，为球应变张量函数；

[0128] 应用中，球应变张量满足(412)式，具体如下：

[0129]

[0130] 其中，αT为热胀系数，Θ为温度变化，εkk()为机械球应变张量函数。

[0131] 应用时，温度变化包括(413)式，具体如下：

[0132] Θ＝T‑T0   (413)，

[0133] 其中，T为推进剂当前温度，T0为推进剂初始温度。

[0134] S412通过(414)式剪切推进剂第一蠕变函数与推进剂蠕变第二函数获得考虑损伤的推进剂蠕变本构模型的蠕变模量和体积蠕变模量，具体如下：

[0135]

[0136] 其中，x()为考虑损伤的推进剂蠕变本构模型的蠕变模量函数，B()为体积蠕变模量函数；

[0137] S413通过(415)式剪切蠕变模量，通过(416)式剪切体积蠕变模量，具体如下：

[0138] χ(t)＝2(1+ν)J(t)                      (415)

[0139] B(t)＝3(1‑2ν)J(t)                      (416)

[0140] 其中，J()为推进剂蠕变柔量函数，ν是推进剂泊松比；

[0141] S414通过(417)式分解偏应变张量，通过418)式分解球应变张量，具体如下：

[0142]

[0143]

[0144] 其中， 为有效应力偏张量函数， 为有效应力球张量函数。

[0145] S415根据Stieltjes卷积积分定义，通过(11)式更新偏应变张量和球应变张量，具体如下：

[0146]

[0147] 应用中，推进剂蠕变柔量函数包括(12)式，具体如下：

[0148]

[0149] 其中，J0为初始蠕变柔量，NJ为蠕变柔量级数的项数，t为加载时间，Jn为第n项蠕变柔量， 为第n项蠕变柔量系数；

[0150] 应用时，有效应力偏张量函数包括(13)式，具体如下：

[0151]

[0152] 其中，Sij()为偏应力张量函数；

[0153] 本实施例应用时，有效应力球张量函数包括(14)式，具体如下：

[0154]

[0155] 其中，σkk()为球应力张量函数。

[0156] (二)推导考虑损伤的推进剂蠕变本构模型的增量形式

[0157] 分别离散偏应变张量和球应变张量，并根据离散后的偏应变张量和球应变张量推导考虑损伤的推进剂蠕变本构模型的增量形式。

[0158] 应用中，将分析时间[0,t]划分为[0,t1]，[t1,t2]，…，[tm‑1,tm]，…，[tM‑1,tM]共M个子时间增量步。针对含损伤的推进剂蠕变型本构方程，推导任意的时间增量步[tm,tm+1]内本构方程的增量形式包括(15)式，具体如下：

[0159]

[0160] 其中，△εij(tm+1)为tm+1时刻应变张量的增量，△eij(tm+1)为tm+1时刻应变偏量的增量， 为tm+1时刻球应变张量的增量。

[0161] (三)推导一致切线刚度阵

[0162] 根据考虑损伤的推进剂蠕变本构模型的增量形式推导有限元材料子程序二次开发需要的一致切线刚度阵，包括以下步骤：

[0163] S431通过(419)‑(420)式简化考虑损伤的推进剂蠕变本构模型的增量形式，具体如下：

[0164]

[0165]

[0166] 其中， 为tm+1时刻等效应力偏量的增量， 为tm+1时刻球应力张量J
的增量， 为本构第一中间变量，γ ()为本构第二中间变量函数， 为本构第三中间变量函数， 为本构第四中间变量函数， 为本构第五中间变量函数。

[0167] 应用中，通过(421)‑(422)式获得γJ，具体如下：

[0168]

[0169]

[0170] 通过(423)式获得

[0171]

[0172] 通过(424)式获得

[0173]

[0174] 其中， 为本构第六中间变量函数；

[0175] 通过(425)式获得

[0176]

[0177] S432通过简化的考虑损伤的推进剂蠕变本构模型的增量形式，推导可知应变增量包括(426)‑(428)式，具体如下：

[0178]

[0179]

[0180]

[0181] 其中，为本构第八中间变量；

[0182] S433根据应变增量，推导一致切线刚度阵，一致切线刚度阵包括(16)‑(17)式，具体如下：

[0183]

[0184]

[0185] 应用中，通过简化(429)式获得C1111(tm+1)，具体如下：

[0186]

[0187] 通过简化(430)式获得C1122(tm+1)，具体如下：

[0188]

[0189] 通过简化(431)式获得C1212(tm+1)，具体如下：

[0190]

[0191] 其中，△γ12(tk+1)为tm+1时刻12方向的工程应变张量增量，Cijkl()为切线刚度张量函数，△σij(tm+1)为tm+1时刻应力增量，△εkl(tm+1)为tm+1时刻应变张量的增量，C2222(tm+1)为tm+1时刻切线刚度张量的第一分量，C2233(tm+1)为tm+1时刻切线刚度张量的第二分量，C2323(tm+1)为tm+1时刻切线刚度张量的第三分量，C3333(tm+1)为tm+1时刻切线刚度张量的第四分量，C1133(tm+1)为tm+1时刻切线刚度张量的第五分量，C1313(tm+1)为tm+1时刻切线刚度张量的第六分量，C1111(tm+1)为tm+1时刻切线刚度张量的第七分量，C1122(tm+1)为tm+1时刻切线刚度张量的第八分量，C1212(tm+1)为tm+1时刻切线刚度张量的第九分量。

[0192] 此外，切线刚度张量函数tm+1时刻切线刚度张量的所有分量中，除了第一分量至第九分量，其他分量均为零。

[0193] 本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD‑ROM、光学存储器等)上实施的计算机程序产品的形式。

[0194] 本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

[0195] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

[0196] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

[0197] 以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。