基于Hyperworks和LSDYNA的挤压仿真

利用有限元技术对动力电池包进行仿真分析主要可以做以下方面的工作：
（1）电池组热管理，可以建立虚拟的电池组和散热通道的三维模型，在此基础上分析散热效果并对不同方案进行对比和优化，取代了试验方法，大大提高了设计效率；



（2）电池的机械性能分析，仿真模拟冲击、碰撞，碾压，针刺对电池的影响；



（3）电池的电性能分析，可研究过充/过放，大电流，充/放，外部短路对电池的影响，也可研究汇流排、动力电缆的大电流发热和温升情况；



（4）电池的结构力学分析，可研究电池组的振动、耐久性和疲劳寿命。



在机械性能方面，挤压仿真是动力电池包必须通过的一项严苛的测试。本文就挤压仿真过程中使用的参数、卡片进行归纳总结。限于计算条件有限，仅以小模型验证仿真思路。

仿真所采用的模型如图所示，一刚性辊子以0~1000N的斜坡载荷挤压两端固定的简支梁。



所采用的仿真流程为:
网格划分—建立材料属性—壳单元属性—实体单元属性—接触—约束—载荷曲线—挤压力载荷—计算控制卡片—k文件的导出—导入ANSYS计算—后处理。
仿真效果如图所示。



该案例只能为电池包的挤压仿真提供思路，并不能代表真实的挤压仿真。实际上，辊子对电池包的挤压速度很缓慢，挤压的过程中可以看成无数个微小时间间隔的静态过程，因此电池包的挤压仿真用Abaqus做准静态仿真更准确。由于水平有限和硬件不足，仅仅以低速碰撞替代准静态挤压，为电池包的挤压仿真探索思路。



在本案例中用到的卡片和关键字总结如下。合理的运用这些关键字，可以对电池包的机械冲击、高速碰撞、自由跌落等进行仿真。

1 定义材料

1.1 MAT24号材料
用于定义弹塑性类型的材料，比如钢、铝等。通过以下参数定义其材料本构：
Rho—密度；E—杨氏模量；Nu—泊松比；SIGY—屈服强度；ETAN—切线模量；TOEL—删除网格单元的最小时间步长；C、P应变率系数（用以表征材料发生膨胀、收缩、剪切变形的速率）；LCSS：材料的应力应变曲线；LCSR—不同应变率下的屈服极限曲线；
注意：定义了LCSS应力应变曲线后，SIGY和ETAN参数将被忽略，因为应力应变曲线中可以推算出；
定义了LCSR后，C和P值也将忽略，因为LCSR指定的不同应变率系数曲线需要不同的C和P值去定义；
所以定义材料的两种方法为：
（1）通过LCSS输入材料真实的应变曲线；
（2）设定应变率系数C和P值。

1.2 MAT20号材料
是一种不发生变形的刚性材料。通过设定Rho、E、Nu三个参数即可完成设定。

2 定义属性

2.1 组件关键字*PART
每一个组件都是一个PART，如下图所示。



其中，EOSID为材料*EOS关键字的状态方程；
HGID：对该组件的沙漏控制，引用一个类型为HourGlass的卡片；
GRAV：重力初始化，0—所有部件的初始化；1—仅对当前部件初始化；
ADPOPT:网格自适应，0—无网格自适应；1—3D网格自适应；2—2D和3D网格均自适应；
TMID：*MAT_THEMAL定义的热力学属性。

2.2 壳单元关键字*Section_shell
该卡片主要用于定义壳单元的属性，比如厚度、积分算法等。



ELFORM：求解的积分算法；2—单点积分算法，计算快，不能精确的处理翘曲，不能在较粗的网格中使用；
SHRF：剪切因子，默认为1，推荐使用5/6。
T1、T2、T3、T4：四个节点处网格的厚度，往往定义T1就可以，T2、T3、T4就默认与T1相一致。
说明：对于壳单元，通常只需要修改积分算法ELFORM和厚度，其余的采用默认即可。

2.3 实体单元关键字* Section_solid
ELFORM：实体单元类型选型；AET为周围环境类型选项。对于碰撞、挤压、冲击类仿真，采取默认值即可。



2.4 刚性墙
刚性墙的建立有两种方法，一种是建立墙的壳单元组件，再对组件赋予MAT20号的刚体材料。另一种是通过Analysis面板，建立RWGeometric类型的刚性墙，如下图所示。



然后点击geom，通过指定一个节点和通过该节点的法向量的方式，定义刚性墙平面。



