复杂焊接结构随机振动疲劳寿命预测方法

摘要

本发明公开了一种复杂焊接结构随机振动疲劳寿命预测方法，包括：建立焊接系统的有限元模型；确定系统的边界条件和动力学方程，将的边界条件引入的动力学方程；在载荷输入点施加不同的外部激励载荷，进行扫频计算，获得该外部激励载荷下的节点力-位移传递函数，通过节点力计算膜应力和弯曲应力，得到通过计算焊缝处等效结构应力的传递函数；对每个载荷输入进行傅氏变换，获得每个外部激励载荷的功率谱和载荷之间的互功率谱；根据每个外部激励载荷的功率谱、互功率谱和等效结构应力传递函数得出等效结构应力功率谱；使用Dirlik法获得等效结构应力概率密度函数，统计单位应力变化范围及发生的频次；利用焊缝结构的主S-N曲线，预测焊缝振动的疲劳寿命。

说明书

技术领域

本发明涉及复杂焊接结构在随机振动条件下的疲劳寿命预测方法。涉及专 利分类号G06计算；推算；计数G06F电数字数据处理G06F17/00特别适用于 特定功能的数字计算设备或数据处理设备或数据处理方法G06F17/50计算机辅 助设计。

背景技术

焊接结构以结构设计灵活、易于结构的变更和改型、焊前准备工作简单等 独特优点，从而在汽车、铁路、航空、船舶等载运工具行业获得了相当广泛的 应用，以轨道车辆为例，没有焊接就没有其任何产品，但是它的缺点同样突出， 这就是它的焊缝上的抗疲劳能力远低于母材，疲劳失效总是从焊缝上开始。随 着载运工具的服役载荷环境越来越复杂，焊接结构上焊缝疲劳开裂的问题也日 益严重，这给载运工具的服役安全带来了巨大的疲劳隐患，因此，有人将焊接 结构比喻为“双刃剑”是有一定道理的。

焊接结构疲劳寿命预测方法主要是“名义应力法”。基于名义应力的评价方法 遇到的困难具有普遍性，因为它只是在载荷简单、接头几何形状也简单的情况 下有效，而工程上，载荷与接头几何形状相当复杂。于是，IIW提出了热点应力 法[3]，它试图通过应力外推的策略获得焊趾上的应力集中，可事实上，这种方 法因人而异的不一致性依然存在，例如，结果对有限元网格的大小非常敏感。 如果事先假设焊趾上缺口的几何形状，然后用断裂力学的方法预测其疲劳寿命， 可是缺口的假设也将因人而异，这是另外一类的不一致性问题。

在焊接结构疲劳寿命预测方面，国内使用的技术主要是从英国引进的标准 (BS7608)，以及BS7608的拓展，它们均是基于名义应力法的：首先，创建计 算对象的有限元模型；然后，计算焊接接头上的静态名义应力；接着，从BS7608 等标准中挑选出合适的S-N曲线数据[4]，并计算该应力水平下的疲劳损伤；最 后，利用线性比例关系，计算其它名义应力水平的疲劳寿命或疲劳损伤。它的 缺陷是：

(1)当内嵌的S-N曲线数据不能对号入座时，计算将难以继续进行，方法 本身有明显的局限性，而事实上，这种情况经常发生；

(2)模型中没焊缝(焊趾、焊根)的定义，因此不能计算出焊缝(焊趾、 焊根)上的结构应力，即应力集中，也不能显示焊缝上的疲劳寿命，而事实上， 这些信息最重要；

(3)由于是准静态的名义应力的计算，没有考虑疲劳载荷的频率影响，因 此，既使是有S-N数据可用，寿命评估偏差也难于判断。

(4)焊接结构服役期间的在线实测过程中很容易获得随机的位移谱、速度 谱、加速度谱，名义应力法、热点应力法对此没有对接的接口。

发明内容

本发明针对以上问题，提出一种复杂焊接结构随机振动疲劳寿命预测方法， 包括如下步骤：

—建立包含焊缝细节的焊接系统的有限元模型；

—确定系统的边界条件，建立动力学方程，将所述的边界条件引入所述的 动力学方程；

—在载荷输入点分别施加不同的外部激励载荷，进行扫频计算，获得该外 部激励载荷下的节点力-位移传递函数，通过节点力计算膜应力和弯曲应力，得 到通过计算焊缝处等效结构应力的传递函数；

—对每个实测的载荷输入进行傅氏变换，获得每个外部激励载荷的功率谱 和载荷之间的互功率谱；

—根据每个外部激励载荷的功率谱、载荷间的互功率谱和等效结构应力传 递函数得出等效结构应力功率谱；

—使用Dirlik法获得等效结构应力概率密度函数，并统计单位时间内应力变 化范围及发生的频次；

—利用焊缝结构的主S-N曲线，预测焊缝振动的疲劳寿命。

所述的边界条件中外部激励载荷至少包含：力、位移、速度和加速度。

所述确定系统的边界条件，建立多载荷输入的动力学方程的步骤中建立的 动力学方程为：

