分布载荷

2024-04-30

分布载荷(精选三篇)

分布载荷 篇1

滚珠丝杠副作为一种新型的螺旋传动元件, 由于其高精度、高效率、优越的耐磨损性和高速特性以及运动可逆性等良好的机械性能, 它在国内外工程技术界已经获得了广泛的应用, 尤其是近年来其在航空工业的各个部门更是获得了高度的重视。而研究滚珠丝杠副的载荷分布情况是研究其一系列机械性能的基础。但是, 目前在研究滚珠丝杠副载荷分布时都是简单地按照滚珠均匀受载来分析, 这与实际的载荷分布情况有较大的偏差。因此, 建立一个求解滚珠丝杠副载荷分布的模型显得十分必要。

1 数学模型的建立

1. 1 滚珠丝杠副滚珠等效受力模型的研究

首先对单螺母滚珠丝杠副在丝杠受拉和螺母受压状态下 ( 以下简称T - C状态) 的受力进行分析, 如图1 所示。图中Fa0为滚珠丝杠副的轴向外载荷; Ds为丝杠的外径; Dn为螺母的外径; Fsi, Fni为丝杠和螺母对第i个球的轴向力。

为了有效地分析滚珠与滚道接触时的受力情况, 将某个滚珠和滚道接触时的状态放大, 如图2 所示。图中 βj为第j个滚珠和滚道间的接触角 ( j = 1, 2, 3, …) ; Pj为第j个滚珠所受的法向载荷db为滚珠的直径; t为曲率比。

为了方便分析, 结合图1 和图2 建立一个如图3 所示滚珠丝杠副滚珠等效受力的模型图, 它清楚地显示了滚珠和滚道的接触状态。

根据平衡受力情况和赫兹理论, 丝杠和螺母对第i个球的轴向力, Fsi和Fni, 可以由式 ( 1) 得来:

式中: Fa0———所受的轴向载荷;

Pj———第j个滚珠所受的接触力;

λ———滚珠丝杠的螺旋角。

由滚珠轴向外载荷和滚珠间的作用力关系可得:

式中: Z———滚珠的个数;

Pjsinβjcosλ———第j个球和滚道间接触力的轴向投影。

式 ( 2) 就是得到的轴向外载荷与每个滚珠法向载荷之间的关系。

1. 2 考虑误差情况下变形协调关系的建立

在分析了滚珠丝杠副滚珠受力状态之后, 来考虑滚珠受力时的变形情况。如图4 所示, 为了简化分析, 把第i个滚珠, i - 1 个滚珠, 螺母和丝杠轴的延螺旋角方向的接触状态放大出来。

接下来研究滚珠受力时的变形协调关系。首先, 给丝杠提供一个轴向力, 这个轴向力使所有的滚珠都保持T-C状态。实际受载后, 滚珠及滚道的几何误差 (主要是滚珠的直径误差、丝杠副的导程误差) 会导致的滚珠接触位置有所偏移。假设由于几何误差导致的所有接触位置误差不会改变初始的接触角。如图4所示, 空间变形量在螺旋角的平面内。在接触点a和b, 几何误差分别为σisinβi, σi-1sinβi-1。

图5 显示的是丝杠单元的受力变形分析。ΔL表示第i个滚珠和i-1 个滚珠间的距离。第i个球和滚道的接触点a和c的变形量 εa和 εc, 它们会随着滚珠丝杠的轴向变形而改变。同样的, 还有第i-1 个球和滚道的接触点b和d的变形量 εb和 εd。

在T-C状态下, 丝杠c点和d点之间的轴向位移变化量是 εsi, 螺母a点和b点之间的轴向位移变化量是 εni。他们可以由式 ( 3) 和式 ( 4) 得到:

结合上述两式, 可以得到滚珠的变形协调关系, 即:

这里 ( εc- εa) 和 ( εd- εb) 分别是第i个球和i-1 个球的轴向变形。当第i个滚珠及i - 1 个滚珠和滚道接触时, 相应轴向变形将抵消丝杠的几何误差 σisinβicosλ, σi - 1sinβi - 1cosλ 和滚珠接触变形的轴向投影 δisinβicosλ 和δi - 1sinβi - 1cosλ, 这就是所要研究的位移协调理论。将他们代入式 ( 5) , 就可以得到考虑误差情况下的变形协调条件:

式 ( 6) 就是滚珠在考虑误差情况下变形协调关系。

1. 3 滚珠受载模型的建立

本小节根据误差情况下的变形协调条件, 建立滚珠受载分布的模型, 假设滚珠丝杠副每个滚珠的接触角在每一时刻都相等。根据式 ( 6) 可得:

式 ( 7) 中, δi、δi - 1是第i个滚珠及i - 1 个滚珠和滚道接触时, 相应接触变形变形量, ( σi - 1- σi) 为误差项。对于间距 ΔL内的丝杠轴和螺母在T - C状态下得变形量 εni和εsi可以由式 ( 8) 和式 ( 9) 得到:

这里En, Es分别为丝杠和螺母的弹性模量。An, As分别为丝杠和螺母的有效截面区。假设丝杠的导程为L, 每个导程内的滚珠数为M, 则两个滚珠间的轴向平均距离ΔL为:

由于主要的接触曲率随着方向变化, 而且接触区是一个椭圆, 所以形变量 δ 可以由赫兹弹性接触理论得到:

这里K为常数, 它是由接触点的主曲率半径, 材料的弹性系数等决定。

将式 ( 1) 、式 ( 8) 、式 ( 9) 和式 ( 11) 代入式 ( 6) 中得:

假设每个滚珠的接触角都相等, 则式 ( 12) 化为:

式 ( 13) 反映了第i个滚珠及第i - 1 个滚珠法向接触载荷Pi与Pi - 1的关系, 其中 ( σi - 1- σi) 为误差项, 接触角 β 与螺旋升角 λ 都是已知的, 而其他参数K, εsi, εni都是可以通过上述相应的计算式得出的。