3 控制接触参数

接触在挤压、碰撞、冲击等仿真的直观重要，面面接触控制卡片如下图所示。



其中，SSID用于指定从面；MSID用于指定主面；SPR、MPR控制是否计算并输出界面力；0—不计算；1—计算，并在RCFORC中输出。
FS、FD分别指定静态和动态摩擦系数；
DC：计算摩擦力系数；
VC：粘性摩擦力系数；
VDC：粘性阻尼系数；
PENCHK：初始接触检查；0—不检查；1—检查；2—开启检查，执行接触的最小对角线搜索；
BT、DT：接触开始和结束的时间；
SFS、SFM：主面和从面的罚函数系数，默认值为1;
FSF:库伦摩擦系数；
VSF：粘性摩擦系数。

4 定义初速度

该卡片可用于设置物体碰撞的初速度、物体自由跌落的初速度等。其关键字卡片如下图所示。



其中，NSID用于通过节点集SET的方式为零件或蔽障施加初始速度；
NSIDEX通过节点集方式，指定不包含初速度的节点；
BOXID以BOX集合方式指定初始速度；
IRIGID为1或者-1时，重新导入重心位于BOXID内的刚性体惯量；为2时重新导入所有刚性体惯量；
VX、VY、VZ、VXR、VYR、VZR为分别施加X、Y、Z方向的平动速度和转动速度，速度的单位为mm/s，通过在速度值前边施加“-”号表示速度方向与坐标轴相反。
Option选型卡里设置不同SET卡片下的初始速度类型。
*intial_velocity_option
*intial_velocity 用于节点集（set_node_list）
*intial_velocity_generation用于组件集（set_part_list）
*intial_velocity_node(用于单个节点)

5 定义加速度

加速度可以作用于整个模型，也可以作用于单个零部件。加速度的定义通过Loadcols载荷集下的LoadBody卡片创建。如下图所示。



其中X、Y、Z、RX、RY、RZ分别表示施加X、Y、Z方向的加速度和角加速度；
LCID用于施加用户自定义的加速度曲线，因此可以用于冲击分析、或者是随机振动分析；
SF是加速度曲线的缩放系数，如果将其设为负值，则加速度方向与坐标轴相反。
PART用于指定加速度所施加在的局部零件。激活PART后，如下图所示。通过PSID选型卡，可以指定施加加速度的零部件。



6 定义挤压力



其中，NSID用于指定挤压力的作用点；DOFX用于指定挤压力施加的方向；LCID用于指定挤压力曲线；SF为挤压力的比例因子，值为-1时表示挤压力的方向与坐标轴正向相反。

7 定义加载曲线

曲线的定义在XYPLOT中，选择Curve Editor即可创建加载曲线。



8 输出参数控制

8.1接触控制（CONTROL_CONTACT）



其中，SLSFAC：接触刚度，默认为1；
RWPNAL:刚性墙惩罚因子。0—不考虑；1—考虑；一般默认设置为1。
ISLCHK：初始穿透检查。0和1—不检查穿透；2：检查穿透；
SHLTHK：STS（面面接触）和NTS（点面接触）中是否考虑壳单元厚度变化的影响；0—不考虑；1—除刚体外考虑；2—全部考虑；
PENOPT：对称刚度检查。当主面和从面之间刚度相差较大时，主面和从面之间的接触力可能相差悬殊，需要用此参数重新分配接触刚度。0或1—取接触主面中的最小刚度；2—取主面的刚度值；3—取从面的刚度值；
THKCHG—单面接触中，是否考虑接触过程中壳单元厚度的变化。
其余的参数，在整车碰撞、挤压、机械冲击仿真中使用较少，采取默认值或者不设置即可。

8.2能量控制（CONTROL_ENERGY）



其中，HGEN为是否计算沙漏能；1—不计算；2—计算；
RWEN为是否计算阻碍耗散能；1—不计算；2—计算；
SLNTEN为是否计算界面滑移能；1—不计算；2—计算；
RYLEN为是否计算Rayleigh能；1—不计算；2计算；

8.3输出控制（CONTROL_OUTPUT）



OPIFS:0.1;
IKEDIT:100;
IFLUSH:5000;
以上参数采取默认即可。

8.4壳单元控制（CONTROL_SHELL）



其中，WRPANG:控制壳单元的翘曲程度，翘曲度大于此值即告警；默认设为20。
ESORT:自动处理退化的单元，碰撞挤压过程中四边形单元退化为三角形单元；默认为1。
IRNXX：单元法线的更新。默认为-1。
ISTUPD:是否考虑单元厚度的变化。0—不考虑；1—考虑；
THEORY：壳单元的计算理论。默认为2。
BWC：默认为2。
MITER：默认为1。
PROJ：默认为0。

8.5实时体单元控制（CONTROL_SOLID）



ESORT：默认为0；
FMATRX：默认为1；
NIPTETS：默认为4；
SWLOCL:默认为2；

8.6计算终止时间控制（CONTROL_TERMINATION）



ENDTI：设置求解的结束时间。

8.7时间步长控制（CONTROL_TIMESTEP）



DTINIT:初始时间长值；
TSSFAC:计算时间步长比例；
ISDO：四节点壳单元的基本计算时间值；
TSLIMT:分配给壳单元最小的时间步长值；
DT2MS：与质量比例结果相关的时间步长值；板壳单元的最大稳定时间步长值计算公式为：



其中，Ls为板壳单元的特征长度；ρ为材料的密度；E为材料的弹性模量。

8.8 沙漏控制
沙漏是在单元发生的零能量变形现象。产生沙漏的单元发生了变形，但是单元中不存在应力和应变。导致沙漏的原因是：显示积分算法为了缩短计算时间，采用了单点单元缩减积分算法导致的。缩减积分算法是指参与积分的点数比节点数目少。例如以一个一阶四边形单元为离子，有1、2、3、4四个节点，但是仅仅有一个中心的节点P参与积分计算。假设该单元受到弯曲或者剪切载荷，会发生如下图所示的变形。在这三种变形当中，单元积分点上的主应力和剪应力状态都相同，因此该单元可以自由的在这三种变形当中转换，而不需要额外的外力作用，这就是所谓的沙漏的零能量现象，即不需要外力的变形。沙漏只存在于六面体单元和四边形单元，其他的三角形单元、梁单元、四面体单元不会引发沙漏模式，这是因为这些单元比较“强硬”，不会像四边形那样肆意的发生变形。如果采用全积分算法，则不会存在沙漏问题，但是会占据较长的计算时间。



当采用缩减积分的单元算法时，必然会引起沙漏现象，当总的沙漏能小于5%时，才认为计算的结果是可靠的。LS-DYNA中减小沙漏能的方法有：
（1）尺寸均匀的网格
（2）避免将集中力施加在单个节点上，因为这会激发沙漏模式。
（3）调整模型的体积粘性参数，一般LSDYAN会自动设定体积粘性参数，用户可以通过调整EDBVIS命令来设置该参数；
（4）小位移变形情况，使用EDHGLS增加沙漏系数（HGCO），以增加模型的刚度，抑制沙漏的变形；大位移变形情况时，沙漏系数HGCO超过0.15会造成不稳定；
使用CONTROL_HOURGLASS控制整体的沙漏，如下图所示。



其中，IHQ为沙漏的控制类型；常用设置有1—LSDYNA默认的标准模型；2—Flanagan-Belytschko积分类型；3—用于实体单元精确的体积积分的Flanagan-Belytschko积分类型；8适用于单元类型为16的全积分算法壳单元；当设置IHQ=8时，单元的翘曲刚度将会被计算，但是这会增加系统25%的计算时间。采取默认设置。

8.9 零部件的沙漏控制
当模型中分别采用实体单元和壳单元混合建模时，有的零部件的沙漏控制可能要严苛一些，这就需要对零部件单独进行沙漏的控制。在Proprty属性选项卡中建立类型为HourGlass的卡片。



激活HourGlass卡片，如下图。



其中，IHQ为沙漏的控制类型，共有1~8种类型。默认采用1，LSDYNA的标准算法；8—使用全积分对壳单元进行计算，可以精确的处理翘曲度，但是会增加25%的计算时间。
IQH和QM为全局刚度系数和局部刚度系数，默认为0.1；
IBH为体积粘性类型，默认为0；
Q1为二次体积粘性系数，默认为1.5；可用于减小沙漏；
Q2为线性体积粘性系数，默认为0.06；可用于减小沙漏；
QB和QW为壳单元翘曲和弯曲的沙漏系数，默认与QM相等；
注意：当零部件的沙漏过大时，可以考虑增大整体或者局部的刚度系数来抑制沙漏，刚度系数不能超过0.15；也可以调整Q1和Q2体积粘性系数来减小沙漏。
打开要控制沙漏的组件，激活HGID选项卡，输出上述所建立的HourGlass卡片。