$\left[M\right]\left\{\stackrel{··}{x}\left(t\right)\right\}+\left[B\right]\left\{\stackrel{·}{x}\left(t\right)\right\}+\left[K\right]\left\{x\left(t\right)\right\}=\left\{f\left(t\right)\right\}---\left(1\right)$

其中[M]为系统的质量矩阵，[B]为系统的阻尼矩阵，[K]为系统的刚度矩阵， {x(t)}为系统的位移向量，为系统的速度矢量，为系统的加速度矢量， {f(t)}表示施加的外部激励载荷为力。

将边界条件和动力学方程同时进行傅氏变换至频域，将变换至频域的边界 条件引入动力学方程；

当施加的激励为力载荷时，经傅氏变换到频域的表达式为：

f(t)＝p(ω)·e^iωt    (2)

将公式(2)带入动力学方程公式(1)，经傅氏变换到频域后，方程为：

-ω²[M]{u(ω)}e^iωt+iω[B]{u(ω)}e^iωt+[K]{u(ω)}e^iωt＝{P(ω)}e^iωt  (3)

方程两边抵消掉复指数e^iωt，得到：

[-ω²M+iωB+K]{u(ω)}＝{p(ω)}   (4)。

当施加的激励载荷为位移时，位移经过傅氏变换到频域的表达式为：

u(t)＝u(ω)·e^iωt   (5)

所动力学方程，公式(3)中的自由度{u(ω)}，按位移激励进行重组，焊接 系统的质量阵、阻尼阵和刚度阵按自由度分块为约束自由度，下标为s和无约束 自由度，下标为f，如式(8)所示：

$({-\omega }^2\left(\begin{array}{cc}{M}_{\mathrm{ff}}& M_{\mathrm{fs}}\\ M_{\mathrm{sf}}& M_{\mathrm{ss}}\end{array}\right)+\mathrm{i\omega }\left(\begin{array}{cc}{B}_{\mathrm{ff}}& B_{\mathrm{fs}}\\ B_{\mathrm{sf}}& B_{\mathrm{ss}}\end{array}\right)+\left(\begin{array}{cc}{K}_{\mathrm{ff}}& K_{\mathrm{fs}}\\ K_{\mathrm{sf}}& K_{\mathrm{ss}}\end{array}\right))\left(\begin{array}{c}{u}_f\left(\omega \right)\\ u_s\left(\omega \right)\end{array}\right)=\left(\begin{array}{c}0\\ q_s\end{array}\right)---\left(8\right)$

其中u_s(ω)为已知的位移激励，q_s为待定的约束反力，即强迫位移产生的激 励力，无约束自由度u_f(ω)通过公式(8)的第一式，上半部分得出方程的第一式 解出：

(-ω²M_ff+iωB_ff+K_ff)u_f(ω)＝-(-ω²M_fs+iωB_fs+K_fs)u_s(ω)   (9)

位移产生的激励的约束力形式表达通过公式(8)的第二式，下半部分下式 求出。

q_s＝(-ω²M_sf+iωB_sf+K_sf)u(ω)_f+(-ω²M_ss+iωB_ss+K_ss)u(ω)_s   (10)

当外部激励载荷为速度时，经傅氏变换的表达式为：

$\stackrel{·}{u}\left(t\right)=\mathrm{i\omega u}\left(\omega \right)·e^{\mathrm{i\omega t}}---\left(6\right)$

所述的公式(8)、(9)和(10)分别为：

公式8-1：$({-\omega }^2\left(\begin{array}{cc}{M}_{\mathrm{ff}}& M_{\mathrm{fs}}\\ M_{\mathrm{sf}}& M_{\mathrm{ss}}\end{array}\right)+\mathrm{i\omega }\left(\begin{array}{cc}{B}_{\mathrm{ff}}& B_{\mathrm{fs}}\\ B_{\mathrm{sf}}& B_{\mathrm{ss}}\end{array}\right)+\left(\begin{array}{cc}{K}_{\mathrm{ff}}& K_{\mathrm{fs}}\\ K_{\mathrm{sf}}& K_{\mathrm{ss}}\end{array}\right))\left(\begin{array}{c}\left(\frac{1}{\mathrm{i\omega }}\right)u_f\left(\omega \right)\\ \left(\frac{1}{\mathrm{i\omega }}\right)u_s\left(\omega \right)\end{array}\right)=\left(\begin{array}{c}0\\ q_s\end{array}\right)$

公式9-1：

$(-{\omega }^2{M}_{\mathrm{ff}}+\mathrm{i\omega }{B}_{\mathrm{ff}}+K_{\mathrm{ff}})\left(\frac{1}{\mathrm{i\omega }}\right)u_f\left(\omega \right)=-(-{\omega }^2{M}_{\mathrm{fs}}+\mathrm{i\omega }{B}_{\mathrm{fs}}+K_{\mathrm{fs}})\left(\frac{1}{\mathrm{i\omega }}\right)u_s\left(\omega \right)$

公式10-1：

