仿真约束

2024-05-26

仿真约束（精选6篇）

仿真约束篇1

摘要：针对非对称网络路段容量约束交通均衡分配模型计算困难,设计了一种带路段容量约束的用户均衡交通分配仿真算法。在算法迭代过程中,将按全有全无法在当前最短路上分配流量与前一轮迭代所得到的流量加权组合,各O-D对的组合系数依Logit模型来确定;并不断自适应调节路段排队延误因子和误差因子来模拟实际路段行驶时间,使路段流量逐步低于路段容量,从而达到广义用户均衡,克服了容量约束均衡分配计算量大及Logit随机分配法要求枚举所有路径的困难。随后证明了算法的收敛性,并对一个小型路网进行了数值试验。

关键词：交通分配,容量约束,变分不等式,仿真算法

交通分配可分为均衡分配与随机均衡分配两类,Beckman最早研究了离散网络路段无容量约束的交通分配问题,提出了满足Wardrop用户均衡的分配模型,Frank-Wolfe算法是求解这个模型的最有效方法。但这个模型常使很多路段的流量异常高,这不符合实际情形,当对每条路段加一个容量约束时,又会导致Frank-Wolfe算法的优势丧失。许多学者采用内罚函数法、外罚函数法或广义拉格朗日乘子法将路段容量约束问题转化为一系列无容量约束子问题[1,2,3],通过对惩罚因子的不断调节,反复采用Frank-Wolfe算法求解,由于Frank-Wolfe算法收敛速度慢,妨碍了容量约束模型在大型路网上的应用。随后Dagazno和Sheffi给出了随机用户均衡(SUE)概念[4],此后许多学者对随机用户均衡问题也进行了深入研究[5,6,7,8],在随机均衡交通分配模型中,Logit分配模型是最重要的模型之一,但该分配模型需要首先对各O-D对之间的所有路径进行枚举,这同样导致了其在复杂大型路网上应用的困难。文献[9]提出一种拟Frank-Wolfe迭代算法,提高了分配效率。以上算法主要针对离散网络,而对于非对称网络中的交通分配问题,对应算法却比较少见[8]。

首先通过图表分析了实际路段行驶时间与无路段容量约束的函数近似行驶时间的差异,指出后者低估了实际路段行驶时间,从而可以通过不断校正排队延误因子和误差因子来模拟实际路段行驶时间。基于仿真思想,对非对称网络,提出了一种带路段容易约束的用户均衡交通分配仿真算法,它基于O-D对间最短路的动态生成,不需要路径枚举,在每轮迭代中,将按全有全无方法进行流量分配,并与前一轮迭代所得到的流量进行加权得到新一轮的路径流量,而各O-D对的加权系数则依据Logit模型来确定,通过不断自适应调节排队延误因子和误差因子,使路段流量逐步低于路段容量,最终达到广义用户均衡。一方面该算法不需要在每轮迭代中求解无约束交通分配模型,另一方面,它又避免了随机均衡分配法在各O-D对之间进行路径枚举。随后证明了算法的收敛性。通过一个交通分配问题算例,说明了算法的可行性。

1容量约束交通均衡分配问题

给定一个交通网络G=(N,A),N为节点集,A为路段集,R为起点集,S为终点集,R与S可以有公共元素;网络中共有|A|条路段;(r,s)是以r为起点,s为终点的O-D对;Prs为O-D对(r,s)间所有路径集合。容量约束交通均衡问题考虑了路段的容量,使该问题更加符合实际路网情形,在非对称网络,容量约束交通均衡模型具有如下变分不等式形式:

t(x*)T(x-x*)≥0; ∀x∈Ω (1)

式(1)中:

式(2)中:x*为均衡路段流量解;f $_{p}^{r s}$ 是O-D对(r,s)间路径p上的流量,组成的向量为f=(…,f,…)T;xa为路段a的流量,组成的向量为x=(…,xa,…)T;ta(x)为路段a的行驶时间函数,t(x)=(…,ta(x),…)T为路段行驶时间函数向量,假定是严格单调可微且雅克比矩阵∇t(x)是严格对角占优的;Ca为路段a的容量,C=(…,Ca,…)T;qrs是O-D对(r,s)间的需求量;如果O-D对(r,s)间路径p经过路段a则δ=1,否则为0。上述模型式(1)—式(2)当可行域Ω非空(Ω≠∅)时,具有唯一路段流量解[10]。

模型式(1)—式(2)的均衡条件为:

${\begin{cases} \sum_{a \in A} \tilde{t}_{a} (x^{*}) δ_{a p}^{r s} - u_{r s}) f_{p}^{r s}^{*} = 0, \\ \sum_{a \in A} \tilde{t}_{a} (x^{*}) δ_{a p}^{r s} - u_{r s} \geq 0, f_{p}^{r s}^{*} \geq 0 \\ d_{a} (x_{a}^{*} - C_{a}) = 0 ； x_{a}^{*} \leq C_{a}, d_{a} \geq 0 \end{cases} (3)$

式(3)中: $\tilde{t}_{a} (x) = t_{a} (x) + d_{a}$ 称为广义路段行驶时间,Lagrange乘子urs,da分别代表O-D对(r,s)间的最小行驶时间与路段a的排队延误时间,因此上述条件经常称为广义用户均衡条件。由于在拥挤网络中,da并不是可以忽略的[10],如图1(a)所示,在容量Ca附近时间函数近似曲线与实际时间曲线误差较大;从图1(b)可以看出,在无路段容量约束模型中,由于t(x)是严格对角占优的,当xa>Ca时,实际延误da要高于函数近似延误d′a,其误差因子Δda=da-d′a>0,因此仅用t(x)并不能模拟路段的实际行驶时间。需要注意的是,虽然在实际路网均衡条件下,排队延误是唯一的,但Lagrange乘子da并不唯一。

2路段容量约束用户均衡交通分配仿真算法

模型式(1)—式(2)是一个约束变分不等式问题,因具有容量约束,相应算法比较少见,因此现设计仿真算法来求解该模型。

2.1算法的迭代步骤

STEP0:输入O-D表及网络信息;设O-D对个数为N,令i表示第i 个O-D对;置各路段初始流量为0,置迭代次数k=1,第i个O-D对最短路径集Ai=∅;取允许误差ε>0,在0流下求各O-D对最短路p,(i=1,2,…,N),按全有全无分配各O-D对流量并计算路段行驶时间tk,令,随机产生Δd>0,∀a∈A;将产生的最短路归入各O-D对的最短路径集A(i=1,2,…,N)中,置mpki=1。

STEP1:确定各路段行驶时间t $_{a}^{k + 1}$ =t $_{a}^{k}$ +d $_{a}^{k}$ ,计算各个O-D对最短路p $_{i}^{k + 1}$ ,(i=1,2,…,N)。如果p∉Ai,置Ai:=Ai∪p $_{i}^{k + 1}$ ,路径p的出现次数mpk+1i=1;否则Ai中路径不发生变化,置mpk+1i=mpki+1。将每个O-D对的需求量按全有全无分配法全部分配到这批新产生的路径上,记所得到的各路径上的流量向量为。按一定规则选取流量加权组合系数α(0≤α≤1),对第i个O-D对,在其最短路径集Ai中分配流量为:。

STEP2:计算输出路段流量xk+1及路行驶时间t $_{a}^{k + 1}$ ,如果 $\max_{a \in A} | x_{a}^{k + 1} - x_{a}^{k} | + \max_{a \in A} {0, x_{a}^{k + 1} - C_{a}} \leq ε$ ,停止;否则转STEP3。

STEP3:令延误因子:

;误差因子:,置k:=k+1,转STEP1。

2.2α $_{i}^{k}$ 的选取方法

由于现实路网中各个O-D对的路径集和交通需求都存在一定差异,不同于Frank-Wolfe算法中各O-D对共用一个加权系数,本算法在迭代过程中各O-D对可选取不同的加权系数。选取方式如下:

