CN109344433B - 基于响应信号的灵敏度数值计算方法 - Google Patents

基于响应信号的灵敏度数值计算方法 Download PDF

Info

Publication number
CN109344433B
CN109344433B CN201810981540.0A CN201810981540A CN109344433B CN 109344433 B CN109344433 B CN 109344433B CN 201810981540 A CN201810981540 A CN 201810981540A CN 109344433 B CN109344433 B CN 109344433B
Authority
CN
China
Prior art keywords
frequency
response function
sensitivity
matrix
rigidity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810981540.0A
Other languages
English (en)
Other versions
CN109344433A (zh
Inventor
费庆国
朱锐
姜东�
陈素芳
杭晓晨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southeast University
Original Assignee
Southeast University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southeast University filed Critical Southeast University
Priority to CN201810981540.0A priority Critical patent/CN109344433B/zh
Publication of CN109344433A publication Critical patent/CN109344433A/zh
Application granted granted Critical
Publication of CN109344433B publication Critical patent/CN109344433B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Complex Calculations (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明提供了一种基于响应信号的灵敏度数值计算方法，构造速度频响函数矩阵，并获得前m阶模态频率，从结构第一个节点开始添加刚度摄动项；基于速度频响函数信息代入矩阵修正公式形式获得摄动后的速度频响函数；提取结构的频率信息，获得结构模态频率对刚度的灵敏度；按照节点顺序改变刚度摄动点位置，从而获得整个结构模态频率对刚度的灵敏度，绘制灵敏度曲线。本发明方法首先通过有限元计算获得结构的速度频响函数信息，当结构刚度发生变化时，利用矩阵变换公式无需进行有限元再次计算，只需要初始的速度频响函数信息进行数值计算即可获得摄动后的速度响应信息，简化计算效率，更加方便，基于速度频响函数实现了频率对刚度的灵敏度快速计算方法，具有实际工程意义。

Description

基于响应信号的灵敏度数值计算方法
技术领域
本发明涉及一种灵敏度分析方法，具体涉及一种结构灵敏度数值计算方法。
背景技术
灵敏度分析是一种评价因设计变量或参数的改变而引起结构响应特性变化率的方法。结构系统灵敏度的研究是一个很特别的领域，它是当前计算力学和结构工程领域的主要研究方向之一。
实际应用中，结构的灵敏度分析在结构优化、可靠性评估和参数识别等占据重要作用，方法按照策略可以分为离散法和变分法。传统的灵敏度分析方法一般需要在摄动下，重新进行有限元计算，特别当结构较大，有限元数量较多时，计算量大，效率低。有些工况需要测量速度，从而基于速度频响函数信息进行如对刚度参数等的灵敏度分析，如果解决传统灵敏度计算效率低，已成为亟待解决的实际工程问题。
发明内容
发明目的：本发明的目的在于针对现有技术的不足，提供一种基于响应信号的灵敏度数值计算方法。
技术方案：本发明提供了一种基于响应信号的灵敏度数值计算方法，包括以下步骤：
(1)基于Matlab软件进行有限元分析，获得速度频响函数，构造频响函数矩阵，并获得前m阶模态频率，从结构第一个节点开始添加刚度摄动项；
(2)基于步骤(1)中速度频响函数矩阵代入矩阵修正公式获得摄动后的速度频响函数矩阵；
(3)提取结构的频率信息，获得结构模态频率对刚度的灵敏度；
(4)按照节点顺序改变刚度摄动点位置，重复步骤(2)(3)获得对应的灵敏度，从而获得整个结构模态频率对刚度的灵敏度，绘制灵敏度曲线。
进一步，步骤(1)包括以下步骤：
(11)结构的自由度为N，获得结构的速度频响函数矩阵为：
其中，sij表示在结构节点j作用单位脉冲下节点i的速度频响函数，i＝1，2…N，j＝1，2…N；
(12)在公式(1)矩阵中，按照从左到右、从上到下的原则做sij幅频图，判断该曲线局部极大值的个数为T，若T小于m，则继续作图；若T等于m，则终止作图；
(13)以(12)中筛选的幅频图为例，该幅频图的局部极大值的横坐标即为频率值，分别为fr(r＝1，2…m)，将fr元素组成m维列向量f。
进一步，步骤(2)包括以下步骤：
(21)在结构i节点处添加刚度摄动量Δki，此时结构的速度频响函数矩阵为SV *，根据速度频响函数矩阵的逆矩阵与动态刚度矩阵的关系，推导获得：
其中，ui∈RN×1，表示在列向量第i个元素为e为虚数单位，即e2＝-1，ω为圆频率；
(22)由矩阵修正公式可知：
即可以建立摄动前后结构的速度响应的联系；
(23)由公式(1)、(2)、(3)联列，化解可以获得刚度摄动后速度频响函数矩阵与初始频响函数矩阵的关系：
进一步，步骤(3)包括以下步骤：
(31)基于摄动后的速度频响函数矩阵重复步骤(12)(13)提取摄动后结构的频率(r＝1，2…m),将元素组成m维列向量f*；
(32)定义获取灵敏度函数gi，表示在节点i处发生刚度摄动后的灵敏度，函数输入量为摄动前后的频率，具体公式如下：
gi(fr *,fr)＝[li,dis] (5)
其中，等式左边括号量即为输入函数gi的输入量，等式右边即为计算输出结果，具体形式为：
进一步确定灵敏度最大处的阶次信息，调用MATALB中find函数：
有益效果：本发明方法首先通过有限元计算获得结构的速度频响函数信息，当结构刚度发生变化时，利用矩阵变换公式无需进行有限元再次计算，只需要初始的速度频响函数信息进行数值计算即可获得摄动后的速度响应信息，简化计算效率，更加方便，基于速度频响函数实现了频率对刚度的灵敏度快速计算方法，具有实际工程意义。
附图说明
图1为实施例中8自由度弹簧-阻尼-质量系统示意图；
图2为实施例中8自由度弹簧-阻尼-质量系统在节点1处加刚度摄动量后的示意图；
图3为系统结构的速度频响函数s11曲线；
图4为系统结构在节点1处添加摄动量后的速度频响函数曲线；
图5基于速度响应信号结构频率对刚度的灵敏度曲线。
具体实施方式
下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
本实施例采用一个简单的8自由度弹簧-阻尼-质量系统来验证，如图1所示，系统的参数分别为：mi＝1kg(i＝1,2…8)，弹簧ki＝10N/m(i＝1,2…9)，阻尼器ci＝0.02Nm/s(i＝1,2…9，添加的刚度摄动量为Δk＝2N/m。具体包括以下步骤：
步骤1，基于Matlab软件进行有限元分析，获得速度频响函数，构造频响函数矩阵，并获得前8阶模态频率，从结构第一个节点开始添加刚度摄动项：
1.1)结构的自由度为8，获得结构的速度频响函数矩阵为：
1.2)在公式(1)矩阵中，按照从左到右，从上到下的原则做sij幅频图，判断该曲线局部最大值的个数为T，若T小于8，则继续作图；若T等于8，则终止作图；首先作图s11，如图3所示，可知共有8个局部最大值，满足上面要求，故终止作图；
1.3)以1.2)中筛选的幅频图s11为例，该幅频图的局部极值的横坐标即为频率值，分别为fr(r＝1，2…8)，将fr元素组成8维列向量f，具体为：
f1＝0.174,f2＝0.344,f3＝0.504,f4＝0.646
f5＝0.772,f6＝0.874,f7＝0.946,f8＝0.994
f＝[0.174 0.344 0.504 0.646 0.772 0.874 0.946 0.994]T
在结构节点1处处添加刚度摄动量Δk1＝2，如图2所示。
步骤2，基于步骤(1)中速度频响函数矩阵代入矩阵变换公式获得摄动后的速度响应信息：
2.1)此时结构的速度频响函数为根据速度频响函数矩阵的逆矩阵与动态刚度矩阵的关系，推导获得：
η＝1,Δki＞0
e为虚数单位，即e2＝-1，ω为圆频率；
2.2)由矩阵修正公式可知：
即可以建立其摄动前后结构的速度响应的联系；
2.3)由公式(1)、(2)、(3)联列化解，可以获得刚度摄动后速度频响函数矩阵与初始频响函数矩阵的关系：
由于等式右边均为已知项，故可以获得摄动后的速度频响函数矩阵
步骤3，提取结构的频率信息，获得结构模态频率对刚度的灵敏度：
3.1)基于摄动后的速度频响函数矩阵，重复步骤(1.2)(1.3)可以作图见图4，由图4可以提取摄动后结构的频率(r＝1，2…8),将元素组成8维列向量f*，具体值如下：
f*＝[0.178 0.35 0.51 0.654 0.778 0.874 0.948 0.996]T
3.2)定义获取灵敏度函数g1，函数输入量为摄动前后的频率，具体公式如下：
其中，等式左边括号量即为输入函数gi的输入量，等式右边即为计算输出结果，具体形式为
进一步确定灵敏度最大处的阶次信息，调用MATALB中find函数
步骤4，按照节点顺序改变刚度摄动点位置，重复步骤(2)(3)获得对应得灵敏度，从而获得整个结构模态频率对刚度的灵敏度，绘制灵敏度曲线，见图5。
图5结果表明，当结构刚度发生改变时，在节点5处灵敏度的最大，故在结构节点5处修改刚度影响最大。本发明的快速灵敏度分析方法，突破了传统灵敏度需要进行多次的计算的局限性，只需要进行一次有限元计算，利用初始结构的速度响应信号进行数值计算获得刚度摄动后的结构的速度响应，更加快捷。