$\left(\begin{array}{c}{q}_s=\\ (-{\omega }^2{M}_{\mathrm{sf}}+\mathrm{i\omega }{B}_{\mathrm{sf}}+K_{\mathrm{sf}})\left(\frac{1}{\mathrm{i\omega }}\right)u{\left(\omega \right)}_f+(-{\omega }^2{M}_{\mathrm{ss}}+\mathrm{i\omega }{B}_{\mathrm{ss}}+K_{\mathrm{ss}})\left(\frac{1}{\mathrm{i\omega }}\right)u{\left(\omega \right)}_s;\end{array}\right)$

当外部激励载荷为加速度时，经傅氏变换后的表达式为：

$\stackrel{··}{u}\left(t\right)=-{\omega }^{2}u\left(\omega \right)·e^{\mathrm{i\omega t}}---\left(7\right)$

所述的公式(8)、(9)和(10)分别为：

公式8-2：$({-\omega }^2\left(\begin{array}{cc}{M}_{\mathrm{ff}}& M_{\mathrm{fs}}\\ M_{\mathrm{sf}}& M_{\mathrm{ss}}\end{array}\right)+\mathrm{i\omega }\left(\begin{array}{cc}{B}_{\mathrm{ff}}& B_{\mathrm{fs}}\\ B_{\mathrm{sf}}& B_{\mathrm{ss}}\end{array}\right)+\left(\begin{array}{cc}{K}_{\mathrm{ff}}& K_{\mathrm{fs}}\\ K_{\mathrm{sf}}& K_{\mathrm{ss}}\end{array}\right))\left(\begin{array}{c}-\left(\frac{1}{\mathrm{i\omega }}\right)u_f\left(\omega \right)\\ -\left(\frac{1}{\mathrm{i\omega }}\right)u_s\left(\omega \right)\end{array}\right)=\left(\begin{array}{c}0\\ q_s\end{array}\right)$

公式9-2：

$(-{\omega }^2{M}_{\mathrm{ff}}+\mathrm{i\omega }{B}_{\mathrm{ff}}+K_{\mathrm{ff}})\left(\frac{1}{i{\omega }^2}\right)u_f\left(\omega \right)=-(-{\omega }^2{M}_{\mathrm{fs}}+\mathrm{i\omega }{B}_{\mathrm{fs}}+K_{\mathrm{fs}})\left(\frac{1}{i{\omega }^2}\right)u_s\left(\omega \right)$

公式10-2：

$\left(\begin{array}{c}{q}_s=(-{\omega }^2{M}_{\mathrm{sf}}+\mathrm{i\omega }{B}_{\mathrm{sf}}+K_{\mathrm{sf}})(-\frac{1}{i{\omega }^2})u{\left(\omega \right)}_f+(-{\omega }^2{M}_{\mathrm{ss}}+\mathrm{i\omega }{B}_{\mathrm{ss}}+\\ K_{\mathrm{ss}})(-\frac{1}{i{\omega }^2})u{\left(\omega \right)}_s.\end{array}\right)$

所述扫频计算，获得该外部激励载荷下的节点位移传递函数的步骤具体如 下：所述的节点为有限元模型焊缝处焊线上的点；

当输入的激励载荷为简谐力时，在所述的公式(4)中输入单位幅值的简谐 力，得到力—位移输出系统的传递函数H_disp(ω)如下：

$H_{\mathrm{disp}}\left(\omega \right)=\frac{u\left(\omega \right)}{p\left(\omega \right)}=\frac{1}{-{\omega }^{2}M+\mathrm{i\omega B}+K}---\left(11\right).$

当输入为单位力时，公式(11)的输出结果即为力—位移传递函数H_disp(ω)；

通过力-位移传递函数，通过节点力计算膜应力和弯曲应力，得到通过计算 焊缝处等效结构应力的传递函数的过程如下：

F'_e(ω)＝B^-1K^eH_disp(ω)

F_e(ω)＝B^-1F_e′(ω)＝B^TK^eB^-1H_disp(ω)_f

[F(ω)]＝[N][F^e(ω)]

其中，ω是频率，K^e单元局部坐标系下的单元刚度矩阵，F_e′(ω)为单元局 部坐标系下的节点力，B为从系统坐标向单元局部坐标转换矩阵，是个常系数 矩阵，[N]为合成矩阵，F_e(ω)为系统坐标系下的节点力；

系统坐标系(x,y,z)下求解的焊趾处节点力矩阵{F^e}i需要

{F(ω)}_i＝{F_ix(ω),F_iy(ω),F_iz(ω),M_ix(ω),M_iy(ω),M_iz(ω)…}   (12)

i＝1,2,3,…n，F_ix其中F代表力，i代表节点号，x，y，z在代表平行于全局坐 标轴的力，M代表力矩；

{F′(ω)}_i＝{T}_i{F(ω)}_i   (13)

{f_{iy`</sub>(ω)}<sup>T</sup>＝{F<sub>iy</sub>′<sub>`}(ω)}^TL^-1   (14)