计算第k次迭代时新产生的路径p $_{i}^{k}$ 的初始行驶时间tpki及A $_{i}^{k}$ 中各条路径的初始行驶时间,取:

$α_{i}^{k} = \max {\frac{\exp (- θ t_{p_{i}^{k}})}{\sum_{q_{i} \in A_{i}^{k}} m_{q_{i}} \exp (- θ t_{q_{i}})}, ζ_{i}} (4)$

式(4)中,θ>0为Logit模型中的分布参数;A为第k次迭代后,已经求出的第i个O-D对的最短路集合;tqi是路径qi的初始行驶时间;mqi表示到目前为止,路径qi作为最短路而被重复选中的次数。式(4)中ζi是一个与算法迭代误差ε有关的充分小的正数,其作用是既要保证算法能够收敛,又要使每一步迭代路径流量都有适当的调整。

2.3算法的收敛性

定理:当Ω≠∅且qi<∞(i=1,2,…,N)时,则对任意给定的迭代误差ε>0,适当选取ζi(i=1,2,…,N),上述迭代算法是有限终止的。

证明:由于路网中只有有限条路段,要证明上述定理,只需证明对路网中任意一条路段a,有

$| x_{a}^{k} - x_{a}^{k - 1} | \leq \frac{ε}{2}, \max {0, x_{a}^{k} - C_{a}} \leq \frac{ε}{2}, k \to \infty$ 即可。

先证。

考察第i个O-D对,取:

及

其中P为各O-D对所有可能路径的集合。有:

$\frac{\exp (- θ t_{p_{i}^{k}})}{\sum_{q_{i} \in A_{i}} m_{q_{i}} \exp (- θ t_{q_{i}})} \leq \frac{Μ_{i}^{k}}{\sum_{q_{i} \in A_{i}} m_{q_{i}} m_{i}^{k}} = \frac{Μ_{i}^{k}}{m_{i}^{k}} \frac{1}{k} \leq \frac{Μ}{m} \frac{1}{k} \to 0, k \to \infty 。$ 可以看出,若取 $ζ_{i} = \frac{ε}{2 Ν | A_{i} | q_{i}}$ ,那么当k充分大时,总有α $_{i}^{k}$ =ζi,因此对A $_{i}^{}$ 中任一条路径pi,由

知 $| f_{p_{i}}^{k} - f_{p_{i}}^{k - 1} | = | α_{i}^{k} | | \tilde{f}_{p_{i}}^{k} - f_{p_{i}}^{k - 1} | \leq q_{i} α_{i}^{k} = q_{i} ζ_{i}, k \to \infty 。$

其中:qi第i个O-D对的交通需求量。由于交通网络中每个O-D对路径是有限的,即当k→∞时Ai中路径不再产生变化,从而对路网中任意一条路段a,其路段流量为 $x_{a}^{k} = \sum_{i = 1, 2, \dots Ν} \sum_{p \in A_{i}} δ_{a p} f_{p}^{k}$ ,即有

$\begin{array}{l} | x_{a}^{k} - x_{a}^{k - 1} | \leq \sum_{i = 1, 2, \dots Ν} \sum_{p \in A_{i}} δ_{a p} | f_{p}^{k} - f_{p}^{k - 1} | \leq \\ \sum_{i = 1, 2, \dots Ν} | A_{i} | q_{i} ζ_{i} \leq \frac{ε}{2}; k \to \infty 。 \end{array}$

再证。

$\begin{array}{l} f_{p}^{k} = α_{i}^{k} \tilde{f}_{p}^{k} + (1 - α_{i}^{k}) f_{p}^{k - 1} = (1 - α_{i}^{k}) f_{p}^{k - 1} = \\ (1 - ζ_{i}) f_{p}^{k - 1} \to 0, k \to \infty 。 \end{array}$

既有 $x_{a}^{k} = \sum_{p \in Ρ} δ_{a p} f_{p}^{k} \to 0 ‚ k \to \infty$ 。与 $x_{a}^{k} - C_{a} > \frac{ε}{2} > 0$ , 矛盾。

本算法是通过最短路算法来产生新路径,且总是在最短路上进行流量分配,而一些较差路径不能被选择或分配交通流量,这符合用户均衡交通分配的原则。随着迭代次数的增加,算法将不会产生新路径,并且通过流量的不断调节,最终达到广义的用户均衡。在STEP1与STEP3中,路段行驶时间采用仿真形式得到,因此该算法为仿真算法。在STEP3中,误差因子更新方式和Lagrange乘子更新方式相似,但算法不需要在每次迭代中进行无约束交通分配运算,并且每次迭代都更新已产生的所有路径流量,更重要的是算法可以求解非对称网络。由于算法虽然不能直接得到准确延误因子,但却是非对称网络交通分配的快速算法。由均衡条件式(3)可知延误因子可在多面体式(5)中得到。