Claims (1)

1.一种基于响应信号的灵敏度数值计算方法，其特征在于：包括以下步骤：
(1)基于Matlab软件进行有限元分析，获得速度频响函数，构造频响函数矩阵，并获得前m阶模态频率，从结构第一个节点开始添加刚度摄动项；
(2)基于步骤(1)中速度频响函数矩阵代入矩阵修正公式获得摄动后的速度频响函数矩阵；
(3)提取结构的频率信息，获得结构模态频率对刚度的灵敏度；
(4)按照节点顺序改变刚度摄动点位置，重复步骤(2)(3)获得对应的灵敏度，从而获得整个结构模态频率对刚度的灵敏度，绘制灵敏度曲线；
步骤(1)包括以下步骤：
(11)结构的自由度为N，获得结构的速度频响函数矩阵为：
其中，sij表示在结构节点j作用单位脉冲下节点i的速度频响函数，i＝1，2…N，j＝1，2…N；
(12)在公式(1)矩阵中，按照从左到右、从上到下的原则做sij幅频图，判断该曲线局部极大值的个数为T，若T小于m，则继续作图；若T等于m，则终止作图；
(13)以(12)中筛选的幅频图为例，该幅频图的局部极大值的横坐标即为频率值，分别为fr，r＝1，2…m，将fr元素组成m维列向量f；
步骤(2)包括以下步骤：
(21)在结构i节点处添加刚度摄动量Δki，此时结构的速度频响函数矩阵为根据速度频响函数矩阵的逆矩阵与动态刚度矩阵的关系，推导获得：
其中，ui∈RN×1，表示在列向量第i个元素为e为虚数单位，即e2＝-1，ω为圆频率；
(22)由矩阵修正公式可知：
即可以建立摄动前后结构的速度响应的联系；
(23)由公式(1)、(2)、(3)联列，化解可以获得刚度摄动后速度频响函数矩阵与初始频响函数矩阵的关系：
步骤(3)包括以下步骤：
(31)基于摄动后的速度频响函数矩阵重复步骤(12)(13)提取摄动后结构的频率r＝1，2…m，将元素组成m维列向量f*；
(32)定义获取灵敏度函数gi，表示在节点i处发生刚度摄动后的灵敏度，函数输入量为摄动前后的频率，具体公式如下：
其中，等式左边括号量即为输入函数gi的输入量，等式右边即为计算输出结果，具体形式为：
进一步确定灵敏度最大处的阶次信息，调用MATALB中find函数：
CN201810981540.0A 2018-08-27 2018-08-27 基于响应信号的灵敏度数值计算方法 Active CN109344433B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810981540.0A CN109344433B (zh) 2018-08-27 2018-08-27 基于响应信号的灵敏度数值计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810981540.0A CN109344433B (zh) 2018-08-27 2018-08-27 基于响应信号的灵敏度数值计算方法

