一种基于ANSYS参数化设计语言的三维模型梯度有限元求解方法

摘要

本发明属于有限元模拟领域，并公开了一种基于ANSYS参数化设计语言的三维模型梯度有限元求解方法，该方法主要包括如下步骤：模型前处理，求解器设置并求解；提取目标变量和节点自然坐标；求解形函数导数矩阵；求解雅各比列矩阵；求解雅各比列矩阵逆矩阵；计算插值点变量梯度；单元梯度值积分求单元解。本发明方法避免了同类问题导出数据、二次求解、费时费力的缺点，实现了后处理过程中梯度值的快速求解，能够大幅降低计算时间。

说明书

技术领域

本发明属于有限元模拟领域，更具体地，涉及一种基于ANSYS参数化设计语言的三维模型梯度有限元求解方法。

背景技术

有限元模拟分析是科学研究和工程设计的一个重要手段，而在众多的有限元软件中，ANSYS的应用又极为广泛，其中的参数化设计语言(ANSYS Parametric DesignLanguage，APDL)为用户基于特定模型的进一步的开发提供了可能性。

有限元的后处理已经提供了丰富的的结果。然而，还是有非常多的数据是后处理没有直接提供的，部分结果可以通过后处理提供数据的简单运算得到，但是依然有一些结果需要通过复杂的有限元理论来求解。在一些特殊的模拟分析中，需要分析目标变量梯度，但是ANSYS的后处理只提供了温度梯度，部分研究人员选择将目标变量数据导出作为温度载荷进行二次运算，得到的温度梯度就是目标变量梯度。但这并未从根本上解决梯度的求解问题，而且在瞬态计算中，频繁的导出导入数据和二次运算非常耗时。

发明内容

针对现有ANSYS后处理中只有温度梯度结果，且常见梯度求解方法需要二次求解，耗时费力的问题，本发明提供了一种基于ANSYS参数化设计语言的三维模型梯度有限元求解方法，其目的在于采用ANSYS参数化设计语言实现后处理过程中一次性快速求解单元梯度值。

为实现上述目的，按照本发明，提供了一种基于ANSYS参数化设计语言的三维模型梯度有限元求解方法，其特征在于，该方法包括以下步骤：

(1)对三维模型进行前处理；

(2)获取求解单元各节点的目标变量数据A_i和自然坐标C_i(x,y,z)，根据各节点的自然坐标C_i(x,y,z)获得节点坐标矩阵(C₁,C₂,……C_n)，其中i为节点编号，n为节点数量，x,y,z为各节点自然坐标系下的坐标；

(3)根据求解单元特性选择形函数N(r,s,t)，将所有节点的局部坐标代入形函数N(r,s,t)，从而得到形函数矩阵，再对形函数矩阵求导，进而得到各节点的形函数导数列矩阵S_i(r,s,t)，所有节点的形函数导数矩阵S_i(r,s,t)构成形函数导数矩阵S(r,s,t)，其中r,s,t为局部坐标；

(4)将形函数导数矩阵和节点坐标矩阵的相乘，得到雅各比矩阵J(r,s,t)＝S(r,s,t)(C₁,C₂,C₃...,C_n)＝S(r,s,t)C(x,y,z)然后再代入插值点坐标，得到各插值点对应的雅各比列矩阵J_j，其中j为插值点编号并且插值点数量与节点数量相同；

(5)获取雅各比列矩阵的逆矩阵

(6)根据下式获得各插值点P_i(r,s,t)的梯度值：

(7)对求解单元积分，以获得单元解P_elem，计算方法如下：

其中H_i为插值点积分权重。

优选地，所述步骤(6)中插值点梯度值求解采用三层循环计算实现。

优选地，所述求解单元特性为力学特性、热学特性和/或电学特性。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1)本发明的方法，提供一种基于ANSYS参数化设计语言的可用于梯度计算的有限元方法，可以一次性在后处理过程中得到目标变量梯度值，不需要数据导出和二次模拟，能够大大提高计算速度。

2)本法明提供的方法，适宜于所有三维结构体单元的应力梯度求解，不同体单元可以结合单元自身特征求解梯度

附图说明

图1是本发明实施例中SOLID5单元的示意图；

图2是本发明实施例中SOLID5单元积分点位置示意图；

图3是本发明实施例中基于SOLID5单元的静水应力梯度求解流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

图1是本发明实施例中SOLID5单元的示意图。单元包括八个节点，编号分别为I,J,K,L,M,N,O,P。对应的局部坐标为(1,-1,-1),(1,1,-1),(-1,1,-1),(-1,-1,-1),(1,-1,1),(1,1,1),(-1,1,1),(-1,-1,1)。其形函数为：

其中i为节点编号。

图2是本发明实施例中SOLID5单元积分点位置的示意图。插值点编号1,2,3,4,5,6,7,8和单元节点I,J,K,L,M,N,O,P一一对应，根据SOLID5单元特性，插值点坐标为对应节点坐标0.577350269189626倍，例如插值点1的坐标为0.577350269189626×(1,-1,-1)，积分权重均为1，即H_i＝1。所以单元积分数值运算简化为：

图3是本发明实施例中基于SOLID5单元求解静水应力梯度的ANSYS参数化设计语言求解流程图。该方法的具体步骤包括：

(1)根据模拟需求进行前处理，对求解器设置，譬如对节点进行赋初值，设置模型的边界条件等。

(2)采用*GET命令提取八个节点的三个主应力值σ_xx,σ_yy,σ_zz，之后计算静水应力σ_H＝1/3×(σ_xx+σ_yy+σ_zz)，存储到变量NHYSTR中，提取所有节点的自然坐标(x_i,y_i,z_i)，存储在矩阵NCOD中；

(3)根据SOLID5单元特征，结合ANSYS中关于该单元的说明，选择形函数为

对每个节点对应的形函数求导，得到单个节点形函数导数，带入八个插值节点坐标，得到八个节点的形函数导数矩阵Dermat1，Dermat2，...,Dermat8,每个矩阵存储八个元素，分别为八个插值点带入后得到的数值；

(4)使用*MOPER命令和MULT命令分别将形函数导数矩阵Dermat1，Dermat2，...,Dermat8和节点坐标矩阵NCOD相乘，得到雅各比矩阵JacMat1,JacMat2,...,JacMat8；

(5)使用*MOPER命令和INVERT命令求解得到雅各比矩阵的逆矩阵Jacinv1,Jacinv2,...,Jacinv8；

(6)使用*DO，*ENDDO命令，将雅各比矩阵逆矩阵Jacinv1,Jacinv2,...,Jacinv8和形函数导数矩阵Dermat1，Dermat2，...,Dermat8,以及节点的静水应力NHYSTR对应元素相乘，计算出所有插值点的静水应力梯度分量SGX,SGY,SGZ，三个变量的平方根就是插值点的静水应力梯度SG，即：

(7)将求解所有插值点的结果相加就得到单元静水应力梯度SGELEM。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。