因此, 式 ( 2) 与式 ( 13) 组成了滚珠丝杠副变形协调条件及几何误差下的载荷分布模型, 即:

这个滚珠受载模型可以运用牛顿迭代法计算得到, 通过计算出每个滚珠的受载情况, 就可以分析出整个滚珠丝杠滚珠的受载分布情况。

此外, 上述模型是滚珠丝杠副在丝杠受拉和螺母受压状态时 ( 即T - C状态) 的建立分析的。当然, 其他外部载荷条件 ( 如假设单螺母和丝杠都承受张力) 的非线性方程也可以根据上述的方法构造建立。

2 滚珠受载实例计算与误差分析

2. 1 无误差时滚珠受载分析

现就滚珠丝杠副滚珠载荷分布问题, 进行了载荷分布的理论计算与推导, 举例对某型号滚珠丝杠副的滚珠受载进行计算。各项参数如下: 丝杠公称直径d0=16 mm, 螺母外径Ds= 30 mm, 滚珠直径db= 3. 175 mm, 导程L = 4 mm, 滚道曲率半径和滚珠半径之比t = 1. 04, 滚珠圈数n = 4, 每个滚珠的接触角 β0= 45° , 螺旋升角 λ= 4. 55° , 工作滚珠总数Z = 52, 轴承钢的弹性模量E = 2. 07 × 105MPa, 泊松比 μ = 0. 29, 且给螺母施加30k N的轴向外载荷。

由式 ( 10) 得两个滚珠间的轴向平均距离 ΔL为:

丝杠和螺母有效截面积为:

计算K, 由查文献[1]可知:

计算滚珠与滚道面出的曲率和∑ ρ

参数ma, mb, mamb以及可以根据曲率函数查表得到。代入可得K:

由Kn的值得则式 ( 12) 中Kp的值:

当滚珠直径无误差时, 式 ( 15) 迭代方程中 σi - 1- σi= 0, 式 ( 14) 模型可变为:

运用牛顿迭代法计算出每个滚珠的受载情况, 求得最大法向接触载荷P1= 0. 928 k N, 最小法向接触载荷P52=0. 711 k N。以滚珠数为横坐标, 第i个滚珠法向载荷Pi与最小载荷P52之比为纵坐标, 得到滚珠丝杠副的载荷分布情况。

为了验算本模型计算结果的准确性, 以滚珠数i为横坐标, 第i个滚珠接触载荷与最小接触载荷之比Pi/ Pmin, 将本模型在无误差情况的计算结果和半螺距模型[2]进行比较, 初始的接触角 β0= 45°。

在T - C状态的情况下, 半导程模型和本文模型的比较结果如图6 所示。由图可知, 第i个滚珠的接触载荷随着i的增大而呈减小的趋势。而且, 本模型滚珠的负载分布曲线介于半螺距模型曲线之中, 这表明本模型计算结果与半导程模型结果是吻合的。

2. 2 滚珠几何误差对载荷分布的影响

2. 2. 1 滚珠几何误差的计算分析

图4 所示的, 滚珠几何误差可以转化成本模型所计算的误差。在本文建立的模型中, 滚珠几何误差对载荷分布的影响是可以分析计算的。

丝杠滚珠的直径误差实际上是一种随机误差, 很难分析, 但是在一定的统计范围内误差对滚珠受载的影响是可以分析的。假设给定一组滚珠的直径误差, 如图7 所示, 该组平均几何误差为- 2 μm。对于滚珠直径误差的分析, Masaomi Tsutsumi等[6]结合数值实验的结果分析指出, 对于一组随机误差的滚珠直径, 当误差的平方差Sd≤5 μm时, 滚珠所受的法向载荷变化量 ΔP≤1 N ( 不到1%的变化量) 。因此当给定一组误差的平方差Sd≤5 μm时, , 由于误差随机性引起的载荷变化量可以忽略, 可以直接用平均几何误差来代替计算。

当滚珠平均几何误差为- 2 μm时, 式 ( 15) 迭代方程可以转化为:

这样联立式 ( 17) 迭代方程和式 ( 2) 即可以求出所有滚珠的接触载荷。

同理可以分析求出滚珠平均几何误差分别为+ 1μm, + 2 μm, - 4 μm时的滚珠受载情况。

2. 2. 2 误差对滚珠受载的影响分析

现研究误差对滚珠载荷分布的影响。首先, 假设位于丝杠中间的第i个滚珠的直径误差是负误差, 对于不同的负误差, 计算出每个滚珠的接触载荷。我们来研究当误差σi - 1- σi= - 2 um和 σi - 1- σi= - 4 um两种负误差情况下的滚珠接触载荷分布情况。首先计算出每个载荷Pi与最小载荷P52之比, 然后再画出载荷分布图, 载荷分布图如图7 所示。由图可知, 第i个滚珠的接触载荷呈减小的趋势, 减少的量与误差的比例有关。从图中还可以得出结论, 滚珠丝杠副工作滚珠中的载荷会随着负误差的增加而增加。

接下来, 再研究误差 σi - 1- σi= + 1 um和 σi - 1- σi= + 2 um两种正误差情况下的滚珠接触载荷分布情况。同样地计算出每个载荷Pi与最小载荷P52之比, 画出载荷分布图, 载荷分布图如图8 所示。从图中可以得出结论, 滚珠丝杠副工作滚珠的载荷随着正误差的增加而减少。

3 结论

根据以上丝杠几何误差对滚珠丝杠载荷分布影响的分析, 可以得出以下结论:

1) 本文提出的滚珠丝杠机构载荷分布计算模型可以用来计算滚珠几何误差在一定分布规律内的滚珠载荷。

2) 接触位置的负误差会导致滚珠在该位置的赫兹接触载荷的减小, 正误差会导致载荷的增加, 增加量和减小量与丝杠其他滚珠共同分担。