Publications (2)

Publication Number Publication Date
CN109344433A CN109344433A (zh) 2019-02-15
CN109344433B true CN109344433B (zh) 2019-06-21

Family

ID=65291594

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810981540.0A Active CN109344433B (zh) 2018-08-27 2018-08-27 基于响应信号的灵敏度数值计算方法

Country Status (1)

Country Link
CN (1) CN109344433B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110569587B (zh) * 2019-08-29 2022-12-02 湖北工业大学 基于频响函数预估结构局部修改后动力学特性的方法
CN110569611B (zh) * 2019-09-12 2023-02-03 南京林业大学 一种基于多复变量法的结构频响函数灵敏度分析方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107220450B (zh) * 2017-06-08 2018-05-11 东南大学 一种非均匀材料连续分布力学参数场间接获取方法
CN107356387B (zh) * 2017-07-21 2018-06-19 东南大学 一种模态试验中多传感器附加质量消除方法
CN108256256B (zh) * 2018-01-31 2019-03-12 东南大学 一种保证阻尼系统特定模态频率不变的方法

Also Published As

Publication number Publication date
CN109344433A (zh) 2019-02-15

Similar Documents

Publication Publication Date Title
CN109344433B (zh) 基于响应信号的灵敏度数值计算方法
CN106709182B (zh) 一种地震作用下顺层岩质边坡稳定可靠性安全评价方法
CN107092759B (zh) 基于重力坝坝基参数反演的坝体位移监测点优化布置方法
CN102081359B (zh) 基于DSP Builder的变时滞超混沌数字电路设计方法及电路
Zhang et al. New results on $ H_ {\infty} $ filtering for fuzzy time-delay systems
CN107066476A (zh) 一种基于物品相似度的实时推荐方法
CN113722966B (zh) 一种集成电路板仿真多级分布式并行计算方法
CN103954295A (zh) 一种基于加速度传感器的计步方法
CN108108559B (zh) 一种基于子结构的结构响应获取方法及灵敏度获取方法
Zi-Li et al. Self-tuning information fusion Kalman predictor weighted by diagonal matrices and its convergence analysis
CN114218875A (zh) 一种用于流场预测的加速方法及装置
CN103365965B (zh) 一种数据的汇总处理方法和装置
Li et al. An automated operational modal analysis algorithm and its application to concrete dams
CN108984976B (zh) 一种基于加速度响应结构灵敏度计算方法
CN104580660B (zh) 一种移动智能终端及其计步方法、系统
CN108493936A (zh) 基于子空间辨识方法的电力系统低频振荡估计的改进方法
CN113221475A (zh) 一种用于高精度流场分析的网格自适应方法
CN107727110A (zh) 一种步数的统计方法及装置
CN109299512B (zh) 一种基于质量影响的快速灵敏度分析方法
CN109684723A (zh) 一种二维结构内部声学性能分析方法
CN109299513B (zh) 一种模态频率对质量的灵敏度分析方法
CN102043887A (zh) 基于误差估计的网格自适应方法
CN105973266A (zh) 一种应用于移动终端的节能计步方法及装置
CN109635452A (zh) 一种高效的多峰随机不确定性分析方法
CN115601404A (zh) 用于复杂地形不同攻角风场结果内插的神经网络构建方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant