CN102542086B - 自适应步长的顺序耦合分析方法 - Google Patents

Abstract

本发明公开一种自适应步长的顺序耦合分析方法。其步骤为：（1）建立耦合场分析模型，设定第一、二个子步的迭代步长，对各耦合状态变量进行第一、二个子步的顺序耦合分析，得到其迭代节点耦合状态变量值；（2）提取各耦合状态变量当前迭代节点及其前两个迭代节点的耦合状态变量值，得到各耦合状态变量曲线当前迭代节点的斜率；设定各耦合状态变量步长调整区间，得到下一子步迭代步长；（3）将所有耦合状态变量下一子步的迭代步长的最小值作为下一子步的系统迭代步长，并采用迭代步长对瞬态场耦合系统的各耦合状态变量进行下一子步的顺序耦合分析，得到下一子步的迭代节点的耦合状态变量值；返回执行步骤（2），直至达到系统分析终止时间。

Description

自适应步长的顺序耦合分析方法

技术领域

本发明涉及多场耦合系统仿真分析方法，尤其是涉及一种瞬态场耦合系统的自适应步长的顺序耦合分析方法。

背景技术

顺序耦合法是对两个或多个不同工程领域相互耦合的物理场按一定顺序迭代分析的方法，对难以建立包含所有耦合场自由度且耦合非线性程度不是很高的多场耦合问题具有良好的分析效果。采用顺序耦合方法处理多场耦合问题，能够降低分析模型维数，减轻系统方程病态，降低分析成本。

顺序耦合法分析稳态场耦合问题时，不涉及迭代步长问题，只需依次对各物理场进行分析，并将前一个分析结果作为载荷施加到第二个分析中，当耦合状态变量收敛时结束迭代，如稳态流固耦合、微机电系统(MEMS)开关的静力结构耦合等典型稳态场耦合问题都能通过顺序耦合分析得到各物理场平衡状态。而对于瞬态场耦合系统问题，由于物理场在时域上变化，解耦分析时，需对各场按瞬态进行分析，各物理场瞬态载荷即为耦合状态变量在迭代节点的解，瞬态步步长即为顺序耦合迭代步长。瞬态分析本质是对非线性载荷线性近似加载并分析的过程，对于耦合问题，因在步长区间内忽略了其它耦合场的作用，各瞬态步内存在非平衡载荷残差，过大的非平衡载荷残差将危及数值积分的收敛性。目前，针对该问题一般采用极小迭代步长来保证精度及稳定性，但极小步长必定造成过多的迭代次数，效率较低。

发明内容

本发明的目的在于提供一种自适应步长的顺序耦合分析方法。该方法通过建立基于耦合状态变量曲线斜率的步长自适应策略，根据非平衡载荷残差动态调整迭代步长，顺序分析各耦合场。

为了实现上述目的，本发明采用技术方案的如下：其自适应步长的顺序耦合分析方法包括如下步骤：

步骤（1）：建立运动场分析模型，然后进行初始条件分析，以便在初始系统模型中各物体的坐标与各种运动学约束之间达成协调，保证系统满足所有的约束条件；

设定一个初始迭代步长作为第一个和第二个子步的迭代步长，对运动场的各耦合状态变量进行第一个和第二个子步的顺序耦合分析，得到各耦合状态变量的第一个和第二个子步的迭代节点的耦合状态变量值；

步骤（2）：提取每个耦合状态变量的当前迭代节点及其前两个迭代节点的耦合状态变量值，利用公式（1）得到各耦合状态变量曲线的当前迭代节点的斜率：

式（1）中，i表示迭代节点的序号且i≥3，j表示耦合状态变量的编号，

分别为各个耦合状态变量在当前迭代节点及其前两个迭代节点的耦合状态变量值，λ_ij、λ_(i-1)j分别为各个耦合状态变量当前子步及其前一子步的迭代步长，L_n′(t_ij)为各个耦合状态变量的状态耦合变量曲线在当前迭代节点的斜率。

步骤（3）：设定每个耦合状态变量的步长调整区间，按以下方法得到各个耦合状态变量的下一子步的迭代步长：

若步骤（2）所述斜率L_n′(t_ij)的绝对值大于对应耦合状态变量的步长调整区间的设定值上限，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（2）：

λ_(i+1)j=λ_ij*τ                        （2）

式（2）中，λ_ij、λ_(i+1)j分别为当前子步与下一子步的迭代步长，τ为步长调整系数且0<τ<1；

若所述斜率L_n′(t_ij)的绝对值小于对应耦合状态变量的步长调整区间的设定值下限，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（3）：

λ_(i+1)j=λ_ij*υ                            （3）

式（3）中，υ为步长调整系数且υ>1；

若所述斜率L_n′(t_ij)的绝对值属于对应耦合状态变量的步长调整区间，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（4）：

λ_(i+1)j=λ_ij                          （4）

步骤（4）：将步骤（3）所获得的所有耦合状态变量的下一子步的迭代步长的最小值作为瞬态场耦合系统的下一子步的系统迭代步长，并采用所述系统迭代步长对瞬态场耦合系统的各耦合状态变量进行下一子步的顺序耦合分析，得到各耦合状态变量的下一子步的迭代节点的耦合状态变量值；

步骤（5）：返回执行步骤（2），直至达到所述瞬态场耦合系统的系统分析终止时间。

与现有技术相比，本发明具有的有益效果是：

1.通过提取耦合状态变量各迭代节点耦合状态变量值，代入构造的三点求导公式以获得到各耦合状态变量曲线的当前迭代节点的斜率，提供了高效跟踪耦合场耦合状态变量在多场耦合分析时变化进程的可行方法。

2.通过设定每个耦合状态变量的步长调整区间，对斜率L_n′(t_ij)的绝对值与对应耦合状态变量的步长调整区间的设定值进行动态判别，并通过步长调整系统动态调整每个耦合状态变量的下一子步的迭代步长。当斜率L_n′(t_ij)的绝对值大于设定值上限时，以小步长迭代分析，避免各瞬态步内积累过大的非平衡载荷残差，以保证多场耦合分析精度；当斜率L_n′(t_ij)的绝对值小于设定值下限时增加迭代步长，以提高多场耦合分析效率。

3.在每个迭代节点取所有耦合状态变量的下一子步的迭代步长的最小值作为瞬态场耦合系统的下一子步的系统迭代步长，对于包含多个耦合状态变量的复杂瞬态场耦合问题，有效避免各个耦合状态变量迭代步长因取值不一而造成各耦合场分析次序冲突，减少冗余的数据交互及迭代工作。

4.构建异构CAE系统动态双向融合联合仿真平台，采用自适应步长的顺序耦合方法对低压塑电器短路分断过程中的运动-电磁场耦合问题进行迭代分析，自适应步长顺序耦合分析方法与采用步长为0.125ms的固定顺序耦合方法精度相当，但电磁场迭代分析次数与运动场迭代分析次数都分别减少了11次，效率显著提高。采用步长较大的（0.5ms）的固定顺序耦合分析法对问题分析时，虽然分析效率也较高，但分析误差达到了21.9%，精度较差。这也表明，传统的固定步长顺序耦合分析方法难以在分析效率与分析精度之间取得平衡的效果，而本发明有效地弥补了这个缺点。

附图说明

图1是本发明在分析低压断路器短路分断过程中的运动场与电磁场耦合问题的流程总图。

图2是本发明异构CAE系统动态双向融合仿真分析平台工作原理图。

图3是本发明对低压断路器短路分断分析所得动触头位移曲线与采用固定步长顺序耦合分析方法结果对比图。

图4是本发明对低压断路器短路分断分析所得动触头速度曲线与采用固定步长顺序耦合分析方法结果对比图。

图5是本发明对低压断路器短路分断分析所得动触头加速度曲线与采用固定步长顺序耦合分析方法结果对比图。

具体实施方式

以下结合附图和实施例对本发明作进一步的说明。

如图1所示，是本发明在分析低压断路器短路分断过程中的运动场与电磁场耦合问题的流程总图。

低压塑壳断路器动、静触头间产生电动斥力，包括电流流经触头时收缩产生的Holm力及回路电流激发磁场产生的Lorentz力。Holm力在动静触头分离后即消失，Lorentz随磁场消失而消失。动触头臂在电动斥力作用下，克服触头弹簧预压力，先于操作机构动作，动静触头间产生电弧，电路电压急剧增大，达到限流目的，之后，电磁短延时控制器触发操作机构动作，带动动触头臂运动，电弧进入灭弧室直至消失。

分断过程各变量的耦合关系可划归为运动场与电磁场两类，因电场与磁场可以通过有限元分析工具直接耦合迭代，故将其划归为一个场。触头区域磁场由流经动、静触头臂电流产生的回路磁场及灭弧栅片被磁化所产生的磁化磁场叠加而成，回路磁场又包含动触头臂电流源磁场及静触头臂电流源磁场，电流源磁场分布受动触头臂位姿及短路电流强度影响，电流源磁场的变化又影响到磁化磁场的磁感应强度，进而影响整个触头区域磁场分布。另一方面，动触头臂所受等效洛伦兹力由触头区域磁场分布及动触头臂位姿决定，等效洛伦兹力的变动会对操作机构多体动力学运动产生影响。建立运动场与电磁场数学模型，构建自适应步长的顺序耦合分析方法，对短路分断中电磁-运动耦合问题迭代分析，分析步骤为：

本发明自适应步长的顺序耦合分析方法包括如下步骤：

步骤1.建立运动场分析模型：采用多体动力学分析软件建立操作机构样机模型，设定各元件质量及转动惯量，施加相关约束副及外力。系统自动组装质量矩阵、外力矢（矩）矩阵、运动约束矩阵等系数矩阵，建立运动分析代数微分方程组（DAE）。

步骤2.进行初始条件分析，以便在初始系统模型中各物体的坐标与各种运动学约束之间达成协调，保证系统满足所有的约束条件。

初始条件分析通过求解相应的位置、速度、加速度的目标函数的最小值得到：

Minimize： $C=\frac{1}{2}{(p-p_0)}^{T}W(p-p_0)$

Subject to：Φ(p)=0

其中，W为权重矩阵，p₀为用户输入的值。

选择积分算法（SI2/SI1/I3）及求解器（GSTIFF/WSTIFF/Dstiff/CONSTANT_BDF）,采用BDF刚性积分法对上述代数微分方程组求解，得到各个时刻相关元件在全局坐标系中位置、速度、加速度及元件间反力等动力学分析指标，提取动触头位姿矩阵。

步骤3.设定一个运动场初始迭代步长作为第一个和第二个子步的迭代步长，初始迭代步长应足够小，一般以满足各耦合场舍入误差为宜，该实例中取初始迭代步长为0.1ms。对运动场的各耦合状态变量进行第一个和第二个子步的顺序耦合分析，得到各耦合状态变量的第一个和第二个子步的迭代节点的耦合状态变量值。

步骤4.提取运动场每个耦合状态变量的当前迭代节点及其前两个迭代节点的耦合状态变量值，利用公式（1）得到各耦合状态变量曲线的当前迭代节点的斜率：

式（1）中，i表示迭代节点的序号，因为公式（1）需要当前迭代节点及其前两个迭代节点的耦合状态变量值，故i≥3，j表示耦合状态变量的编号，设第一个子步开始时刻为初始迭代节点，序号为0，之后各个子步的节点序号依次类推。

分别为各个耦合状态变量在当前迭代节点及其前两个迭代节点的耦合状态变量值，λ_ij、λ_(i-1)j分别为各个耦合状态变量当前子步及其前一子步的迭代步长，L_n′(t_ij)为各个耦合状态变量的状态耦合变量曲线在当前迭代节点的斜率。

步骤5.设定运动场每个耦合状态变量的步长调整区间，按以下方法得到各个耦合状态变量的下一子步的迭代步长：

若所述斜率L_n′(t_ij)的绝对值大于对应耦合状态变量的步长调整区间的设定值上限，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（2）：

λ_(i+1)j=λ_ij*τ                      （2）

式（2）中，λ_ij、λ_(i+1)j分别为当前子步与下一子步的迭代步长，τ为步长调整系数且0<τ<1；

若所述斜率L_n′(t_ij)的绝对值大于对应耦合状态变量的步长调整区间的设定值下限，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（3）：

λ_(i+1)j=λ_ij*υ                      （3）

式（3）中，υ为步长调整系数且υ>1；

若所述斜率L_n′(t_ij)的绝对值大于对应耦合状态变量的步长处在调整区间中，则每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（4）：

λ_(i+1)j=λ_ij                        （4）

步骤6.将步骤5所获得的运动场的所有耦合状态变量的下一子步的迭代步长的最小值作为瞬态场耦合系统的下一子步的系统迭代步长。

步骤7.建立电磁场分析模型，进行一个子步的电磁场分析，并提取耦合状态变量值：进行单元类型、材料属性设置、网格划分等相关前处理工作。分析励磁电流。由于动、静触头臂形状不规则，需对动静触头臂进行电流传导分析。将电流远端电压设置为零，并耦合电流流入端面电流(CURR)自由度，施加短路总电流，所得各单元电流体密度J_q（q∈N，为有限单元个数）作为磁场分析励磁电流。分析叠加磁场磁感应强度。将触头臂电流密度分析结果作为励磁电流重载，设定灭弧栅片B-H曲线，在动静触头臂及灭弧栅片周围建立磁导率为1的空气单元，综合考虑电流源磁场在空气中的发散，耦合分析电流源磁场与磁化磁场的叠加磁场磁感应强度B_q，由

求解动触头臂所受洛伦兹力，并通过力矩等效将动触头臂各体单元所受洛伦兹力之和等效于动触头中心处。

步骤8.采用步骤7所述系统迭代步长运动场系统的各耦合状态变量进行下一子步的顺序耦合分析，得到各耦合状态变量的下一子步的迭代节点的耦合状态变量值。

步骤9.返回执行步骤4，直至达到运动场与电磁场耦合系统的系统分析终止时间。

采用本发明自适应步长的顺序耦合分析方法和步长分别为0.125ms与0.5ms的固定步长顺序耦合法对低压塑壳断路器的运动场与电磁场耦合问题进行分析，得到表1所示的结果：

由表1断路器分断时间及分析效率可知，各方法对应动触头臂转动30°的分断时间分别为14.4ms、14.2ms、11.8ms，采用相同工况对实物样机进行实验，测得其触头臂转过30°耗时15.1ms，误差计算可得，自适应步长的顺序耦合方法和0.125ms的固定步长顺序耦合方法分析结构与实验结果较接近，各个耦合场的分析次数比较可得，自适应步长顺序耦合方法在分析精度相近的情况下，耦合场迭代次数显著少于步长为0.125ms的固定步长顺序耦合方法。

本发明相应的异构CAE系统动态双向融合仿真平台工作原理如图2所示。动态双向融合仿真平台通过重启动分析实现仿真过程连续性，每个子过程都需要传输仿真文件和数据文件，包含几何模型、网格模型、载荷工况、边界条件、脚本、计算结果、仿真报告等。通过二次开发编制本地的数据交互平台，直接对仿真文件M_i(t)、N_i(t)和数据文件V_i(t)、W_i(t)交换，提高仿真效率。