{m_{ix`</sub>(ω)}<sup>T</sup>＝{m<sub>i</sub>′<sub>x`}(ω}^TL^-1   (15)

f代表线力，就是节点力平均到焊线；

结构应力的传递函数为：

${\sigma }_s\left(\omega \right)={\sigma }_m\left(\omega \right)+{\sigma }_b\left(\omega \right)=\frac{f_{{\mathrm{ty}}^{`}}\left(\omega \right)}{d}+\frac{6m_{{\mathrm{ix}}^{`}}\left(\omega \right)}{d^2}---\left(17\right)$

σ_m(ω)是膜应力，σ_b(ω)是弯曲应力，σ_s为结构应力；

等效结构应力幅△Ss

$S_S\left(\omega \right)=\frac{{\sigma }_S\left(\omega \right)}{d^{(2-m)/2m}·I{\left(r\left(\omega \right)\right)}^{1/m}}---\left(18\right)$

其中I(r)是弯曲度比r的无量纲函数，常数m＝3.6，d为板厚；

当输入单位简谐载荷时，获得等效结构应力的传递函数为：

$H_i^s\left(\omega \right)=S_s\left(\omega \right)---\left(19\right).$

多载荷同时作用下的等效结构应力互功率谱的计算应用：

${\mathrm{PSD}}_s\left(f\right)=\underset{i,j=1}{\overset{i,j=n}{\Sigma }}{H}_i^s\left(f\right)H_j^s{\left(f\right)}^{*}{G}_{\mathrm{ij}}\left(f\right)---\left(23\right)$

其中i和j分别表示两个不同的载荷点输入的载荷作用下的等效结构应力 传递函数，相乘代表他们之间的耦合响应；Gi_j(f)是实测输入载荷的功率谱统计 结果与前面的S_S(ω)截然不同的，i和j分别代表两个不同输入载荷之间的互功 率谱统计，如i和j取相同值时，G_ii(f)就是载荷i的自功率谱统计。

为了计算结构应力变化范围的概率密度函数，需要用到PSD矩函数，其定 义如下：

$m_n=\underset{0}{\overset{\infty }{\int }}{f}^n·\mathrm{PS}{D}_s\left(f\right)\mathrm{df}=\Sigma f^n·{\mathrm{PSD}}_s\left(f\right)·\mathrm{\delta f}---\left(24\right)$

其中PSD_S(f)为等效结构应力单边功率谱密度函数；

基于Dirlik法的统计效结构应力的概率密度函数；

$p{\left(S\right)}_D=\frac{\frac{D_1}{Q}·e^{\frac{-z}{Q}}+\frac{D_2·Z}{R^2}·e^{\frac{{-Z}^2}{2·R^2}}+D_3·Z·e^{\frac{-Z^2}{2}}}{2·\sqrt{m_0}}---\left(25\right)$

其中：

单位时间内峰值点统计数单位时间内过零点统计数单位时间内峰值数与过零点数的比例关系为其它中间变量表达式如下：

$R=\frac{\gamma -x_m-D_1^2}{1-\gamma -D_1+D_1^2},$$D_1=\frac{2·(x_m-{\gamma }^2)}{1+{\gamma }^2},$$x_m=\frac{m_1}{m_0}·\sqrt{\frac{m_2}{m_4}},$$z=\frac{S}{2·\sqrt{m_0}}$

$Q=\frac{1.25·(\gamma -D_3-D_2·R)}{D_1},$$D_2=\frac{1-\gamma -D_1+D_1^2}{1-R},$D₃＝1-D₁-D₂

单位时间内等效结构应力变化范围及频次统计

n_i(S_i)＝p(S_i)dS   (26)

单位时间内疲劳损伤统计及累积

$E\left[D\right]=\underset{i}{\Sigma }\frac{n_i\left(S_i\right)}{N\left(S_i\right)}=\frac{\underset{0}{\overset{\infty }{\int }}{S}^{1/h}P\left(S\right)\mathrm{dS}}{C_d^{1/h}}---\left(27\right)$

其中Cd，h为材料相关的主S-N曲线参数[11]，当损伤达到1时，疲劳寿命 结束[12]，疲劳寿命结果为时间(单位：秒)，表明该结构在这种振动条件下能 够存活的时间：

附图说明

为了更清楚的说明本发明的实施例或现有技术的技术方案，下面将对实施 例或现有技术描述中所需要使用的附图做一简单地介绍，显而易见地，下面描 述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不 付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图

图2-1为实施例1中T型焊接接头的有限元模型一示意图

图2-2为实施例1中T型接头的第一阶模态示意图

图2-3为实施例1垂向载荷工况示意图

图2-4为实施例1单位载荷结构应力的对比示意图

图2-5为实施例1单位载荷等效结构应力对比示意图

图2-6为实施例1实测力载荷时间里程示意图

图2-7为实施例1中随机振动理论统计的寿命与准静态法估算的理论的对比 示意图

图3-1为实施例2中T型焊接接头的有限元模型二示意图

图3-2为实施例2中T型接头的第一阶模态示意图

图3-3为实施例2中实测力载荷时间的历程示意图

图3-4为实施例2中实测力载荷功率谱密度示意图

图3-5为实施例2中焊接点序列的结构应力传递函数示意图

图3-6为实施例2中焊接结构中心点处应力传递函数

图3-7为实施例2中在频率为8.55Hz单位载荷下焊缝结构中心点处结构应 力与静态单位载荷下对比示意图

图3-8为实施例2中焊缝结点序列的等效结构应力传递函数示意图

图3-9为实施例2中焊缝结点序列的等效结构应力响应功率密度示意图

图3-10为实施例2焊缝结点序列的等效结构应力概率密度示意图

图3-11为实施例2中焊缝结点序列的疲劳寿命对比示意图

具体实施方式

为使本发明的实施例的目的、技术方案和优点更加清楚，下面结合本发明 实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：

如图1所示：一种复杂焊接结构随机振动疲劳寿命预测方法，包括如下步 骤：

首先要建立包含焊缝细节的完整系统有限元模型，因为考虑了焊缝细节的 模型才能反映焊接结构的应力集中，模型应该尽量贴近实际结构。其中焊缝局 部结构可以根据具体需要使用典型的8节点六面体实体单元或四节点四边形薄 壳单元模拟。

然后，确定系统的边界条件，在本发明的实施例中边界条件重点描述力、 位移、速度和加速度的外部激励载荷，边界条件中的连接关系，本领域技术人 员可根据实际情况自行调节，在这里不再赘述。

确定边界条件后，建立动力学方程：

$\left[M\right]\left\{\stackrel{··}{x}\left(t\right)\right\}+\left[B\right]\left\{\stackrel{·}{x}\left(t\right)\right\}+\left[K\right]\left\{x\left(t\right)\right\}=\left\{f\left(t\right)\right\}---\left(1\right)$

其中[M]为系统的质量矩阵，[B]为系统的阻尼矩阵，[K]为系统的刚度矩阵， {x(t)}为系统的位移向量，为系统的速度矢量，为系统的加速度矢量， {f(t)}表示施加的外部激励载荷为力。

然后将所述的边界条件引入所述的动力学方程。求解动力学方程方法有多 种，由于工程结构自由度庞大，载荷一般为随机载荷，且时间周期比较长，这 就给动力学方程的求解带来困难，所以目前能够运用到解决实际工程问题的方 法主要有的基于时域的模态叠加法和频域的功率谱密度法。

由于时域的模态叠加法能够缩减动力学方程自由度，但随机的外载荷很难 用固定的函数来描述，结构所受的随机载荷作用时间一般很长，复杂结构的时 域法会因需要计算的时间步态太多，计算时间太长、结果文件太大等困难导致 在现有计算机条件下无法进行，且模态的截断会导致焊接结构局部应力集中计 算精度下降，因此本专利采用频域的功率谱密度法来计算。

当施加的激励为力载荷时，经傅氏变换到频域的表达式为：

f(t)＝p(ω)·e^iωt   (2)

将公式(2)带入动力学方程公式(1)，经傅氏变换到频域后，方程为：

-ω²[M]{u(ω)}e^iωt+iω[B]{u(ω)}e^iωt+[K]{u(ω)}e^iωt＝{P(ω)}e^iωt   (3)

方程两边抵消掉复指数e^iωt，得到公式4：

[-ω²M+iωB+K]{u(ω)}＝{p(ω)}   (4)

考虑到动力学方程，即公式1是受力平衡方程，方程的{f(t)}为力载荷，对 于位移，速度和加速度等的位移驱动载荷，在引入动力学方程时，需要转化成 力载荷。当施加的激励载荷为位移时，位移经过傅氏变换到频域的表达式为：

u(t)＝u(ω)·e^iωt   (5)

所动力学方程，公式(3)中的自由度{u(ω)}，按位移激励，进行重组，焊 接系统的质量阵、阻尼阵和刚度阵按自由度分块为约束自由度，下标为s和无约 束自由度，下标为f，如式(8)所示：

$({-\omega }^2\left(\begin{array}{cc}{M}_{\mathrm{ff}}& M_{\mathrm{fs}}\\ M_{\mathrm{sf}}& M_{\mathrm{ss}}\end{array}\right)+\mathrm{i\omega }\left(\begin{array}{cc}{B}_{\mathrm{ff}}& B_{\mathrm{fs}}\\ B_{\mathrm{sf}}& B_{\mathrm{ss}}\end{array}\right)+\left(\begin{array}{cc}{K}_{\mathrm{ff}}& K_{\mathrm{fs}}\\ K_{\mathrm{sf}}& K_{\mathrm{ss}}\end{array}\right))\left(\begin{array}{c}{u}_f\left(\omega \right)\\ u_s\left(\omega \right)\end{array}\right)=\left(\begin{array}{c}0\\ q_s\end{array}\right)---\left(8\right)$

其中u_s(ω)为已知的位移激励，q_s为待定的约束反力，即强迫位移产生的激 励力，无约束自由度u_f(ω)通过公式(8)的第一式，上半部分得出方程的第一式 解出：(-ω²M_ff+iωB_ff+K_ff)u_f(ω)＝-(-ω²M_fs+iωB_fs+K_fs)u_s(ω)   (9)

位移产生的激励的约束力形式表达通过公式(8)的第二式，即下半部分下 式求出。q_s＝(-ω²M_sf+iωB_sf+K_sf)u(ω)_f+(-ω²M_ss+iωB_ss+K_ss)u(ω)_s   (10)

以上的过程实际上是把有限元模型中存在位移激励的节点自由度和没有位 移激励的自由度分开。

更进一步的，当外部激励载荷为速度或加速度时，傅氏变换的结果如下：

当外部激励载荷为速度时，经傅氏变换的表达式为：

$\stackrel{·}{u}\left(t\right)=\mathrm{i\omega u}\left(\omega \right)·e^{\mathrm{i\omega t}}---\left(6\right)$

所述的公式(8)、(9)和(10)分别为：

公式8-1：$({-\omega }^2\left(\begin{array}{cc}{M}_{\mathrm{ff}}& M_{\mathrm{fs}}\\ M_{\mathrm{sf}}& M_{\mathrm{ss}}\end{array}\right)+\mathrm{i\omega }\left(\begin{array}{cc}{B}_{\mathrm{ff}}& B_{\mathrm{fs}}\\ B_{\mathrm{sf}}& B_{\mathrm{ss}}\end{array}\right)+\left(\begin{array}{cc}{K}_{\mathrm{ff}}& K_{\mathrm{fs}}\\ K_{\mathrm{sf}}& K_{\mathrm{ss}}\end{array}\right))\left(\begin{array}{c}\left(\frac{1}{\mathrm{i\omega }}\right)u_f\left(\omega \right)\\ \left(\frac{1}{\mathrm{i\omega }}\right)u_s\left(\omega \right)\end{array}\right)=\left(\begin{array}{c}0\\ q_s\end{array}\right)$

公式9-1：

$(-{\omega }^2{M}_{\mathrm{ff}}+\mathrm{i\omega }{B}_{\mathrm{ff}}+K_{\mathrm{ff}})\left(\frac{1}{\mathrm{i\omega }}\right)u_f\left(\omega \right)=-(-{\omega }^2{M}_{\mathrm{fs}}+\mathrm{i\omega }{B}_{\mathrm{fs}}+K_{\mathrm{fs}})\left(\frac{1}{\mathrm{i\omega }}\right)u_s\left(\omega \right)$

公式10-1：

$\left(\begin{array}{c}{q}_s=\\ (-{\omega }^2{M}_{\mathrm{sf}}+\mathrm{i\omega }{B}_{\mathrm{sf}}+K_{\mathrm{sf}})\left(\frac{1}{\mathrm{i\omega }}\right)u{\left(\omega \right)}_f+(-{\omega }^2{M}_{\mathrm{ss}}+\mathrm{i\omega }{B}_{\mathrm{ss}}+K_{\mathrm{ss}})\left(\frac{1}{\mathrm{i\omega }}\right)u{\left(\omega \right)}_s;\end{array}\right)$

当外部激励载荷为加速度时，经傅氏变换后的表达式为：

$\stackrel{··}{u}\left(t\right)=-{\omega }^{2}u\left(\omega \right)·e^{\mathrm{i\omega t}}---\left(7\right)$

所述的公式(8)、(9)和(10)分别为：

公式8-2：$({-\omega }^2\left(\begin{array}{cc}{M}_{\mathrm{ff}}& M_{\mathrm{fs}}\\ M_{\mathrm{sf}}& M_{\mathrm{ss}}\end{array}\right)+\mathrm{i\omega }\left(\begin{array}{cc}{B}_{\mathrm{ff}}& B_{\mathrm{fs}}\\ B_{\mathrm{sf}}& B_{\mathrm{ss}}\end{array}\right)+\left(\begin{array}{cc}{K}_{\mathrm{ff}}& K_{\mathrm{fs}}\\ K_{\mathrm{sf}}& K_{\mathrm{ss}}\end{array}\right))\left(\begin{array}{c}-\left(\frac{1}{i{\omega }^2}\right)u_f\left(\omega \right)\\ -\left(\frac{1}{i{\omega }^2}\right)u_s\left(\omega \right)\end{array}\right)=\left(\begin{array}{c}0\\ q_s\end{array}\right)$

公式9-2：

$(-{\omega }^2{M}_{\mathrm{ff}}+\mathrm{i\omega }{B}_{\mathrm{ff}}+K_{\mathrm{ff}})\left(\frac{1}{i{\omega }^2}\right)u_f\left(\omega \right)=-(-{\omega }^2{M}_{\mathrm{fs}}+\mathrm{i\omega }{B}_{\mathrm{fs}}+K_{\mathrm{fs}})(-\frac{1}{i{\omega }^2})u_s\left(\omega \right)$

公式10-2：

$\left(\begin{array}{c}{q}_s=(-{\omega }^2{M}_{\mathrm{sf}}+\mathrm{i\omega }{B}_{\mathrm{sf}}+K_{\mathrm{sf}})(-\frac{1}{i{\omega }^2})u{\left(\omega \right)}_f+(-{\omega }^2{M}_{\mathrm{ss}}+\mathrm{i\omega }{B}_{\mathrm{ss}}+\\ K_{\mathrm{ss}})(-\frac{1}{i{\omega }^2})u{\left(\omega \right)}_s.\end{array}\right)$

将对应的外部激励载荷引入动力学方程后，在焊接结构的载荷输入点分别 施加不同的外部激励载荷，进行扫频计算。

当输入的激励载荷为简谐力时，在所述的公式(4)中输入单位幅值的简谐 力，得到力—位移输出系统的传递函数H_disp(ω)如下：

$H_{\mathrm{disp}}\left(\omega \right)=\frac{u\left(\omega \right)}{p\left(\omega \right)}=\frac{1}{-{\omega }^{2}M+\mathrm{i\omega B}+K}---\left(11\right).$

当输入为单位力时，公式(11)的输出结果即为力—位移传递函数H_disp(ω)。

力-位移传递函数，通过节点力计算膜应力和弯曲应力，得到通过计算焊缝 处等效结构应力的传递函数的过程如下：

F'_e(ω)＝B^-1K^eH_disp(ω)

F_e(ω)＝B^-1F_e′(ω)＝B^TK^eB^-1H_disp(ω)_f

[F(ω)]＝[N][F^e(ω)]

基本标准公式是不带参数的是一个结点位移与结构应力之间的关系式，而 本计算，每一步都是与频率参数相关的，并且本计算的所有的结点力和结构应 力都是复数计算(即都包含实部和虚数两个部分)，这些都与传统公式不同。

如果焊线节点数为4个，则合成矩阵如下：

$N=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0\\ 0& 0& 1& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)$