3) 上述中一些运用于单螺母结构的结论也适用于双螺母结构。

摘要:滚珠丝杠副受外载荷作用时滚珠的受力是比较复杂的, 但目前在研究滚珠丝杠副滚珠受力时都是简单地按照均载来分析, 这与实际的载荷分布情况有较大的偏差。建立一个变形协调条件下的滚珠丝杠副载荷分布的模型, 并分析了几何误差对滚珠受载的影响。

关键词:滚珠丝杠副,载荷分布,几何误差

参考文献

[1]T.A.Harris, 滚动轴承分析[M] (1、2卷) , 北京:机械工业出版社, 1967.

[2]M.Izawa, H.Shimoda, Study on the load distribution in the ball screw, Japan Journal of Precision Machine 42 (11) (1976) 1021-1028.

[3]程光仁, 施祖康, 张超鹏.滚珠螺旋传动基础[M].北京:机械工业出版社, 1987.

[4]万长森.滚动轴承的分析方法[M].北京:机械工业出版社, 1987.

[5]饶振纲, 王勇卫.滚珠丝杆副及自锁装置[M].北京:国防工业出版社, 1990.

基于载荷二维分布的可靠性分析方法 篇2

机械零件及系统的载荷一般可以用复杂的时间-载荷历程来描述, 也即可以用随机过程来对载荷历程进行数学描述。因此, 当零件承受载荷多次作用时, 则不能简单地应用载荷-强度干涉模型进行计算, 而应该从载荷历程的统计特性出发, 考虑载荷作用次数对于零件及系统可靠度的影响。近年来, 有学者对载荷多次作用下零件的可靠度计算进行了研究[1,2,3,4,5]。文献[4]假设强度不退化, 分析了载荷作用次数对可靠度的影响。文献[5]通过递推公式得到多级载荷下考虑强度退化的可靠度计算。但是, 这些模型中没有考虑载荷的随机性, 并且没有对载荷的随机历程进行详细阐述, 导致模型中随机变量物理意义不明确, 无法全面表达随机载荷多次作用时的可靠度变化规律。

本文首先对载荷历程进行了分析, 提出随机载荷的二维分布:横向分布和纵向分布;然后, 根据组成系统的零件间的结构关系, 提出考虑共因失效时机械系统的可靠度计算公式。

1 基于载荷纵向分布的零件可靠性分析

载荷-强度干涉模型是机械零件静态可靠度分析的基础。用fs (s) 、fr (r) 分别表示载荷和强度的概率密度函数, Fr (·) 表示强度的分布函数, 则载荷小于强度的概率也即可靠度为

这就是Freudenthal于1947年提出的著名的载荷-强度干涉模型[6]。其中, 载荷和强度是广义的, 载荷可以是应力、温度、腐蚀、载荷的作用次数等, 强度可以是疲劳强度、抗热性、抗腐蚀性、零件的失效循环数等。该模型适用于随机载荷作用一次或者考虑个体差异的恒幅载荷作用下的可靠度计算。实际上, 零件往往承受随机载荷的多次作用, 这时, 载荷需要用一个随机过程来描述。工程上, 经常通过多次记录载荷的整个时间历程来获得多个样本, 并通过对样本的统计分析得到载荷的统计特性。

对于确定时刻 (或者确定载荷作用次数) 通过统计得到的载荷概率密度函数, 称为载荷的纵向概率密度函数, 而不同时刻的纵向概率密度函数一般是不同的。用fsi (si) 、Fri (ri) 分别表示第i次作用时载荷概率密度函数和强度的分布函数, Ri表示载荷作用第i次时零件不发生失效的概率, 假设零件的初始可靠度为R0, 则随机载荷作用n次 (n≥1) 后的可靠度为

R (n) =i=0nRi=R0i=1n[-fsi (si) (1-Fri (si) ) dsi] (2)

这就是基于载荷纵向分布的零件可靠度模型。

2 强度不退化时基于载荷横向分布的零件可靠性分析

对于载荷整个时间轴的统计分析得到的概率密度函数称为载荷的横向概率密度函数。通常使用的各种计数法 (如雨流计数法等) 得到的分布是这种分布。值得注意的是, 每个样本得到的概率密度函数一般都不相同, 通常使用的概率密度函数可以看作是载荷横向分布的一次实现 (只统计一个样本) , 或者载荷横向分布的平均分布。为了对此进行描述, 用Xj表示载荷横向分布某一具体实现的统计特征参数 (例如对于正态分布, 其统计特征参数为均值和方差;对于三参数威布尔分布, 其特征参数为位置参数、形状参数、尺度参数) , 定义F (X) 为载荷特征分布函数, fX (Xj) dX表示载荷横向分布为某一特定分布的概率, 并且该分布的统计特征参数为Xj, fX (Xj) 为载荷特征概率密度函数。

对于某一确定样本统计得到的横向分布, 假设该确定分布的概率密度函数为fsj (s) (注意:该概率密度函数不是特征概率密度函数) , 特征参数为Xj, 载荷作用次数为n, 则某一确定载荷si出现的次数为ni=nfsj (si) ds。令强度的概率密度函数为fr (r) , 则si作用ni次后零件的可靠度为

Rsi= (∫∞sifr (r) dr) ni= (1-Fr (si) ) ni

那么, 随机载荷作用n次后的可靠度为

RXj (n) =iRsi=i (1-Fr (si) ) ni=i (1-Fr (si) ) nfsj (si) ds=iexp[nfsj (si) dsln (1-Fr (si) ) ]=R0exp[nfsj (si) dsln (1-Fr (si) ) ]=R0exp[n-fsj (s) ln (1-Fr (s) ) ds] (3)

由式 (3) 可知, 载荷每次作用的平均可靠度为

RXj (1) =R0exp[∫∞-∞fsj (s) ln (1-Fr (s) ) ds] (4)

考虑到载荷横向分布的不确定性, 得到随机载荷作用n次后零件的可靠度为