如图3、4、5所示，是本发明自适应步长的顺序耦合方法与固定步长的顺序分析方法分析所得动力学曲线结果对比图，采用步长为0.125ms和0.5ms的固定步长顺序耦合分析方法结果作为对比。图3、4、5分别为自适应步长顺序耦合法、步长为0.125ms的固定步长顺序耦合法、步长为0.5ms的固定步长顺序耦合法对断路器短路分断分析的位移、速度及加速度曲线。采用相同工况对实物样机进行实验，自适应步长的顺序耦合方法基本能很好地反应HSM1-160在20ka短路电流下的实际分断情况。而步长为0.5ms的固定步长顺序耦合法由于迭代步长过大，与实验结果存在较大误差。自适应步长顺序耦合分析方法与采用步长为0.125ms的固定顺序耦合方法精度相当，但自适应步长方法通过动态调整迭代步长但电磁场迭代分析次数与运动场迭代分析次数都分别减少了11次，效率显著提高。采用步长为0.5ms的固定顺序耦合法分析对问题分析时，因为迭代步长取得过大，远远超过了非平衡载荷产生临界步长，导致各瞬态步非平衡载荷残差过大，虽然分析耗时较短，但精度较差，这也表明，传统的固定步长顺序耦合分析方法难以在分析效率与分析精度之间取得平衡的效果，而本发明自适应步长的顺序耦合分析方法有效地弥补了这个缺点。

Claims (1)

1.一种自适应步长的顺序耦合分析方法，其特征在于，包括如下步骤： 

步骤（1）：建立运动场分析模型，然后进行初始条件分析，以便在初始系统模型中各物体的坐标与各种运动学约束之间达成协调，保证系统满足所有的约束条件； 

设定一个初始迭代步长作为第一个和第二个子步的迭代步长，对运动场的各耦合状态变量进行第一个和第二个子步的顺序耦合分析，得到各耦合状态变量的第一个和第二个子步的迭代节点的耦合状态变量值； 

步骤（2）：提取每个耦合状态变量的当前迭代节点及其前两个迭代节点的耦合状态变量值，利用公式（1）得到各耦合状态变量曲线的当前迭代节点的斜率： 

式（1）中，i表示迭代节点的序号且i≥3，j表示耦合状态变量的编号，

分别为各个耦合状态变量在当前迭代节点及其前两个迭代节点的耦合状态变量值，λ_ij、λ_(i-1)j分别为各个耦合状态变量当前子步及其前一子步的迭代步长，L_n′(t_ij)为各个耦合状态变量的状态耦合变量曲线在当前迭代节点的斜率;

步骤（3）：设定每个耦合状态变量的步长调整区间，按以下方法得到各个耦合状态变量的下一子步的迭代步长： 

若步骤（2）所述斜率L_n′(t_ij)的绝对值大于对应耦合状态变量的步长调整区间的设定值上限，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（2）： 

λ_(i+1)j=λ_ij*τ              （2） 

式（2）中，λ_ij、λ_(i+1)j分别为当前子步与下一子步的迭代步长，τ为步长调整系数且0<τ<1； 

若所述斜率L_n′(t_ij)的绝对值小于对应耦合状态变量的步长调整区间的设定值下限，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（3）： 

λ_(i+1)j=λ_ij*υ                        （3） 

式（3）中，υ为步长调整系数且υ>1； 

若所述斜率L_n′(t_ij)的绝对值属于对应耦合状态变量的步长调整区间，则使每个耦合状态变量的下一子步的迭代步长与当前子步的迭代步长满足如下关系式（4）： 

λ_(i+1)j=λ_ij                      （4） 

步骤（4）：将步骤（3）所获得的所有耦合状态变量的下一子步的迭代步长的最小值作为瞬态场耦合系统的下一子步的系统迭代步长，并采用所述系统迭代步长对瞬态场耦合系统的各耦合状态变量进行下一子步的顺序耦合分析，得到各耦合状态变量的下一子步的迭代节点的耦合状态变量值； 

步骤（5）：返回执行步骤（2），直至达到所述瞬态场耦合系统的系统分析终止时间。 

Images