如果焊线节点数为5个(节点数为n，则为n×6矩阵)，则合成矩阵如下：

其中，ω是频率，K^e单元局部坐标系下的单元刚度矩阵，F_e′(ω)为单元局 部坐标系下的节点力，B为从系统坐标向单元局部坐标转换矩阵，是个常系数 矩阵，[N]为合成矩阵，F_e(ω)为系统坐标系下的节点力；

系统坐标系(x,y,z)下求解的焊趾处节点力矩阵{F^e}_i需要

{F(ω)}_i＝{F_ix(ω),F_iy(ω),F_iz(ω),M_ix(ω),M_iy(ω),M_iz(ω)…}   (12)

i＝1,2,3,…n，F_ix其中F代表力，i代表节点号，x，y，z在代表平行于全局坐 标轴的力，M代表力矩；

{F′(ω)}_i＝{T}_i{F(ω)}_i   (13)

{f_{iy`</sub>(ω)}<sup>T</sup>＝{F′<sub>iy`}(ω)}^TL^-1   (14)

{m_{ix`</sub>(ω)}<sup>T</sup>＝{m′<sub>ix`}(ω}^TL^-1   (15)

f代表线力，就是节点力平均到焊线；

结构应力的传递函数为：

${\sigma }_s\left(\omega \right)={\sigma }_m\left(\omega \right)+{\sigma }_b\left(\omega \right)=\frac{f_{{\mathrm{ty}}^{`}}\left(\omega \right)}{d}+\frac{6m_{{\mathrm{ix}}^{`}}\left(\omega \right)}{d^2}---\left(17\right)$

σ_m(ω)是膜应力，σ_b(ω)是弯曲应力，σ_s为结构应力；

等效结构应力幅△Ss

$S_S\left(\omega \right)=\frac{{\sigma }_S\left(\omega \right)}{d^{(2-m)/2m}·I{\left(r\left(\omega \right)\right)}^{1/m}}---\left(18\right)$

其中I(r)是弯曲度比r的无量纲函数，常数m＝3.6，d为板厚；

当输入单位简谐载荷时，获得等效结构应力的传递函数为：

$H_i^s\left(\omega \right)=S_s\left(\omega \right)---\left(19\right).$

所述的对每个实测的载荷输入进行傅氏变换，获得每个外部激励载荷的功 率谱和载荷之间的互功率谱的方法如下：

以两个载荷输入为例，首先计算输入载荷p_i(t)和p_j(t)的互相关函数:

$R_{\mathrm{ij}}\left(\tau \right)=\underset{\tau \to \infty }{\lim }\frac{1}{T}{\int }_0^T{p}_i\left(t\right)p_j(t-\tau )\mathrm{dt}---\left(20\right)$

载荷p_i(t)和p_j(t)之间的的互功率谱密度函数为:

$S_{\mathrm{ij}}\left(\omega \right)=\underset{\tau \to \infty }{\lim }\frac{2}{T}\left({\int }_0^T{p}_a\left(t\right)e^{-\mathrm{i\omega t}}\mathrm{dt}\right)\left({\int }_0^T{p}_b\left(t\right)e^{\mathrm{i\omega t}}\mathrm{dt}\right)---\left(21\right)$

其中为S_ij(ω)的复共轭函数；当i,j相同时公式(21)计算的结果为自功 率谱密度函数；当i,j不同时公式(21)计算的结果为互功率谱密度函数。

多载荷同时作用下的等效结构应力互功率谱的计算应用：

${\mathrm{PSD}}_s\left(f\right)=\underset{i,j=1}{\overset{i,j=n}{\Sigma }}{H}_i^s\left(f\right)H_j^s{\left(f\right)}^{*}{G}_{\mathrm{ij}}\left(f\right)---\left(23\right)$

其中i和j分别表示两个不同的载荷点输入的载荷作用下的等效结构应力 传递函数，相乘代表他们之间的耦合响应；G_ij(f)是实测输入载荷的功率谱统计 结果与前面的S_S(ω)截然不同的，i和j分别代表两个不同输入载荷之间的互功 率谱统计，如i和j取相同值时，G_ii(f)就是载荷i的自功率谱统计。

为了计算结构应力变化范围的概率密度函数，需要用到PSD矩函数，其定 义如下：

$m_n=\underset{0}{\overset{\infty }{\int }}{f}^n·\mathrm{PS}{D}_s\left(f\right)\mathrm{df}=\Sigma f^n·{\mathrm{PSD}}_s\left(f\right)·\mathrm{\delta f}---\left(24\right)$

其中PSD_S(f)为等效结构应力单边功率谱密度函数；

由于Dirlik法可以直接从应力响应功率谱统计获得应力变化范围的概率密 度函数，避免了传统把频域结果经过傅氏逆变换到时域的处理过程中人为添 加一些伪随机数来填补功率谱密度缺失的相位信息对疲劳预测结果的影响， 所以采用Dirlik法来统计等效结构应力。

基于Dirlik法的统计效结构应力的概率密度函数；

$p{\left(S\right)}_D=\frac{\frac{D_1}{Q}·e^{\frac{-z}{Q}}+\frac{D_2·Z}{R^2}·e^{\frac{{-Z}^2}{2·R^2}}+D_3·Z·e^{\frac{-Z^2}{2}}}{2·\sqrt{m_0}}---\left(25\right)$

其中：

单位时间内峰值点统计数单位时间内过零点统计数单位时间内峰值数与过零点数的比例关系为其它中间变量表达式如下：

$R=\frac{\gamma -x_m-D_1^2}{1-\gamma -D_1+D_1^2},$$D_1=\frac{2·(x_m-{\gamma }^2)}{1+{\gamma }^2},$$x_m=\frac{m_1}{m_0}·\sqrt{\frac{m_2}{m_4}},$$z=\frac{S}{2·\sqrt{m_0}}$

$Q=\frac{1.25·(\gamma -D_3-D_2·R)}{D_1},$$D_2=\frac{1-\gamma -D_1+D_1^2}{1-R},$D₃＝1-D₁-D₂

单位时间内等效结构应力变化范围及频次统计

n_i(S_i)＝p(S_i)dS   (26)

单位时间内疲劳损伤统计及累积

$E\left[D\right]=\underset{i}{\Sigma }\frac{n_i\left(S_i\right)}{N\left(S_i\right)}=\frac{\underset{0}{\overset{\infty }{\int }}{S}^{1/h}P\left(S\right)\mathrm{dS}}{C_d^{1/h}}---\left(27\right)$

其中Cd，h为材料相关的主S-N曲线参数[11]，当损伤达到1时，疲劳寿命 结束，疲劳寿命结果为时间(单位：秒)，表明该结构在这种振动条件下能够存 活的时间：

实施例1，准静态算例验证

如图2-1所示：T型焊接接头结构一，该结构的最低一阶模态的频率为 479.4HZ，如图2-2所示，而激扰力载荷频率范围为(0-40HZ)，激励频率与固 有频率相差较远，结构不会发生共振。分别用准静态法siatic和随机振动法 random计算该焊接结构的疲劳寿命，准静态法的计算结构应该与该随机振动疲 劳算法的结果理论上应该相近。

施加如图2-3的垂向单位力(牛顿)作用，上边与下边施加固定约束，利 用商用有限元软件分别计算静态和垂向单位力扫频的节点位移解，并相关公式 计算结构应力与等效结构应力。单位载荷结构应力的对比如图2-4所示：代表 准静态法的static曲线与代表随机振动法的红色和绿色曲线基本重合。如图2-5 所示：等效结构应力的对比结果，两种方法也基本相同。图2-7为实测力载荷 的时间历程示意图。如图2-7所示：该专利新方法和传统方法计算的焊缝的寿 命结果非常接近。

实施例2动态算例验证

为了突显随机振动疲劳算法的优点，特设计载荷频率于固有频率存在交集 的小算例(如图3-1T型焊接接头结构二)，该结构的最低一阶模态的频率为 8.55Hz，如图3-2～3-4所示，而激扰力载荷频率范围为(0-40HZ)，激励频率 在8.55Hz处有较多的能量分布(见图3-5)，结构将会发生比较强烈共振。分 别用准静态法和随机振动法计算该焊接结构的疲劳寿命，理论上随机振动疲劳 算法寿命预测结果要比准静态法的计算得寿命要短很多。

载荷的施加位置及大小和实施例一致。其功率谱密度见图3-4。

通过对比图3-6～图3-11可以看出，当激励载荷频率与结构最低阶固有频 率存在交集，结构将会发生较强烈共振，结构应力水平在固有频率附近急剧增 加，寿命结果下降近百倍，此时用准静态算法计算就会导致结果与实际偏差较 大，同时也证明了该焊接结构随机振动疲劳寿命预测法的优点。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局 限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本 发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护 范围之内。