R (n) =R0∫fX (Xj) exp[n∫∞-∞fsj (s) ln (1-Fr (s) ) ds]dX (5)

式 (5) 为基于载荷横向分布的零件可靠度计算公式。

为了便于工程应用, 对于已经得到的多个样本, 可以利用式 (5) 的离散化形式进行可靠性评估, 也即先通过计数拟合出每个样本的概率密度函数, 再按照下式进行计算:

R (n) =R0jΡjexp[n-fsj (s) ln (1-Fr (s) ) ds] (6)

式中, Pj为经过拟合后载荷概率密度函数为fsj (s) 的概率, 并且可以近似地认为该概率等于概率密度函数为fsj (s) 的样本个数占全体样本的比例。

这样就得到了基于载荷横向分布的、适合实际应用的零件可靠度计算公式。

3 考虑强度退化时基于载荷横向分布的零件可靠性分析

以上分析中假设了强度不退化, 当考虑强度退化时, 可以采用以下的近似计算方法:

对于某一确定的横向概率密度函数fsj (s) , 当载荷作用n次时, 大小为si的载荷的出现次数为ni=nfsj (si) ds次, 其造成的损伤为nfsj (si) ds/Nsi, Nsisi作用下的零件使用寿命。由线性损伤等效原理 (这里只是为了叙述方便, 实际上可以采用其他非线性等效原理) , si作用ni次, 相当于大小为sj0的载荷作用ni0次, 并且

nfsj (si) dsΝsi=ni0Νsj0

式中, Nsj0为sj0作用下的零件使用寿命。

又由零件的S-N曲线, 近似得到Nsm=C, 其中mC近似为常数, 得到

ni0=nsimfsj (si) dssj0m

因此随机载荷作用n次相当于载荷sj0作用nsj0次, 并且

nsj0=-nsimfsj (si) dssj0m (7)

于是, 可以认为强度符合确定载荷sj0作用下的变化规律, 而且sj0作用第q-1次时的剩余强度分布函数为Frq-1|sj0 (r) 。于是得到

R (n) =R0fX (Xj) q=1-nsimfsj (si) dssj0m (1-Frq-1|sj0 (sj0) ) dX (8)

特别地, 当sj0= (-simfsj (si) ds) 1m时, nsj0=n, 也就是说, 随机载荷作用n次时, 相当于确定载荷sj0作用n次, 也即可以认为零件的强度近似遵循确定载荷sj0作用下的变化规律, 并且该等效载荷的形式与英国PD6493 1991推荐使用的用于非恒幅载荷平面缺陷的评定的、基于Paris公式的等效载荷完全相同, 从而证明了公式的正确性。令S= (-simfsj (si) ds) 1m, 则式 (8) 变成

R (n) =R0fX (Xj) q=1n (1-Frq-1|S (S) ) dX (9)

同样, 式 (8) 的离散化形式为

R (n) =R0jΡjq=1-nsimfsj (si) dssj0m (1-Frq-1|sj0 (sj0) ) (10)

4 横向概率密度函数与纵向概率密度函数的关系

以上分析了基于载荷横向分布和纵向分布零件可靠度的计算模型。从模型的建立过程中可以发现, 基于零件纵向分布的可靠度模型需要对各个时刻 (或者不同载荷作用次数) 的载荷概率密度函数进行统计。然而这在实际上是做不到的, 因为这需要记录大量的载荷-时间历程样本, 并且进行大量的统计。为便于工程应用, 往往假设载荷历程为各态历经过程, 也即认为各个时刻的载荷概率密度函数相同。这对载荷-时间历程样本数目和时刻的选择提出了较高的要求, 样本应该具有一定的代表性。更重要的是, 当得到的载荷-时间历程样本数量不足, 不能通过各个时刻的统计得到各态历经假设下的载荷概率密度函数, 并且每个样本计数后拟合的概率密度函数相差比较大时, 基于载荷横向分布的可靠度模型则可以充分地利用有限的样本直接对零件可靠度进行近似计算。

下面通过分析载荷横向分布和纵向分布概率密度函数的关系, 给出在有限样本条件下, 根据载荷的横向分布统计特性近似得到载荷纵向分布概率密度函数的方法, 从而可以方便地使用基于载荷纵向分布的可靠度模型。这种方法避免了对样本进行各个时刻的统计, 更适于利用传统可靠性模型进行计算。假设载荷历程为各态历经过程, 每个时刻的概率密度函数相同, 均为fs (s) 。那么每次载荷大小为s0的概率为P0=fs (s0) ds, 该载荷不出现的概率为1-P0。如果载荷作用了n次 (假设n足够大) , 那么该载荷出现的次数n0的期望值为E (n0) =nP0=nfs (s0) ds。另外, 在载荷作用的n次中, 载荷s0出现次数为x的概率, 即

P (n0=x) =CxnP0x (1-P0) n-x

P0为定值, 对P (n0=x) 取极值, 得到x=nP0=nfs (s0) ds。由此可以近似认为各个载荷-时间历程样本中, 通过计数统计得到的出现概率最大的概率密度函数就是各态历经假设下载荷的纵向概率密度函数。或者分别求出各个样本中相同级别大小的载荷的平均出现次数, 拟合出载荷的平均横向概率密度函数, 并将其作为纵向概率密度函数。这样就可以结合载荷-强度干涉模型, 使用基于载荷纵向分布的可靠度模型计算零件及系统的可靠度。

假设不考虑强度退化, 零件强度服从N (600, 302) MPa的正态分布, 载荷分布服从N (500, 302) MPa的正态分布, 初始可靠度为1。当载荷的概率密度函数为纵向概率密度函数时, 利用式 (2) 可以得到可靠度随载荷作用次数的关系。当载荷的概率密度函数为横向概率密度函数时, 应用式 (3) 得到可靠度随载荷作用次数的变化规律。二者的对比如图1所示。

1.纵向概率密度函数 2.横向概率密度函数

由图1可以看出, 两种情况下零件可靠度随载荷作用次数的增加而下降。对于相同的载荷作用次数, 当载荷的概率密度函数为横向概率密度函数时, 零件可靠度小于概率密度函数为纵向概率密度函数时的可靠度, 并且开始随着载荷作用次数的增加, 其差值逐渐增大, 当载荷作用次数较大, 零件趋于失效时, 其差值逐渐缩小。产生差值的原因是, 纵向概率密度函数是横向分布中出现概率最大的概率密度函数或者是平均横向概率密度函数, 但是出现概率并非绝对等于1。从图1中可以看出, 二者在一定载荷作用次数范围内差值并不是很大, 可以进行近似计算, 并且能够反映出零件可靠度的基本变化规律。所以, 当获得的载荷-时间历程样本比较少时, 也可以采用以上提出的近似方法, 使用基于纵向载荷分布的可靠性模型进行零件的可靠度计算。从另外一个角度讲, 在本例中, 如果假设基于载荷纵向分布的计算结果是准确值, 那么基于载荷横向分布的计算结果就是估算值。而由之前的推导过程可以看出, 与纵向分布概率密度函数相同的横向概率密度函数出现的概率最大。因此更加说明了在有限数目样本情况下使用基于载荷横向分布的可靠性模型进行近似估算的合理性。

5 考虑共因失效时基于载荷横向分布的系统可靠度分析

以上着重分析了基于载荷横向分布的零件可靠度计算模型, 对于由零件组成的系统, 可以根据系统结构功能函数求解系统的可靠度。本文着重分析考虑共因失效的系统可靠度。当系统多个零件受到相同载荷作用时, 如果作零件相互独立的假设将引起较大的误差。目前, 研究人员提出了很多考虑共因失效的可靠度计算模型, 如α因子模型、β因子模型、BFR模型等。但是这些方法仍然基于元件层面, 物理意义不明确[7,8,9]。本文直接应用之前建立的可靠度模型, 分析当系统元件承受同一载荷多次作用时的可靠度变化规律。

为了叙述方便, 以下均假设组成系统的零件相同, 并且不考虑强度退化, 每个零件具有相同的强度概率密度函数和分布函数fr (r) 、Fr (r) 。对于由a个零件组成的串联系统, 当大小为确定值s的载荷作用一次时, 系统可靠度为

Rseries1= (1-Fr (s) ) a

假设载荷的横向分布概率密度函数为fsj (s) , 特征概率密度函数为fX (Xj) , 则利用与式 (3) 类似的推导过程, 可以得到随机载荷作用n次后串联系统的可靠度:

离散化模型为

同理, 由a个零件组成的并联系统, 当大小为确定值s的载荷作用一次时, 系统可靠度为

Rparallel1=1- (Fr (s) ) a

随机载荷作用n次后并联系统的可靠度为

离散化模型为

R¯parallel=R0jΡjexp (n-fsj (s) ln[1- (Fr (s) ) a]ds) (14)

a个零件组成的k/a系统, 当大小为确定值s的载荷作用一次时, 系统可靠度为

Rk/a1=i=kaCai (1-Fr (s) ) i (Fr (s) ) n-i

随机载荷作用n次后k/a系统的可靠度为

Rk/a=R0fX (Xj) exp (n-fsj (s) ln (i=kaCai (1-Fr (s) ) i (Fr (s) ) n-i) ds) dX (15)

离散化模型为

R¯k/a=R0jΡjexp (n-fsj (s) ln (i=kaCai (1-Fr (s) ) i (Fr (s) ) n-i) ds) (16)

以上模型在建立的过程中考虑了零件承受的载荷为同一载荷这个信息, 并且没有作零件相互独立的假设, 因此反映了共因失效这种失效相关性。

以上分析了随机载荷作用下系统可靠度随载荷作用次数的变化规律。如果考虑载荷出现次数的随机性, 则可以得到当强度退化时零件的时变可靠度计算模型。如果载荷出现的次数满足以下条件[10]:①当t=0时, n (0) =0;②在互不重叠的时间段内载荷出现的次数相互独立;③在时刻t和足够小的时间段Δt>0, 有

{Ρ (n (t+Δt) -n (t) =1) =λ (t) Δt+o (Δt) Ρ (n (t+Δt) -n (t) 2) =o (Δt)

则时间t内出现n次载荷的概率为

Ρ (n (t) -n (0) =n) = (0tλ (t) dt) nn!exp (-0tλ (t) dt)

也即将载荷出现的次数看作是强度为λ (t) 的泊松过程。这时, 由全概率公式, 可以得到系统可靠度随时间的变化规律。假设系统初始可靠度为1。对于串联系统, 由式 (11) 得到的系统可靠度为

Rseries (t) =exp (-0tλ (t) dt) +n=1 (0tλ (t) dt) nn!exp (-0tλ (t) dt) (fX (Xj) exp (n-fsj (s) ln ( (1-Fr (s) ) a) ds) dX) (17)

由式 (13) , 得到并联系统的可靠度为

R (t) =exp (-0tλ (t) dt) +n=1 (0tλ (t) dt) nn!exp (-0tλ (t) dt) (fX (Xj) exp (n-fsj (s) ln (1- (Fr (s) ) a) ds) dX) (18)

由式 (15) , 得k/a系统的可靠度为

Rk/n (t) =exp (-0tλ (t) dt) +n=1 (0tλ (t) dt) nn!exp (-0tλ (t) dt) (fX (Xj) exp (n-fsj (s) ln (i=kaCai (1-Fr (s) ) i (Fr (s) ) n-i) ds) dX) (19)

假设有分别由3个相同零件组成的串联系统、并联系统及2/3系统。对已得到的载荷-时间样本进行统计, 得到的各个样本均服从双参数威布尔分布。特征参数的分布律如表1所示。λ (t) =0.6h-1, 初始可靠度为1, 零件强度服从均值为650MPa、标准差为20MPa的正态分布, 则零件可靠度随时间的变化规律如图2所示。

1.并联系统 2.2/3系统 3.串联系统

从图3可以看出, 系统可靠度随时间逐渐降低, 并且在初始阶段下降速度最快, 随后下降速度趋于平缓。其中串联系统下降速度最快, 可靠度最低, 2/3系统其次, 并联系统下降最慢。实际上该例中串联系统相当于1/3系统, 并联系统相当于3/3系统, 结果与实际相吻合。

6 结论

(1) 机械零件及系统往往承受随机载荷的多次作用, 在分析随机载荷下的零件可靠度时应该首先对随机载荷的统计特征进行分析。本文提出载荷的二维分布, 也即载荷的横向分布和纵向分布。传统载荷-强度干涉模型适用于载荷作用一次的可靠度计算。在计算随机载荷指定作用次数后的可靠度时, 应该使用载荷的纵向概率密度函数。但是, 工程上得到的是载荷-时间历程样本。由于样本数量的限制以及载荷纵向分布统计特性获得的困难, 使得基于载荷纵向分布的可靠性模型难以应用, 因此, 本文重点分析了基于载荷横向分布的可靠度计算模型。该模型可以直接根据获得的样本进行可靠度计算, 并从理论上分析了该计算方法的合理性。

(2) 本文在建立基于载荷的横向分布的计算模型时, 分别考虑了强度退化以及强度不退化两种情况, 可以应用于零件早期失效期 (强度不退化或者退化不明显) 以及偶然失效期和耗损失效期 (需要考虑强度退化) 。此外, 通过分析载荷纵向概率密度函数与横向概率密度函数的关系, 提出了根据载荷横向分布统计特征, 近似得到各态历经假设下载荷纵向概率密度函数的方法, 这样便可以结合载荷-强度干涉模型计算零件可靠度与载荷作用次数的关系。通过实例验证了该方法可以近似反映零件可靠度的变化规律, 并同时从另一个角度验证了基于载荷横向分布的可靠性模型的合理性和正确性。

(3) 由基于载荷横向分布的零件可靠度计算模型分析了几种典型机械系统的可靠度计算方法。传统的模型需要假设零件间相互独立, 但实际上由于载荷的相关性, 使得独立假设往往不能反映系统的实际可靠度。本文没有作零件相互独立的假设, 直接根据系统中零件的组成结构推导出串联系统、并联系统以及表决系统的可靠度计算方法。本文进一步考虑了载荷作用次数的随机性, 通过假设载荷作用次数为泊松过程, 建立了系统的时变可靠度计算模型。

参考文献

[1]Huang W, Askin R G.A Generalized SSI ReliabilityModel Considering Stochastic Loading and StrengthAging Degradation[J].IEEE Transaction on Relia-bility, 2004, 53 (1) :77-82.

[2]Koh C G, Ang K K, Zhang L.Effects of RepeatedLoading on Creep Deflection of Reinforced ConcreteBeams[J].Engineering Structures, 1997, 19 (1) :2-18.

[3]Thayalan P, Aly T, Patnaikuni I.Behaviour of Con-crete-filled Steel Tubes Under Static and VariableRepeated Loading[J].Journal of ConstructionalSteel Research, 2009, 65:900-908.

[4]王正, 谢里阳, 李兵.随机载荷作用下的零件动态可靠性模型[J].机械工程学报, 2007, 43 (12) :20-25.

[5]王新刚, 张义民, 王宝艳.载荷作用次数和强度退化对SSI模型的影响[J].机械设计与制造, 2008, (12) :205-207.

[6]Freudenthal A M.Safety of Structures[J].Transac-tion, ASCE, 1947, 112:125-180.

[7]Levitin G.Incorporating Common-cause Failuresinto Nonrepairable Multistate Series-parallel Sys-tem Analysis[J].IEEE Transaction on Reliability, 2001, 50 (4) :380-388.

[8]Xie L Y, Zhou J Y, Hao C Z.System-level Load-strength Interference Based Reliability Modeling ofk-out-ofnSystem[J].Reliability Engineering andSystem Safety, 2004, 84:311-317.

[9]Xie L Y, Zhou J Y, Wang Y Y, et al.Load-strengthOrder Statistics Interference Models for System Re-liability Evaluation[J].International Journal of Per-formance Engineering, 2005, 1 (1) :22-36.

分布载荷 篇3

骨骼的微损伤是指由骨骼疲劳诱发的骨显微结构的进行性变化,骨微损伤的累积可以导致应力性骨折的发生[1]。应力性骨折又称疲劳性骨折,是由于肌肉过度疲劳后无法及时吸收反复的应力过载荷刺激,而使骨骼表面发生微损伤的积累从而诱发裂纹。应力性骨折是由微损伤不断积累超过机体骨自身修复能力而产生的。应力性骨折常见于军事训练和运动员训练中,已成为部队多发病,直接影响部队官兵身体健康,削弱了部队战斗力,逐渐成为军事医学重点研究的内容之一。另一方面,骨质疏松骨折的发生也与微损伤直接相关[2]。研究证实骨骼微裂纹的自我修复与机体的骨重塑密切相关[3]:健康人的骨骼能够通过骨重塑修复微裂纹,而骨质疏松患者骨重塑能力降低,对骨骼微裂纹的修复能力衰退,造成骨骼微损伤的累积,并最终诱发骨质疏松骨折的发生。

为了更系统、深入地解析骨微损伤的发生机制,探索更有效的预防和治疗措施,首先需要精确地定量分析骨表面应变的空间分布以及其与骨微损伤发生的相关关系。1972年,Brekelmans和Rybicki等人首次将有限元法应用于矫形外科的生物力学研究领域中。在探讨股骨内部的应力分布尤其在人体的生物力学研究中,有限元法显示了其极大的优越性。20世纪80年代以后,有限元法的应用范围逐渐扩展到了股骨、颅面骨、关节、颈椎、腰椎以及其附属结构等具有骨性结构的生物力学研究中[4,5]。因此,本研究通过建立大鼠胫骨三维有限元模型,分析计算在生理水平的轴向压缩应力载荷作用下胫骨表面的应变空间分布,定量分析骨表面应变与载荷间的相关关系,并基于此确定胫骨表面最大张应变和压应变的解剖位置,为骨微损伤的研究提供重要的理论依据。

1. 方法

1.1 实验仪器和材料

雄性SD大鼠(8周龄,体质量244 g),由第四军医大学动物实验中心提供,取右侧胫骨,肌肉、骨膜等软组织充分剥离,固定于85%酒精中;显微CT(Micro CT,美国GE公司),扫描基本参数设置为:电压80 k V,电流80μA,曝光时间3 s,扫描分辨率为29μm;所用分析软件包括Mimics17软件(Materialise公司)、Geomagic Studio12软件(Geomagic公司)、Soildworks2016(Soildworks公司)、ABAQUS 6.14(ABA-QUS公司);进行有限元模型分析的软件环境为Windows 8.1 Professional 64 bit;进行有限元模型分析的计算机硬件环境包括CPU为Intel I7-4720HQ,内存为8 GB,图形加速卡为Nvidia GTX 860 M,硬盘为WDC WD10SPCX-24HWST1。

1.2 设计思路及分析流程

建立模型和有限元分析计算的流程如图1所示。首先经Micro CT扫描得到DICOM格式的大鼠胫骨断层图像,将Micro CT图像导入Mimics17中进行三维重建,保存为STL格式文件后导入Geomagic Studio12中。在Geomagic Studio中对三维重建后的模型进行去噪和曲面优化,进一步划分曲面片转化为CAD模型,存储为IGS格式[6]。运用Soildworks将IGS文件中的曲面缝合并生成实体,保存为STP格式文件。将STP导入ABAQUS 6.14中进行网格划分,存储为INP文件后导入Mimics中进行材料属性赋值[7]。有限元模型建立完成后导入ABAQUS中进行求解参数、载荷和边界条件的设置,再进行分析计算。最后,在后处理中得到应变云图并进一步采集数据进行分析。

1.3 有限元模型建立

有限元分析法作为一种求解偏微分方程的数值计算方法,具有良好的通用性和实用性,特别是对大型工程问题和一些复杂材料问题等具有很好的求解结果。但对于结构复杂和几何形状不规则的骨骼而言,直接在有限元软件中重建其三维模型难度很大。因此,需要使用到医学软件、逆向工程软件、三维CAD系统及有限元分析软件等一系列软件进行有限元模型的建立。

1.3.1 三维重建

Micro CT扫描所得DICOM格式图像输入Mimics17中进行图像的组织和编辑。扫描所得图像共1 341张,为获得更高的处理效率,将组织图像中的skip images设置为1。随后,根据大鼠胫骨的Micro CT值进行阈值选择。本实验中为使蒙罩轮廓清晰,将阈值设定为2 729~7 115,为后续蒙罩编辑做好准备工作。为填补蒙罩中的空洞,首先通过轮廓线计算,用轮廓线来填补空洞,余者使用蒙罩编辑功能对每一层面的蒙罩进行Erase或Draw编辑。蒙罩编辑完成后进行三维模型重建,重建参数采用“Optimal”(最佳的)。重建完成后进行表面光滑处理并检测是否有相交三角形的存在,经处理后将三维模型保存为由多个三角形面片定义组成的STL文件格式的模型。

1.3.2 生成实体

经Mimics三维重建得到的STL格式文件导入Geomagic Studio12中进行曲面优化。经过去噪、去除特征和网格诊断等处理进行曲面片划分,划分时注意曲面片数量不应过多。将划分好曲面片的模型转为CAD模型并存储为IGS格式,而后将其导入Soildworks 2016中进行曲面片检查、曲面缝合并生成实体,存储为STP格式。

1.3.3 网格划分

划分好曲面片的STP格式实体模型导入ABAQUS6.14中进行面网格和体网格的划分。本实验中将网格类型设置为四面体(Tet)实体单元。网格划分完成的模型节点数为56 420,单元数为36 247,而后保存模型为INP格式。

1.3.4 材料属性

将划分好网格的INP格式模型导入Mimics中进行材料属性的赋值。根据骨骼的解剖结构和生物力学特性将其化简为各向同性材料,并用3种材料来描述大鼠胫骨模型。本实验中用松质骨、密质骨及髓质3种材料来描述。采取Mimics中的均分法对胫骨进行材料属性的赋值。根据如下的经验公式(1)和(2)分别对胫骨的密质骨、松质骨和髓质进行赋值[8]。其中密度(D)单位为kg/m3,弹性模量(E)单位为MPa,泊松比均设置为0.35,A为CT值。设置完成后保存为INP文件。

1.4 有限元计算分析

有限元模型建立完成后即可导入ABAQUS 6.14中进行有限元计算分析。求解器的关键参数设置如下:设置静力通用分析步;定义时间增量增长方式,最大增量步数为100,增量步大小为0.01~0.1;应用直接稀疏矩阵求解器;仿真前利用动态循环应力加载装置对SD大鼠进行动物在体实验,进行应力性骨折动物模型的构建。由于斜坡加载模式更符合运动时胫骨的受力情况,故在动物实验和整个仿真分析步内,载荷的施加方式均采用线性斜坡加载方式[9],且采用压缩胫骨的加载方向使胫骨产生微裂纹,以此模拟生理状态下的运动状态。

所加载的位移载荷大小根据前人所做动物实验[10]及仿真前预实验来确定。轴向载荷的大小确定需保证大鼠进行在体实验时不会发生完全性骨折且实验后大鼠仍可存活。本实验中预压力5 N,加载频率为0.67 Hz,其峰值压力为35 N,分别在位移载荷为1 200、1 400、1 600、1 800和2 000μm的实验条件下进行累计3 000个周期的载荷加载。为使有限元分析结果与动物实验结果相对应,计算机仿真时将胫骨下端边界条件设置为完全固定,选取胫骨上端施加位移载荷5组,分别为1 200、1 400、1 600、1 800和2 000μm[11];并用集中力代替实际载荷情况,载荷加载的方向选择为使胫骨产生压缩形变的方向。分析完成后,在后处理中,查看应变云图及截面应变云图。通过查询信息功能来查看每一个单元及相应节点的应变等参数,统计大鼠胫骨表面张应变与压应变的最大值。

2 结果

经过基于Micro CT的三维有限元计算分析得到在位移载荷为1 400μm时胫骨表面应变空间分布的云图如图2所示。其中,正值表示张应变,负值表示压应变。从应变云图中可观察发现,在位移载荷为1 400μm时,胫骨最大张应变发生于胫骨内侧中下段,其大小为2 592με,在其两侧,向上延伸至胫骨内侧中段,张应变逐渐减小;胫骨最大压应变发生在胫骨外侧中段,其大小为5 961με,在其两侧,呈空间对称分布,压应变大小逐渐减小,且胫骨外侧压应变大小普遍大于胫骨内侧张应变大小。

经有限元分析计算所得到的最大主应变(拉伸)即表面张应变,在位移载荷为1 200、1 400、1 600、1 800和2 000μm时分别为2 209、2 592、2 978、3 338和3 704με。不同位移载荷大小与胫骨表面最大张应变大小之间关系的统计结果如图3所示。经线性回归分析可发现,位移载荷大小与胫骨表面最大张应变之间呈高度线性正相关关系,其R2为0.999 8。

经有限元分析计算所得到的最大主应变(压缩)即表面压应变,在位移载荷为1 200、1 400、1 600、1 800和2 000μm时分别为-5 067、-5 961、-6 830、-7 681和-8 563με。不同位移载荷大小与胫骨表面最大压应变大小之间关系的统计结果如图4所示。经线性回归分析发现,位移载荷大小同胫骨表面最大压应变间呈高度线性正相关关系,其R2为0.999 9。

3 讨论

本文以医学软件Mimics为载体,将Micro CT扫描获得的三维图像导入该软件中,在Mimics中利用阈值区分开皮质骨、松质骨和髓质,利于建立模型的精确性,随后利用Geomagic Studio、Soildworks等软件对其进行处理,最终建立大鼠胫骨三维模型;然后,使用有限元软件ABAQUS对该三维模型的体网格划分、单元属性设置、条件约束和实际情况模拟进行加载分析,统计分析不同载荷下的胫骨表面最大张应变和压应变的大小,并观察应变的空间分布情况。最后,将计算机仿真结果和动物实验研究的情况与参考文献的研究成果做了比对,验证了利用该方法建立模型的正确性,为微损伤发生率最高部位的预测及治疗应力性骨折等后续研究提供了有力的理论支撑。

在本研究中,我们所施加的位移载荷大小在先前的动物实验研究基础上确定,均为骨骼能够承受的生理水平的应力载荷强度[10]。而课题组通过前期的动物实验研究也证实,该应力载荷强度能够在骨骼表面产生显著的微损伤累积。应力载荷的加载方向采取胫骨轴向周期压缩载荷形式,这种应力加载形式与动物日常活动所施加在骨骼上的应力载荷形式更符合。而对于人体而言,这种应力载荷形式与人体日常的走路、跑步、深蹲等活动施加在骨骼的应力形式极为接近[12]。因此,为了更好地模拟产生骨骼微损伤的生理运动形态,故而采取使胫骨产生轴向周期压缩的载荷形式。

根据先前的临床及动物实验研究统计,胫骨65%的裂纹发生在其外侧,只有15%会发生在其内侧,而20%会发生在双踝部位[13]。在本研究中,我们观察计算机仿真结果所得的空间应变分布云图发现,最大应变值出现在胫骨外侧平台上,并且外侧平台的应变值大于内侧,这一结论与先前的临床和动物实验研究结果基本一致。此外,最大张应变位于胫骨内侧中下段,最大压应变位于胫骨外侧中段,此部位也是胫骨应力骨折的多发部位,这与临床及动物实验统计结果也相一致。而其中,胫骨平台下侧也出现了应力集中现象,这是由于下侧是约束施加的位置。这进一步说明了有限元分析的结果和实际情况是一致的,同时也证明了本研究的胫骨有限元计算模型的可靠性和准确性。故此建模方法可建立出合理可定量的三维有限元模型,并且明确了轴向压缩载荷作用下,胫骨表面张应变和压应变的空间分布,明确了最大张应变和压应变的发生位置,为微损伤发生率最高部位的预测提供了重要的理论参考。

有限元法对一些求解偏微分方程边值问题、几何形状复杂、不规则边界或厚度突变以及材料非线性、几何非线性等问题的求解表现出非常简单和精确的优点,故而在骨骼生物力学研究中,有限元法显示了其极大的优越性[14]。本实验中采用有限元分析法分析大鼠胫骨表面应变分布。先前研究应力载荷作用下大鼠股骨表面的应变分布多采用微应变片测量的方式[15]。而本研究中我们采用仿真的方法与前人使用应变片的测量结果相比较,在相同应力载荷下大鼠胫骨表面应变大小结果相匹配[15]。这进一步证明了本研究中采用有限元法定量胫骨表面应变分布的有效性和准确性。

在下一步的研究中,我们将在以下2个方面进一步完善我们的仿真计算工作。首先,在建模过程中考虑到半月板、韧带等非骨性结构的作用。其次,我们将以人体胫骨为研究对象,研究在人体运动过程中(如步态周期中),尽量使载荷和约束条件与生理环境更加接近,以此为基础进一步定量分析胫骨或膝关节的应力分布情况。

4 结语

上一篇:渠道防渗新方法下一篇:理解化学知识