${\begin{matrix} f_{p}^{r s}^{*} \sum_{a \in A} t_{a} (x^{*}) δ_{a p}^{r s} + f_{p}^{r s}^{*} \sum_{a \in A} d_{a} δ_{a p}^{r s} = \\ u_{r s} f_{p}^{r s}^{*}, & \forall r, s, p \\ \sum_{a \in A} t_{a} (x^{*}) δ_{a p}^{r s} + \sum_{a \in A} d_{a} δ_{a p}^{r s} \geq u_{r s}, & \forall r, s, p \\ d_{a} (x_{a}^{*} - C_{a}) = 0, & \forall a \\ d_{a} \geq 0, & \forall a \end{matrix} (5)$

虽然多面体式(5)中只含有一组变量da,∀a,但da并不唯一,这与Lagrange乘子的不唯一性是一致的。

3计算实例

下面通过一个实例来检验本文提出的分配算法,采用文献[12]的初始数据和道路交通网络:

在如图2所示的网络中有4个O-D对:1-12、 9-4、12-1、4-9,O-D需求量分别为q1,12=460、q9,4=400、q12,1=440、q4,9=420,网络共有34条路段; 路段旁数字表示路段流量为0时的行驶时间, 设每条路段的容量均为300,采用非对称路段行驶时间函数: $t_{a} (x) = t_{a} (0) [1 + 0.5 (\frac{x_{a} + 0.2 x_{\hat{a}}}{1.2 C_{a}})^{4}]$ 。 Ca为路段a的容量, $x_{\hat{a}}$ 为路段a的反向路段 $\hat{a}$ 的流量。

x $_{a}^{C}$ 表示算法所得到的路段流量, 表1给出了相应的分配结果(取ε=1, 保留小数点后两位)。

从表1可以发现,算法所得的结果中,最高路段流量仅为300.06,是允许范围内的误差,意味着该路段达到饱和。图3给出了算法的收敛曲线,可以看出本算法收敛的速度非常快;

图4、图5分别给出了路段(3,2)与路段(7,6)在算法迭代过程中流量的变化趋势,在100次迭代以后趋于稳定。

4结论

对非对称网络,提出了一种带路段容量约束的用户均衡交通分配仿真迭代算法,通过不断修正排队延误因子和误差因子来模拟实际路段行驶时间,再利用最短路算法产生有效路径,在当前迭代过程中,将全有全无分配法得到的各路径流量与上一轮迭代后各路径分配流量进行加权组合,得到新的路径分配流量,各O-D对加权组合系数通过Logit模型得到,随后给出了算法的基本思想与迭代步骤,并证明了算法的收敛性。计算实例还显示:该算法所得分配结果在误差范围内均不超过路段容量,收敛速度非常快,各路段流量在有限次迭代后便趋于稳定。由于分配流量的路径均是各O-D当前迭代下的最优路径,而未被选中的路径(相对较差的路径)不分配流量,这恰好符合用户均衡交通分

配原则,因此算法收敛时可达到广义用户均衡。本算法不需要路径枚举,而比求解容量约束下均衡交通分配模型计算量又小得多。因而它是一个有效的容量约束交通分配算法,适合在非对称道路交通网络中应用。

参考文献

[1] Larsson T,Patriksson M.An augmented lagrangian dual algorithm forlink capacity side constrained traffic assignment problems.Transpor-tation Research,1995;29:433—455

[2] Nie Y,Zhang H M.Models and algorithms for the traffic assignmentProblem with link capacity constraints.Transportation Research,2004;38:285—312

[3] Xu Y D,Li S J,Teo K L.Vector network equilibrium problems withcapacity constraints of arcs.Transportation Research Part E,2012;48:567—577

[4] Sheffi Y.Urban transportation networks:Equilibrium analysis withmathematical program methods.New Jersey,Prentice-Hall,INC,En-glewood Cliffs,1985

[5] Chiou Suh-wen.A novel algorithm for area traffic capacity control withelastic travel demands.Applied Mathematical Model,2011;35:650—666

[6] Meng Qiang,Khoo Hooiling.A computational model for the probit-based dynamic stochastic user optimal traffic assignment problem.JAdv Transp,2012;46:80—94

[7] Szeto W Y,Solayappan M,Yu Jiang.Reliability-Based transit as-signment for congested stochastic transit networks.Computer-AidedCivil and Infrastructure Engineering,2011;26:311—326

[8] D’Acierno L,Gallo M,Montella B.An ant colony optimisation algo-rithm for solving the asymmetric traffic assignment problem.EuropeanJournal of Operational Research,2012;217:459—469

[9]黄崇超,刘炳全.非平衡交通分配的Frank-Wolfe迭代算法.数学的实践与认识,2009;39(8):19—26

[10] Yang Hai,Huang Haijun.Mathematica and economic theory of roadpricing.Elsevier,2005

仿真约束篇2

对于一个并联机构，当其分支含有公共约束或机构中存在冗余约束时我们将其称为过约束并联机构，这类过约束也称为被动过约束[1]，本文提到的过约束机构指被动过约束并联机构。

当过约束并联机构中需要求解的未知力个数超过了可列出的独立平衡方程数时，该类机构的受力分析就非常困难，属于静不定问题。国内外关于过约束并联机构受力分析的文献相对比较少，董军等[2]对具有虚约束的平行四边形机构进行了受力分析，Huang等[3]基于反螺旋理论分别对含有过约束线矢力的4-R(CRR)(R和C分别表示转动副和圆柱副）对称四自由度并联机构进行了动力学分析，赵燕等[4]对含过约束力偶的4-UPU(U和P分别表示虎克铰和移动副）对称四自由度并联机构进行了动力学分析，文献[5,6]对每个支撑分支只提供单维轴向约束力的过约束并联机构进行了静力学分析。但是，这些关于过约束并联机构受力方面的研究文献大多是理论分析，其计算过程往往较为复杂，不易于快速掌握，限制了其应用。

虽然目前应用仿真软件对并联机构进行受力分析仿真验证或其他方面的研究已经不少，但对象大多是非过约束并联机构[7,8,9]，对此只需建立刚体模型。然而，对于过约束并联机构，在其刚体模型上添加运动副造成过约束后，ADAMS软件在仿真过程中会自动添加多余自由度将公共约束或冗余约束解除[10]，使其转化为非过约束机构。因此，非过约束并联机构的受力仿真研究方法已不适用于过约束并联机构。求解过约束并联机构的关节约束反力和驱动力需要结合变形协调方程，故应建立这类机构的刚柔混合模型对其进行受力仿真研究。目前已有的关于建立刚柔混合模型的研究仅针对一些非过约束机构[11,12]，或者是并非以过约束机构的受力分析为研究内容的机构学研究[13]，鲜有建立过约束并联机构刚柔混合模型并开展受力仿真分析的文献报道。

ANSYS软件虽然能建立机构的刚柔混合模型，但其静力分析模块中无法直接驱动刚柔混合模型，只能在确定位形下分析机构的受力情况，若分析该机构在其他位形的受力，则需要重建模型，使分析过程变得繁琐。ADAMS仿真软件不仅可以通过添加驱动使机构运动到任意指定位形，而且可直接测量关节约束反力在任意坐标轴上的分量[14]，对过约束并联机构的受力分析非常方便。因此，本文将联合三维建模软件、有限元分析软件ANSYS和多体动力学仿真软件ADAMS，提出一种建立过约束并联机构刚柔混合模型的方法，并对三种不同类型的过约束并联机构进行受力数值仿真，求解机构各关节约束反力和驱动力大小。

1 一种过约束并联机构受力的理论分析方法简介

为了便于验证本文提出的数值仿真方法的正确性，在此先引用一种过约束并联机构受力分析的理论方法[1]。

该理论方法的步骤如下：

(1）基于反螺旋理论分析各支撑分支提供给动平台的约束力螺旋系。

(2）依据力递推法，求解支撑分支在其约束力螺旋系作用下所有运动副的约束反力；利用材料力学基础知识，从支撑分支末端到基座依次求解各分支杆在运动副约束反力作用下引起的分支末端变形。

(3）将各分支杆引起的分支末端变形叠加即可得到分支末端的总弹性变形，分支约束力螺旋系幅值大小与该总弹性变形在分支各约束力螺旋轴线方向的投影大小之间的映射矩阵即为该支撑分支的约束力螺旋系刚度矩阵。

(4）结合各分支末端变形与动平台微位移之间的变形协调关系以及动平台的受力平衡方程即可求得机构的整体刚度矩阵。

(5）利用机构整体刚度矩阵和各分支约束力螺旋系刚度矩阵可求得各分支提供给动平台的约束力螺旋系幅值（包括驱动力）大小。

2 过约束并联机构受力分析的数值仿真方法

该数值仿真方法的一般步骤为：利用三维建模软件建立过约束并联机构的简易刚体模型；利用ANSYS软件建立所有分支杆的柔性体；利用ADAMS软件建立机构的刚柔混合模型并添加驱动及外力；进行静力学仿真并测量各关节约束反力；改变驱动值分析其他位形下机构的受力情况。其流程如图1所示。

2.1 建立三维模型

利用三维建模软件建立过约束并联机构的动平台、定平台以及所有分支杆模型，由于在整个受力分析过程中主要考虑的是分支杆的变形，所有运动副被视作刚体，又因在ADAMS仿真软件中用虚拟运动副代替真实运动副连接两构件并不影响机构运动和力的传递，即使不存在真实的运动副模型，亦能通过测量虚拟运动副得到所有关节力大小，因而本文在三维建模时不再建立外形复杂的运动副模型。零部件建模完成后按照D-H法规定的机构各运动副轴线位置、方向等条件进行装配，保证装配体具有机构正确的自由度。对于装配中出现的如图2所示的干涉现象无需处理，这种干涉在ADAMS软件中不影响过约束机构的静力分析结果。以3-RRC(R表示驱动加在此转动副处）三平移过约束并联机构[15]为例，其三维简易刚体模型如图3所示。最后将机构的装配体和所有分支杆保存为ANSYS软件和AD-AMS软件能识别的格式。

需要指出的是，后续仿真分析结果说明本文提出的方法在三维建模时对运动副的简化处理可行，它使得该方法建模过程简单快速，最终所建模型与一般理论分析模型一致，易于和理论方法之间互相验证。

2.2 建立支撑分支中各杆件的柔性体

按照ANSYS软件中生成有限元模型的步骤将需要柔性化的各分支杆进行网格划分，并进行模态计算，最终输出ADAMS软件能识别的模态中性文件。在此需要说明两点：(1)本文提出的数值仿真方法分析的是过约束并联机构在受到六维外力作用时，机构的驱动力和各关节约束反力大小，故建立柔性体分支杆时，给杆件赋予非常小的密度值，即忽略杆件的重力；(2)该方法中将过约束并联机构的动平台和定平台以及所有运动副均视作刚体，在柔性杆与柔性杆之间、柔性杆与动定平台之间添加运动副的地方需对柔性杆进行刚化处理，图4所示为对3-RRC机构的各分支杆两端需加转动副（圆柱副）的地方进行刚化处理的结果。

2.3 建立刚柔混合模型

首先，将三维建模软件保存的机构装配体导入ADAMS软件并添加虚拟运动副，此时完成的过约束并联机构模型为全刚性的，图5所示为3-RRC过约束并联机构的刚体模型。其次，用AN-SYS软件中生成的柔性杆替换刚体模型中对应的分支杆，并重新在柔性杆与柔性杆、柔性杆与动定平台之间添加虚拟关节。最后，对过约束并联机构的刚柔混合模型添加驱动与外载荷。图6所示为添加了虚拟关节、驱动及外载荷的3-RRC过约束并联机构的刚柔混合模型。

2.4 过约束并联机构关节约束反力/力矩和驱动力的测量

ADAMS软件中marker点起到局部坐标系的作用，为方便测量机构某方向的关节约束反力大小，先在该关节处添加一个marker坐标系oxyz，如图7所示。

静力学仿真结束后，通过图8所示的测量方法即可得到各关节约束反力（包括驱动力）在任意参考坐标系下的分量大小或总幅值大小。需要强调的是：分支驱动力与其他关节约束反力解耦时，沿驱动力方向的测量结果即为驱动力大小，反之，对测量结果进一步计算才能得知驱动力大小。例如3-RRC机构各分支提供给动平台2个约束力偶作用，分支中驱动力线矢沿靠近动平台的分支杆轴线方向，它与2个约束力偶解耦，故沿图7所示marker坐标系z轴方向的测量结果即为该分支驱动力线矢的大小。

2.5 任意位形下机构的约束力分析

在ADAMS软件中只需改变驱动副的输入量，驱动过约束并联机构的刚柔混合模型运动到其他位形，便可实现对过约束并联机构在任意位形下的静力分析。ANSYS软件的静力分析模块只能在已完成的装配位形下对刚柔混合模型进行受力分析，若要分析机构在其他位形下的受力情况还需重新调整原始装配体，该软件的静力分析模块中不能直接驱动刚柔混合模型运动，故本文提出的过约束并联机构受力的数值仿真方法采用在ADAMS软件中建立其刚柔混合模型的方法。例如，改变3-RRC过约束机构的驱动输入值，使动平台发生水平移动后其刚柔混合模型如图9所示，进而可分析此位形下3-RRC机构的驱动力和约束力大小。

3 基于数值仿真法的过约束并联机构受力分析实例

为说明采用该方法对过约束并联机构受力进行分析的通用性，选取三种不同类型的过约束并联机构进行数值仿真分析并应用前面引用的理论方法进行验证。这三种机构的不同体现在过约束情况不同，它们分别是各分支对动平台仅提供约束力偶的空间3-RRC机构、仅提供约束力的空间7-SS结构以及既提供约束力又提供约束力偶的平面五杆四边形机构。

3.1 3?RRC过约束并联机构的受力分析

空间3-RRC三平移并联机构如图10所示，该机构由三个相同的分支各通过一个C副和两个R副连接动、定平台组成，与定平台相连的转动副R为驱动副。每个分支的三个运动副轴线相互平行，且都平行于定平台B1B2B3，分支二和分支三的运动副轴线之间也相互平行。根据螺旋理论基本知识可判断出该机构每个分支提供给动平台2个垂直于分支运动副轴线的约束力偶，这6个约束力偶线性相关，其最大线性无关数为3，故该机构有3个冗余约束，属于过约束并联机构。根据文献[16]中驱动力螺旋的定义可知，每个分支还提供给动平台一个过分支中第二个转动副中心且沿靠近动平台的分支杆轴线方向的驱动力线矢，从物理意义上讲，该驱动力线矢就是由分支驱动副施加给动平台的1个约束力线矢，其大小与驱动副的驱动力矩大小之间满足一不定的关系[17]。因此，该机构的每个圆柱副共受到2个约束力偶和1个驱动力线矢作用。接下来用本文提出的数值仿真方法分析该机构的6个约束力偶和3个驱动力线矢的大小。

给定该机构的一组结构尺寸：6根圆柱杆横截面直径均为0.016m，杆长均为0.15m，图10中b1=0.15m,b2=0.02m，其中，b1、b2分别表示转动副中心B1、B2到定平台顶点A的距离。建立该机构的刚柔混合模型，如图5所示。建立图10所示的定坐标系OXYZ和动坐标系oxyz，取各分支转角α1=α2=α3=57°为该机构的初始位形。动坐标系下施加6维外力F=(10N,15N,-20N,9N·m,-10N·m,15N·m)T，静力学仿真结束后，定坐标系下测量可得到各分支提供的2个约束力偶大小，在图7所示marker坐标系z轴方向测量即可得到驱动力线矢大小。

为了验证仿真结果的正确性，用本文第1节中引用的理论分析方法[1]对3-RRC过约束并联机构进行受力计算，并将得到的理论值与仿真值进行误差分析。各分支提供给动平台的约束力偶和驱动力线矢幅值的仿真值、理论值及误差分析结果如表1所示，其中误差为仿真值与理论值的差值与理论值的比值。

注：表中数据仅表示各力大小，不包括方向。

从表1可以看出，本文提出的数值仿真方法对3-RRC过约束并联机构的受力分析结果与理论分析结果之间误差小于1%，充分验证了该方法的正确性。

本文提出的过约束并联机构受力的数值仿真分析方法不仅能获得机构处于某一位形时各力的大小，而且可以很方便地分析机构在运动过程中各驱动力和约束力的变化情况及具体大小。以3-RRC三平移过约束并联机构为例，在ADAMS仿真软件中调整该机构刚柔混合模型的三个驱动输入量，使动平台按照轨迹y(m)=0.01 t沿定坐标系Y轴运动2s（在这2s的运动过程中该机构未出现奇异位形），运动结束时刻机构所处位形如图9所示，根据本文提出的数值仿真方法即可获得在上述2s的运动过程中3-RRC机构各分支驱动力线矢和约束力偶幅值的变化情况及大小，其仿真结果如图11所示。

为了验证上述运动过程中得到的各力仿真曲线的正确性，可以针对3-RRC机构在整个运动过程中所处的一系列位形按照本文引用的理论方法[1]进行计算，得到理论曲线后与仿真曲线进行误差比较。这里我们给出t=2s时刻机构各分支的驱动力线矢和约束力偶幅值的仿真值与理论值及误差分析结果，如表2所示。

从表2可以看出，该位形下受力分析的仿真结果与理论计算结果之间误差也很小。

同理，可分析该机构处于任意位形时的受力情况。

3.2 7?SS过约束并联结构的受力分析

预紧式6-1/3-3-1型Stewart结构力传感器结构简图见图12，该机构由下平台、上平台和7个分支DiUi(i=1,2，…，7）组成，各分支均通过S副与上下平台连接，每个分支施加给上平台1个沿分支轴线方向的约束力，上平台共受到7个约束力作用，这7个约束力线性相关，其最大线性无关数为6，限制了上平台所有的自由度，该并联结构有1个冗余约束，属于过约束并联结构。

建立7-SS过约束6维力传感器并联结构的刚柔混合模型如图13所示，7个圆柱分支杆具有相同的横截面尺寸，除中间分支外其余6个分支杆长均为l。如图12所示，在上平台中心建立坐标系OXYZ,Ru为上平台球铰点的分布半径，Rd1、Rd2分别为下平台内圈和外圈球铰点的分布半径，H为上下平台间的距离，α1、β1分别表示下、上平台第一个球铰点和坐标原点连线在XY平面内的投影与坐标系X轴的夹角，α2、β2分别表示下、上平台第4个球铰点和坐标原点连线在XY平面内的投影与坐标系X轴的夹角。

根据文献[18]给定该机构的一组结构参数如表3所示。在坐标系OXYZ下给动平台施加6维外力F=(150 N,200 N,240 N,120 N·m,100N·m,160N·m)T，同样地，先按照上文介绍的理论方法[1]得到各分支提供给上平台的约束力的理论值，进而验证仿真结果的正确性。沿各分支轴线方向的约束力FDiUi(i=1,2，…，7）的仿真值与理论值如表4所示。

从表4可以看出，在六维外力F作用下，分支6提供给上平台的约束力与其他分支约束力相比非常小，其仿真值与理论值误差较大，其余分支误差均小于2%。

3.3 平面五杆平行四边形机构的受力分析

如图14所示的具有虚约束的平面五杆平行四边形机构也可看作是具有三个分支的并联机构，分支AB、FC、ED均通过转动副连接动平台BD杆和定平台AE杆。在A点建立图13所示的坐标系OXYZ，每个分支提供给BD杆1个过A点垂直于机构平面的约束力、1个X轴方向的约束力偶和1个Y方向的约束力偶。此外，每个分支还提供给BD杆1个沿分支轴线方向的约束力，这3个分支轴线方向的约束力线性相关，最大线性无关数为2，则该机构除具有3个公共约束外还具有1个冗余约束，故平面五杆平行四边形机构也属于过约束并联机构。若驱动加在转动副A处，则分支AB还提供了1个驱动力线矢作用。驱动转矩的大小通过驱动力线矢大小与其之间的关系计算可得。图14中α表示分支AB与定坐标系X轴之间的夹角。

令五根杆具有相同的横截面直径D=0.016m，三根分支杆长均为l=0.2m，动定平台等长，长度k=0.36m，建立该机构的刚柔混合模型如图15所示。

如图14所示，在BD杆质心建立动坐标系oxyz，动平台受到六维外力时，显然Fx、Fy、Mz方向的外力由该机构3个分支杆轴线方向的约束力和1个驱动力平衡，Fz、Mx、My方向的外力由该机构的3个公共约束平衡，故可分两种情况施加外力对该机构进行受力仿真分析：施加六维外力F=(0N,0N,12N,10N·m,8N·m,0N·m)T并对其进行静力学仿真，分支公共约束的仿真值与理论值如表5所示；施加六维外力F=(10N,8N,0N,0N·m,0N·m,5N·m)T并对其进行静力学仿真，分支AB提供的驱动力线矢大小和沿各分支杆轴线方向的约束力的仿真值与理论值如表6所示。

表5和表6的分析结果也表明了本文提出的过约束并联机构受力分析的数值仿真方法的正确性。上述分析结果对应的是该平面机构处于α=60°时的位形受力情况，同理，改变转动副A处的驱动值使机构运动到其他位形，可求解任意位形下该机构的驱动力和各关节约束反力大小。

4 结论

(1）提出了一种过约束并联机构受力的数值仿真分析方法，该方法结合三维建模软件、AN-SYS软件和ADAMS软件建立了过约束并联机构的刚柔混合模型，采用虚拟关节代替形状复杂的运动副模型，最终所建模型与理论分析模型完全一致。而且，该方法能够直接驱动模型使机构运动到任意位形对其进行受力分析，不需重建其模型，大大提高了工作效率。

(2）提出的仿真方法建模过程简单，不需要繁琐的理论计算即可快速得到过约束并联机构在任意位形下的驱动力和各关节约束反力大小，为过约束并联机构的受力分析提供了一种便捷途径。

(3）本文提出的方法不仅适用于上述三种类型的过约束并联机构，对其他类型的过约束并联机构也适用，通用性较强。由于该方法未考虑机构各杆件的重力和运动副的变形，后续研究将考虑这些因素。

摘要：鉴于目前对过约束并联机构进行受力仿真分析的研究非常少,提出了一种对该类机构受力的数值仿真分析方法。将机构的动平台、定平台及所有运动副视作刚体,将易产生空间复合弹性变形的分支杆视为柔性体,用虚拟关节代替外形复杂的运动副模型,结合三维建模软件、ANSYS软件及ADAMS软件建立与理论模型完全一致的刚柔混合模型;仅改变驱动量使机构运动到指定位形,便能实现对机构在任意位形下的受力分析,避免了对机构的不同位形重建模型。采用该方法对三种类型的过约束并联机构进行了受力分析,分析结果验证了该方法的正确性,同时也表明该方法具有一定的通用性,且无需复杂的理论计算,便能获得过约束并联机构在任意位形下的驱动力和各关节约束反力大小,为此类机构的受力分析提供了便捷途径。

仿真约束篇3

机械系统中,通过添加各种约束或运动副限制系统中活动部件的运动。为了更好地了解机械系统的运动特性,需要对机械系统进行动力学分析仿真。在进行动力学仿真时,需要建立相应的动力学模型。机械系统动力学模型主要有多刚体动力学模型、刚柔混合体动力学模型等。由于动力学模型是理想模型,有别于实际系统模型,所以,如果按照实际系统对理想模型进行约束时,通常会出现多余约束,以致分析仿真时得不到正确的结果甚至无法进行分析仿真。目前世界范围内最广泛使用的机械系统动力学分析仿真软件主要有Adams(Automatic Dynamic Analysis of Mechanical System)和DADS((Dynamic Analysis and Design System)等,它们都是基于欧拉-拉格朗日方程开发的常用动力学分析软件[1]。针对使用它们进行分析仿真时遇到的多余约束问题,本文分析了约束(即:约束方程)对求解欧拉-拉格朗日方程的影响,给出了基于欧拉-拉格朗日方程开发的动力学分析软件中多余约束的分析处理方法,并以Adams2005为例进行了工程验证。

1基于欧拉-拉格朗日方程的约束分析

欧拉-拉格朗日方程是动力学普遍方程基于广义坐标的表示形式。由于广义坐标数和多刚体系统的自由度数相等,并且广义坐标的选择可以视情况而定。因此,和动力学普遍方程相比,可以将动力学方程数目减为最少,使动力学方程组更为简单。其描述形式[2]如下:

式(1)中:M(q,t)是广义质量矩阵,是正定对称矩阵;Φ(q,t)是约束方程向量;Φq(q,t)是Φ(q,t)关于q的雅克比矩阵;F(q,6)q,t)是广义外力。

欧拉-拉格朗日方程有方程组有且仅有唯一解的充要条件是Φq(q,t)满秩,即:Φ(q,t)中的约束方程是独立的,不相关的。常用的代数约束方程见图1。例如:旋转副对应的代数约束方程为其中的1、2、3、4和5。

在理想模型中,Φq(q,t)经常出现不满秩的情况。如在运动的某个时段,由于难以确定约束方程组中约束方程的最小个数,方程组上经常出现多余的约束方程,造成Φq(q,t)不满秩。也就是说,多刚体系统的理想模型中通常会出现多余约束,造成求解欧拉-拉格朗日方程有方程组时得不到正确的结果甚至无法求解。随着科学技术的发展,出现了很多动力学数值计算方法来消除多余约束对动力学分析的影响,如:文献[2]中数值分析方法、文献[4]中的传统增广法和缩并法、文献[5]中的正交分解法和文献[6]中拉格朗日乘子法等。但是在使用动力学分析软件进行分析计算时,软件所采用的方法是固定的,不能人为进行更改,需要通过其它方式来解决多余约束问题

2约束处理的方法原则

在多刚体系统的物理模型中通常增加一些构件和运动副来改善机构的受力及力的传递情况。但动力学仿真时用的是理想模型这些构件和运动副通常是不需要的,属于多余构件和多余约束。因此在约束模型前,首先要对实际多刚体系统进行分析,找出虚约束和局部自由度并进行处理,处理方法如下:

(1)如果两构件的有多处相同的运动副连接并且运动副的方向相同,则在约束模型时只选择其中一个运动副进行连接。

(2)在机构中如果有两构件相联接,并当此两构件在该联接处拆开后,两构件上的联接点和轨迹是互相重合的,则该联接引入多余约束。那么在动力学分析软件中把引入多余约束的构件和运动副失效掉。

(3)在有些机构中,其某些构件所产生的局部运动并不影响其它构件的运动,那么该局部运动的自由度是多余[7]。可以去掉这些运动副和构件并用相应的运动副来约束该构件连接的两构件。

一般情况下,去掉虚约束和局部自由度后的多刚体系统仍然存在多余约束,即:系统的Φq(q,t)不满秩。因此,在用多刚体动力学分析软件进行动力学仿真分析时,需要变换其中一些运动副的类型或者高副低带来改变Φq(q,t)中的元素使Φq(q,t)满秩。如:如果旋转副的两个旋转约束是多余约束,那么可以球副代替旋转副;如滑动副的三个旋转约束是多余约束,那么可以把滑动副改为点线副;旋转副一个旋转方向是多余约束,那么可以用点线副和平行副代替。

综上所述,多余约束处理方法原则为:合理地简化约束模型、正确地变换约束类型和合理地进行高副低带。

3工程验证

以Adams2005为例,结合空间某装置的部分结构对本文所提方法进行了工程验证。。其中:Part2为机架,Part3只能沿Z轴方向移动,part3和part4以铰链连接,part4和part5及part6分别以铰链连接,part5和part2以铰链连接,part6和part7以铰链连接,part7和part＿9以两个铰链连接,两个part＿9分别在两侧part7沿part2导槽运动,part7和part＿8以两个铰链连接,两个part＿8分别在两侧part7沿part2导槽做直线运动。详见表1。

首先对模型进行分析: part7两侧的约束和构件对运动来说起相同的作用,在用模型中其中一侧的约束和构件是多余的,可以去掉。假定去掉part7负y一侧的约束和构件,在剩下一侧,part7和part_9间的运动不影响其它构件的运动,为局部自由度,所以将其去掉,把part7和part_9固接一起看为一个构件。

然后,以常用的动力学分析软件Adams为例说明本文所述方法的正确性。把三维模型导入Adams2005并根据上一步的分析对模型进行如下约束和处理:用旋转副代替铰链、用滑动副代替直线导轨、用点线副代替高副, part2为机架,用固定副把part2固定在地面(ground)上,另外,由于part8的作用是约束part7沿Z轴平动和沿Y轴转动,所以用点线副和平行副代替连接part8的导轨和铰链。约束完验证模型,验证结果为:系统的实际自由度为1,存在7个多余约束。多余约束如表2

正如前面所说,如果在存在多余约束的情况下,进行动力学仿真, 则会出现不正确的结果。比如,在正常情况下part2与part5间的旋转副JIONT_27的Y方向受力是连续的,而存在多余约束时却出现突变如图3所示。

JIONT_27的Y方向受力曲线说明多余约束(即:Φq(q,t)不满秩)会造成结果异常。为得到正确的结果,需要对Φ(q,t)进行变换,使Φq(q,t)满秩。由于Part2为机架,固定不动,所以不能对JIONT_1的代数约束方程(即:Φ(q,t)对应的元素)进行变换。也就是说不能按照表2中多余约束对Φ(q,t)中对应的元素进行变换。Φ(q,t)的变换方式有很多种,可以根据实际情况选择其中一种进行变换。结合多余约束的处理方法原则,本文给出如下变换:把Part3与part4间的旋转副、Part4与part5间的旋转副、Part6与part7间的旋转副换成点点副,把平行副换成垂直副。变换完成后验证模型,验证结果:系统的实际自由度为1,存在0个多余约束。

此时再次测量旋转副JIONT_27的Y方向受力,受力曲线见图4,从测量结果图看出,受力曲线是一条连续的曲线,此前曲线中存在的突变现在已经消除了。

最后,对该动力学分析软件给出的结果进行验证。按照第一步的分析对模型进行处理,处理后的模型中有7个构件(包括机架part2),6个铰链,2个直线导轨,一个高副,同时有3个公共约束。按文献[3]中自由度公式

W=(6-m)n-(5-m)P5-(4-m)P4-...

或 $W = (6 - m) n - \sum_{k = m + 1}^{5} (k - m) Ρ_{k} (2)$

计算系统自由度如下:

系统自由度为:6×(6-3)-6×(5-3)-2×(5-3)-1×(4-3)=1 (3)

由此可见Adams2005给出的结果与理论计算结果完全一致,说明了本文所述方法的正确性。另外,通过该实例也说明了消除模型中多余约束的必要性。

4 结论

本文通过分析约束对求解欧拉-拉格朗日方程的影响机理,给出了使用基于欧拉-拉格朗日方程的动力学分析软件时模型的约束方法及出现多余约束时的分析处理方法,并以Adams2005为例结合工程实例说明了本方法的可行性和正确性。本文给出的多余约束问题的理论分析及处理方法对使用基于欧拉-拉格朗日方程开发的动力学分析软件进行分析仿真时如何消除多余约束有一定的参考和借鉴意义。

摘要：多余约束问题是多刚体动力学分析仿真时很棘手的常见问题。针对使用基于欧拉-拉格朗日方程开发的动力学分析软件(如:Adams、DADS等)进行分析仿真时遇到的多余约束问题,进行了理论分析,即:分析约束对求解欧拉-拉格朗日方程的影响(约束向量的雅克比矩阵对求解欧拉-拉格朗日方程的影响),并给出了这一类问题的解决方法。合理地约束模型及正确地变换约束类型和方式。最后,以Adams2005为例,结合工程实例验证了该解决方法的正确性。

关键词：欧拉-拉格朗日方程,多余约束,动力学分析

参考文献

[1]张越今,宋健,张云清,等.多体系统动力学分析的两大软件.汽车技术,1997;3(3):16—20

[2]杨东武,段宝岩,狄杰建.多体系统Euler_Lagrange方程组的一类新数值分析方法.应用力学学报,2006;23(1):26—30

[3]耿德旭,欧阳富,刘彦华.空间机构自由度计算公式的重新建立.机械设计与研究,2003;19(2):21—28

[4]潘振宽,赵维加,洪嘉振,等.多体系统动力学微分代数方程组数值方法.力学进展,1996;26(1):29—39

[5]邵先喜.多体力学系统多余约束条件的数值处理方法.青岛大学学报,1997;12(4):87—89

[6]胡其彪.空间可伸展机构的设计与动力学分析研究.杭州:浙江大学,2001

仿真约束篇4

1 C型框架仿真前处理

设计的铆接机C型铆接框架结构, 材料为灰口铸铁, 牌号为HT300, σb=250MPa, 弹性模量E=1.2e5MPa, 泊松比μ=0.25。材料特性是能承受高弯曲应力, 抗拉力。

选用单元类型为SOLID45, 最小精度为0.9mm, 采用自由网格划分的法则对设计的C型铆接框架进行网格划分, 网格划分结果, 如图1所示。

对设计的C型铆接框架的固定端采用刚性全约束, 框架的下端施加均布载荷, 根据设计最大载荷与受力面积可以计算框架下端承受约4.6MPa的压力, 刚性全约束的C型铆接框架, 如图2所示。为了实现铆接框架固定端的弹性约束, 需要将液压缸等效为弹性元件 (等效弹簧) , 弹性约束的C型铆接框架, 如图3所示。据虎克定律:

式中N为液压缸受到的拉力;

l为液压缸的长度;

E为弹性模量 (206GPa) ;

A为液压缸的截面面积。

根据等效刚度的转换公式:

计算液压缸的等效刚度系数ke=2.28×106 N/m, 如图2、图3所示。

2 C型框架有限元仿真结果分析

采用ANSYS有限元CAE软件对刚性约束和弹性约束分别进行求解, 主要考核指标是C型铆接框架在X/Y/Z坐标方向上的变形情况。刚性约束的仿真结果, 如图4所示。结果表明, Y向位移最大, 变形量约为0.0 524m m, X/Z方向基本没有变形。最大压应力约为186.533MPa, 小于HT300的强度极限σb=250MPa。

弹性约束的变形结果, 如图5所示。弹性约束的固定端也有少量下移, 这种下移将会与Y向的位移叠加, 叠加后的最大变形量∆S=0.0647-0.0033=0.0614mm, X/Z向无明显变形。结论表明, C型框架应用弹性约束时变形比刚性约束的变形略大, 但这个变形量并不影响铆接机的铆接精度。基于弹性约束的C型铆接框架的最大压应力约为197.53MPa, 强度满足要求。

分析C型铆接框架刚性/弹性约束仿真结果可知, 刚性约束是的变形量和最大压应力都小于弹性约束, 而弹性约束的结果与实际实验结果更加接近。仿真结果表明, 传统的刚性约束仿真结果误差较大。

3 结语

(1) 采用CAE仿真方法对C型铆接框架进行了刚度和强度校核, 校核结果与实验结果基本一致, 解决了变截面C型铆接框架的校核难题。

(2) 采用弹性支承对C型铆接框架的固定端进行约束, 并将仿真结果与传统的刚性全约束对比。结果表明, 弹性支承结果更接近实际情况。

(3) 弹性支承的作用点设置在网格节点的方法切实, 可行。

(4) 采用自由划分的SOLID四边型网格进行仿真分析, 精度和运算速度都符合工程计算要求。

摘要：为了解决C型铆接框架刚度和强度的校核难题, 提出了采用ANSYS有限元CAE软件对其进行强度和刚度的仿真校核。为了提高仿真结果的精度, 本文提出了采用弹性支承代替传统的刚性支承约束C型框架固定端的方法, 采用自由划分的SOLID四边形网格。仿真结果表明, 采用弹性支承约束C型框架固定端的方法可行、精度高。

关键词：C型铆接框架,弹性支承,校核,刚度,强度

参考文献

[1]XU Gang, LU Jie, HUANG Caiyuan.Actuality&Progress in metal sheetstamping forming technique&equipment[J].China Metal Forming Equipment&Manufacturing Technology, 2004, (4) :16～22. (In Chinese)

[2]XU Weili, LIN Zhongqin, LIU Gang, et al.State and trends of auto bodypanel stamping simulation[J].ChineseJournal of Mechanical Engineering, 2000 (7) . (In Chinese) .

[3]Taylor L, Cao J, Karafillis A P.Nu-merical simulations of sheet metalforming[J].Journal of Materials Pro-cessing Technology, 1995 (50) :168～179.

[4]XU Bingkun, SHI Fazhong, CHENZhongkui.Several key problems in sheetmetal forming numerial simulation[J].Journal of Plasticity Engineering, 2001 (7) :32～35. (In Chinese) .

仿真约束篇5

在飞机应急着陆情况下, 乘员的安全分析、研究是适航安全验证中不可或缺的环节之一。传统的实验方法费时费力、操作复杂, 且成本巨大, 因此近年来, 利用数值仿真方法进行含约束系统成员的安全研究得到了日益广泛的应用。在汽车行业, 有限元仿真技术的应用已经相对比较成熟[1]。然而在民用航空领域, 应用却相对较少, 多数仍依赖实验研究[2]。

因此, 本文以某轻型飞机为研究对象, 建立了飞机座舱的有限元模型, 并在此基础上添加了包括假人和安全带的乘员约束系统有限元模型。按照CCAR-23中的规定的试验条件[3], 对飞机的应急着陆过程进行了数值模拟分析, 对飞机的应急着陆安全性能进行了评估。

2 有限元模型的建立

如图1所示, 为包含假人及其约束系统的飞机座舱有限元模型。假人模型基于实验所用的HybridⅢ, 第50百分位男性假人。假人各部分的关节用球形铰链和旋转铰链单元定义, 并利用非线性的扭矩弹簧和阻尼单元连接来模拟人体各部分之间的弹性[4]。约束系统模型主要包括安全带单元、卷收器单元, 以及触发卷收器工作的加速度传感器单元[5]。图2给出了座椅与三点式安全带系统有限元模型, 以及假人模型的初始姿态。

3 仿真结果与分析

给整个飞机结构以12.8m/s的水平初速度, 同时施加峰值为26g的三角形减加速度脉冲来模拟水平冲击, 并考虑飞机向左偏航10度。如图3所示为水平冲击过程中, 不同时刻假人的运动姿态。可以看出, 在整个冲击过程中, 安全带始终保持在假人的身体上, 并且座椅和约束系统可以很好的约束住假人。人体头部伤害判据值HIC (Head Injury Criterion) 是冲击过程中假人安全的重要指标之一, 反应了人体头部的损伤积累情况, 其计算公式如下:

式中, R (x) 为头部质心处的合成加速度, 用重力加速度g的倍数表示, T0, TE分别为冲击起始、终止时间, t1, t2为使HIC达到最大值的时间段的起始和终止时间, t2-t1=36ms。

图4给出了经过滤波处理的假人头部重心的合成加速度时程曲线。通过计算得到HIC=895, 小于1000, 满足相关规定要求。

4 结论

本文针对含乘员及其约束系统的某轻型飞机应急着陆过程, 利用大型动力有限元软件LS-DYNA进行了数值仿真分析。重点考察了假人的动态响应, 及其头部伤害判据值 (HIC) 等各项安全指标, 满足相关规范标准的规定, 说明该型飞机能够较好的满足应急着陆动态要求。本文的研究工作为有限元仿真技术在民用航空领域的应用积累了丰富的经验, 同时, 分析结果也为飞机应急着陆被动安全设计提供了一定的参考依据。

参考文献

[1]Stephen Kang, Paul Xiao.Comparison of HybridⅢRigid Body Dummy Models[R].10th International LS-DYNA Users Conference, 2008.

[2]Ala Tabiei, Charles Lawrence, Edwin L.Fasanella.Validation of Finite Element Crash Test Dummy Models for the Prediction of Orion Crew Member Injuries during a Simulated Vehicle Landing[R].10th International LS-DYNA Users Conference, 2008

[3]CCAR 23 R3, 正常类、实用类、特技类和通勤类飞机适航规定[S].2004.

[4]沈永峰, 邸建卫.含乘员约束系统的轿车正面碰撞计算机模拟计算[J].现代制造工程, 2008, 10:66-70.

仿真约束篇6

包边是连接车身封闭件(如四门两盖的内外板)最常用的一种成形工艺。包边作为车身封闭件制造的最后一道工序,其成形质量(尺寸准确性和表面光整度)直接影响车身的装配精度和外观质量。

近年来,多品种小批量生产正在成为汽车的主要生产方式,传统专用模具包边的生产方式不再适用。采用机器人操作辊轮沿零件轮廓进行辊压包边,不需要专用的模具,对不同的车型具有很强的适应性。近几年来,机器人辊压包边技术引起了汽车制造业的广泛关注,世界各大汽车公司都在积极开发和使用机器人辊压包边设备代替传统模压设备[1]。

机器人辊压包边属于约束操作,其环境力约束是车门板料的变形抗力,而在机器人末端的辊轮开始进给,与板料接触的初期,材料尚未完全屈服,此时变形抗力最大。这一环境力约束有可能超过包边机器人的关节驱动能力而导致作业中断。研究包边机器人运动、进给方式与板料变形抗力的大小,对机器人的轨迹规划、顺应等后续研究至关重要。

人们对传统模压包边过程中,工具的受力及设计优化进行了研究[2,3,4,5,6,7]:Bay等[2]的研究表明,在光滑硬质成形工具作用下,材料塑性变形中的摩擦力主要与正压力和滑移变形长度有关;Ford汽车公司研究了通过加载外部压力来改善卷边工艺过程中工件和工具的受力[3];Chu等[4]分析了包边过程中所需要的成形压力,其中考虑了摩擦力和金属板厚度的影响;Attanasio等[5]通过使用特殊的非对称渐进成形设备,优化了成形工具的受力和轨迹。王珏等[8]总结了国内外车门包边技术,揭示了包边设备设计应考虑的因素。然而,对机器人辊压包边这一新方式的研究尚不多。

本文研究机器人辊压包边的初始压入运动方式对包边辊轮所受最大材料变形抗力的影响。

1 机器人辊压包边工艺过程

车身门盖外板冲压成形后,其周边留有一圈高6~12mm、与外板本体成90°的翻折边,这条边就是待包边的部位。包边成形一般分为两步: 首先将待包边从90°翻折至45°, 称之为预包边; 然后再将其从45°翻折至0°并压实,称为终包边,如图1所示。

(a)翻边 (b)预包边 (c)终包边

机器人辊压包边系统如图2所示[1]。机器人2带动辊压机构3,沿工作台及夹具5的工件周边进行辊压,实现包边过程。通过规划机器人的运动轨迹和工作台及夹具5的配合,可以对不同形状的零件进行包边。

1.机座 2.机器人 3.辊压机构 4.工件 5.工作台及夹具

机器人辊压包边时,辊轮先在合适的位置以某种压入方式将翻折边从90°压至45°或从45°压至0°,然后在工作台及夹钳的配合下,机器人带动辊轮沿工件周边进行辊压完成预包边或终包边,同时,辊轮绕自身轴线转动以减小摩擦阻力。

对于通用压力机模压包边,压力机的行程一般为竖直向下,这就决定了它只能以固定的竖直压入方式进行包边。而机器人辊压包边则与之不同,机器人具有六个空间自由度,包边时可以灵活调整辊轮的空间位姿,于是可以采用多种压入方式。可能的压入方式有竖直压入、水平压入和环向压入三种,如图3所示。竖直压入方式对应于模压包边的一个截面状态,在压入过程中,辊轮的运动方向始终沿竖直方向;水平压入方式与竖直压入方式类似,辊轮的运动方向为水平方向;环向压入时,辊轮向外板弯曲方向转动,最终转角与预包边角或包边角相等。

(a) 竖直压入 (b)水平压入 (c)环向压入

包边过程中,辊轮受到来自工件的反力,该反力通过操作空间与关节空间的映射关系,与关节驱动力矩相对应;同时,不同的压入方式下辊轮受到的最大反力也不同。因此,为保证机器人各驱动关节正常工作,准确获得不同压入方式下包边辊轮所受的最大反力十分关键,这对包边机器人的选用和包边工艺的选择是非常重要的。

2 板料变形过程的有限元仿真

2.1 压入过程有限元模型的建立

对辊压包边的开始压入过程(包含竖直、水平、环向三种压入方式)进行了有限元模拟。外板采用S4R单元,材料为铝合金AA6111-T4,假设材料为各向同性,考虑材料的摩擦力。包边辊轮视为刚体。模型具体几何参数为:外板长L1=100mm,宽W=50mm,翻折边高度H=10mm,翻边模半径Rd=0.5mm,板厚t=1.25mm;辊轮直径D=16mm,长L2=20mm。三种压入方式下辊轮与铝板装配关系的有限元模型如图4所示。

(c)环向压入

2.2 板料变形仿真结果与分析

竖直压入方式的仿真结果如图5所示。

(b)45°压至0°

由图5a知,由90°压至45°时,辊轮受到沿x方向和y方向的反力Fx、Fy较大,且Fx、Fy均随着辊轮竖直压入位移d的增大而增大,变化趋势相同;沿z方向的反力Fz很小。实际上,随着辊轮向下位移的增大,材料由弹性变形过渡到塑性变形,由于材料强化的影响,在弯曲过程中需要的弯矩M不断增大,于是辊轮继续下压所受的反力不断增大。辊轮中心轴线与x方向和y方向均成45°,所以理论上沿这两个方向上的分力相等。在z方向,辊轮沿板长方向所受的来自工件两侧的反力基本相等,因而合力为零。在这个过程中辊轮所受的最大反力(合力)为1862N。

而由图5b知,由45°压至0°时,辊轮受到沿y方向的反力较大,远远大于x方向和z方向的反力分量,x方向的反力分量基本为零。这一过程中辊轮轴线与x轴夹角为零,因而辊轮受到的反力主要是沿y方向的反力。另一方面,反力随着竖直位移的增大而增大,在压入过程临近结束时反力迅速增大。根据薄板弯曲理论,弯曲半径越小,弯曲越困难。当弯曲半径达到一定值时,进一步弯曲所需要的弯矩比前期弯曲相同角度所需要的弯矩大得多,并且呈非线性增长。此过程辊轮所受的最大反力为2211N。

环向压入方式的仿真结果如图6所示。

(b)45°压至0°

由图6a知,环向压入时,辊轮在初始极小角度范围内与工件不接触,反力为零。当达到一定的环向位移量后,辊轮开始受力。由90°压至45°的压入过程中,辊轮所受的反力主要沿x方向和y方向,沿z方向的反力分量较小。此过程最大反力为1214N。

而由图6b知,45°压至0°的压入过程中,辊轮反力主要沿y方向,x、z方向的反力分量很小。同样,根据薄板弯曲理论,弯曲半径越小,所需的弯矩越大,所以最大反力出现在包边过程结束时,约为1394N。

在三种压入方式下,工件的最终变形见图7,辊轮受到的最大反力(合力)见表1。由图7可知,三种压入方式都能使工件产生相同的变形效果,达到预包边的目的。

(c)环向压入

3 铝合金板料压弯试验

为了验证仿真结果,在万能材料试验机上进行了与竖直压入方式类似的压弯试验,结果如图8所示。

(b) 45°压至0°

由图8a知,从90°压至45°时,辊轮受到的最大反力为1588N。由于试验增加了内板,当外板翻折边被压至与内板贴合时,材料的受力状态已由弯曲转变为单向压缩,继续下压会使辊轮反力显著增大。实际上,当位移量大约为3mm时外板已弯至0°,此时辊轮的最大反力约为2100N,见图8b。

为了深入了解包边机实际工作情况,我们多次到上海大众某厂参观,并与现场工程师进行了交流。该厂的四门两盖车间均采用了机器人辊压包边技术,整个工艺由多个机械手流水操作。在辊压包边一环,机器人操作辊轮以环向压入的方式开始整个包边过程。

4 结论

(1)包边机器人采用不同的压入运动方式,末端辊轮所受的最大反力不同,竖直压入时反力最大,环向压入时反力最小。在某种固定的压入方式下,辊轮所受的反力随着弯曲半径的减小而增大,最大反力出现在压至0°时。

(2)由于环向压入时辊轮受到的反力最小,机器人辊压包边采用这种压入动作方式是比较合理的。这一点在某汽车厂调研中得到了证实。

参考文献

[1]TOYOTA.Roll Hemming Method for Joining Panels,Involves Performing Pri mary Bending of Outer PanelFlange,Using Pri mary Roller Having Outer DiameterLarger than Outer Diameter of Secondary Roller:JP,2002263756-A[P].2002-01-09.

[2]Bay N,Wanhei m T.Real Area of Contact and FrictionStress at High Pressure Sliding Contact[J].Wear,1976,38:201-209.

[3]Golovashchenko S F,Bessonov N M.Development ofSharp Flanging Technology for AluminumPanels[C]//NUMISHEET 2005:Proceedings of the 6th Interna-tional Conference and Workshop on Numerical Si mula-tion of 3D Sheet Metal Forming Process.Detroit:American Institute of Physics.[2009-02-15].http://adsabs.harvard.edu/abs/2005AIPC..778..687G.

[4]Chu Seok-Jae,Lee Won-Kyu.Determination of theRequired Punch Force in Hemming[C]//ProceedingsKORUS 2003.7th Korea-Russia International Sympo-siumon Science and Technology.Soeul,2003:290-295.

[5]Attanasio A,Ceretti E,Giardini C.Opti mization ofTool Path in Two Points Incremental Forming[J].Journal of Materials Processing Technology,2006,177(1/3):409-412.

[6]Zhang Guohua.Analysis and Opti mization of SheetMetal Flanging and Hemming Process[D].Ann Arbor,Michigan:University of Michigan,2001.

[7]Livatyali H,Laxhuber T,Altan T.Experi mental In-vestigation of Forming Defects in Flat Surface-convexEdge Hemming[J].Journal of Materials ProcessingTechnology,2004,146:20-27.

【仿真约束】推荐阅读：

法律约束07-14

安全约束07-21

目标约束05-12

运行约束06-10

有效约束06-25

条件约束06-28

约束理论07-07

约束资源07-08

参数约束07-08

金融约束08-02

>> 查看更多相关文档