一种基于自相关函数的结构风振响应高效频域估计方法与流程

本发明涉及基于自相关函数的风振响应高效频域估计方法。
背景技术：
：风振响应分析是大型复杂结构(例如大跨度屋盖结构、超高层建筑等)抗风设计的重要环节之一。大型复杂结构的振动特性难以用简单的振型表示，要考虑高阶振型的贡献以及振型之间的耦合，因此，脉动风荷载作用下的结构动力响应极为复杂。传统的时域风振响应分析，需要输入大量风荷载时程信息，采用Newmark-β积分法，计算量较大，计算时所需存储空间较大，导致计算效率低。而传统的频域风振响应分析，同样需要输入大量信息，包括风荷载的协方差矩阵，采用数值积分法，计算耗时较长，导致计算效率低。此外，传统频域计算方法难以考虑振型耦合效应，计算误差大。技术实现要素：本发明的目的是为了解决现有风振响应分析方法计算效率低以及计算误差大的缺点，而提出一种基于自相关函数的结构风振响应高效频域估计方法。一种基于自相关函数的结构风振响应高效频域估计方法具体过程为：步骤一：基于风洞试验测得的结构表面脉动风压时程数据pi(t)，计算各风荷载加载点的风压时程数据的平均值风压时程数据的标准差σpi、风压时程数据的自相关函数Rpi(τ)和风压时程数据的相干函数Cohpij(ω)；其中，i、j为风荷载加载点总数Q中任意两个点，1≤i≤Q，1≤j≤Q，Q为风荷载加载点总数，为正整数；步骤二：采用指数函数拟合自相关函数，得到i点的风荷载加载点的无量纲风压谱曲线的峰值频率ωmi；由自功率谱的定义可知Spi(ω)σpi2=1π∫-∞∞Rpi(τ)exp(-iωτ)dτ=2π·ωmiωmi2+ω2]]>式中，Spi(ω)为各加载点的风压自功率谱密度函数；ω为频率；τ为时差；则无量纲风压谱表示为可得出，式中，Si(ω)为无量纲风压谱；步骤三：采用指数函数拟合相干函数，得到该结构的相干指数kc，计算i、j两点的风荷载相干系数；步骤四：提取结构的质量矩阵[M]、刚度矩阵[K]、阻尼矩阵[C]、风荷载加载点与结构自由度的转换矩阵[R]，响应的影响面矩阵[I]；对结构进行模态分析，求解特征方程得到结构的各阶自振频率ωnk，k＝1,2,…,N，k为自振频率所对应的模态阶数，N为结构自由度数，为正整数；取结构振型矩阵其中为第k、r阶振型，k,r＝1,2,…W；W为计算所选取的模态阶数，W为正整数，10≤W≤N；k、r为自振频率ωnk、ωnr所对应的模态阶数；使得计算结构k阶振型的广义阻尼比式中，T为转置矩阵；步骤五：根据步骤一、步骤二、步骤三和步骤四计算解析积分；得到[Σkr]；步骤六：计算模态响应协方差形成模态响应协方差矩阵；步骤七：计算位移响应协方差矩阵[Σx]＝[Ψ]T[Σy][Ψ]，进而可求得其任意响应的协方差矩阵。本发明的有益效果为：本方法设计基于自相关函数的大型复杂结构风振响应高效频域估计方法，解决传统风振响应分析方法在大型复杂结构分析中计算效率低以及计算精度较差的缺点。该方法首先将风洞试验得到的脉动风荷载处理为供计算用的统计量和频谱参数信息，大大简化了风荷载的输入。然后，结合振型分解法，对每个振型，在每个自由度上对风荷载频谱进行简化的解析积分，得到模态响应协方差。最后，通过对耦合振型的组合得到风振响应的协方差。该方法在保证模态耦合计算精度的同时大大提高了计算效率。本发明在保证计算精度的前提下，大大提升了大型复杂结构风振计算的效率，降低了计算所占用的存储空间。充分利用了解析积分的计算原理，将数值积分转化为解析积分的代数运算，避免了数值积分带来的误差，同时也大大提高了计算效率。在实施过程中，步骤明确，可操作性较强，通过实例分析发现，本发明较传统算法效率大为提升，较为实用。采用该发明对某扇形平面大跨度悬挑网架结构进行风振响应分析，与传统分析方法相比，最大误差不超过1％，远优于传统的简化频域算法误差(16.7％)，效率较传统算法提高200倍。附图说明图1为本发明的流程图；图2a为取-6.42，σpi取2.48的结构表面某点风荷载时程按步骤一统计建模的结果示意图，纵坐标为为归一化的风荷载时程，横坐标为t，为时间(秒)，pi(t)为风荷载时程，pi(t)为平均风荷载，σpi为风荷载标准差；图2b为Rpi(τ)＝exp(-ωmi|τ|)，ωmi＝2.40rad/s的结构表面某点风荷载时程按步骤一统计建模的结果示意图，纵坐标为Rpi(τ)，Rpi(τ)为风荷载的自相关函数，ωmi为风荷载频率，横坐标为τ，τ为时差(秒)；图2c为的结构表面某点风荷载时程按步骤一统计建模的结果示意图，纵坐标为为无量纲风荷载功率谱，横坐标为ω，ω为频率(弧度/秒)，Spi(ω)为风荷载功率谱密度函数；图3a为扇形平面悬挑网架结构的有限元模型示意图；图3b为扇形平面悬挑网架结构按步骤二进行模态分析给出的一阶振型(ωn1＝7.4rad/s)结果示意图；图3c为扇形平面悬挑网架结构按步骤二进行模态分析给出的二阶振型(ωn2＝7.5rad/s)结果示意图；图3d为扇形平面悬挑网架结构按步骤二进行模态分析给出的二阶振型(ωn3＝9.9rad/s)结果示意图；图4a为步骤三计算的模态响应均方根结果示意图，RMSmodaldisp为模态响应均方根值，modalorder为模态阶数，frequency为结构自振频率，条形图为模态响应标准差，圆圈为结构各阶模态的频率；图4b为步骤四计算的模态响应协方差矩阵结果示意图；图5是本发明的算法与传统算法在精度和效率上的比较。具体实施方式具体实施方式一：结合图1说明本实施方式，本实施方式的一种基于自相关函数的结构风振响应高效频域估计方法具体过程为：步骤一：基于风洞试验测得的结构(例如大跨度屋盖结构(如两端支承钢结构屋盖跨度大于36m，悬挑长度大于10m的结构，例如多数普通体育场、机场候机楼、会展中心展览馆等的屋盖结构)、超高层建筑(如高度超过100m的高层结构)，)表面脉动风压时程数据pi(t)，计算各风荷载加载点的风压时程数据的平均值风压时程数据的标准差σpi、风压时程数据的自相关函数Rpi(τ)和风压时程数据的相干函数Cohpij(ω)；其中，i、j为风荷载加载点总数Q中任意两个点，1≤i≤Q，1≤j≤Q，Q为风荷载加载点总数，为正整数；步骤二：采用指数函数拟合自相关函数，得到i点的风荷载加载点的无量纲风压谱曲线的峰值频率ωmi；由自功率谱的定义可知Spi(ω)σpi2=1π∫-∞∞Rpi(τ)exp(-iωτ)dτ=2π·ωmiωmi2+ω2]]>式中，Spi(ω)为各加载点的风压自功率谱密度函数；ω为频率；τ为时差；则无量纲风压谱表示为可得出，式中，Si(ω)为无量纲风压谱；步骤三：采用指数函数拟合相干函数，得到该结构的相干指数kc，计算i、j两点的风荷载相干系数；步骤四：提取结构的质量矩阵[M]、刚度矩阵[K]、阻尼矩阵[C]、风荷载加载点与结构自由度的转换矩阵[R]，响应的影响面矩阵[I]；对结构进行模态分析，求解特征方程得到结构的各阶自振频率ωnk，k＝1,2,…,N，k为自振频率所对应的模态阶数，N为结构自由度数，为正整数；取结构振型矩阵其中为第k、r阶振型，k,r＝1,2,…W；W为计算所选取的模态阶数，W为正整数，根据工程经验确定，通常，10≤W≤N；k、r为自振频率ωnk、ωnr所对应的模态阶数；使得计算结构k阶振型的广义阻尼比式中，T为转置矩阵；步骤五：根据步骤一、步骤二、步骤三和步骤四计算解析积分；得到[Σkr]；步骤六：计算模态响应协方差形成模态响应协方差矩阵；步骤七：计算位移响应协方差矩阵[Σx]＝[Ψ]T[Σy][Ψ]，进而可求得其任意响应的协方差矩阵。具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中基于风洞试验测得的结构表面脉动风压时程数据pi(t)，计算各风荷载加载点的风压时程数据的平均值风压时程数据的标准差σpi、风压时程数据的自相关函数Rpi(τ)和风压时程数据的相干函数Cohpij(ω)；具体过程为：由MATLAB中的mean函数求解；σpi由MATLAB中的std函数求解；Rpi(τ)由MATLAB中的xcorr函数求解；Cohpij(ω)由MATLAB中的mscohere函数求解。其它步骤及参数与具体实施方式一相同。具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中采用指数函数拟合自相关函数，得到i点的风荷载加载点的无量纲风压谱曲线的峰值频率ωmi；具体过程为：采用指数函数拟合自相关函数Rpi(τ)＝exp(-ωmi|τ|)，得到i点的风荷载加载点的无量纲风压谱曲线的峰值频率ωmi。其它步骤及参数与具体实施方式一或二相同。具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三中采用指数函数拟合相干函数，得到该结构的相干指数kc，计算i、j两点的风荷载相干系数；具体过程为：采用指数函数拟合相干函数得到该结构的相干指数kc，计算i、j两点的风荷载相干系数式中，Dij表示两风荷载加载点i、j间距离，U表示参考风速，cij为i、j两点的风荷载相干系数，为非负实数。其它步骤及参数与具体实施方式一至三之一相同。具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤五中根据步骤一、步骤二、步骤三和步骤四计算解析积分；得到[Σkr]；具体过程为：根据步骤一、步骤二、步骤三和步骤四计算解析积分；σijkr=2πσpiσpjωmiωmj∫0∞[(ω2-ωnk2)(ω2-ωnr2)-4ξnkξnrωnkωnrω2]·(ω2+ωmiωmj)[(ω2-ωnk2)2+4ξnk2ωnk2ω2]·[(ω2-ωnr2)2+4ξnr2ωnr2ω2]·(ωmi2+ω2)·(ωmj2+ω2)·[1+(cijω)2]dω=2πσpiσpjωmiωmj∫0∞Θ(ω)Λ(--1ω)·Λ(-1ω)dω=σpiσpjωmiωmj·I1I0]]>式中，σijkr为k、r阶耦合模态在i、j两点间的耦合风振响应均方根；σpi为i点的风荷载标准差，σpj为j点的风荷载标准差；ωmi为i点的风荷载加载点的无量纲风压谱曲线的峰值频率；ωmj为j点的风荷载加载点的无量纲风压谱曲线的峰值频率；ωnk为结构k阶自振频率；ωnr为结构r阶自振频率；ξnk为结构k阶阻尼比；ξnr为为结构r阶阻尼比；cij为i、j两点的风荷载相干系数；Θ(ω)为被积函数的分子多项式；θq为分子多项式的系数；ω2q为频率的2q次幂，q＝0,1,2,3；为被积函数的分母复多项式；为复频率；λs为分母多项式的系数；为复频率的s次幂，s＝1～7；；I1为积分的分子行列式；I0为积分的分母行列式；形成W2个Q阶方阵，表示为式中，i,j＝1,2,…Q；k,r＝1,2,…W；W为计算所选取的模态阶数，W为正整数，根据工程经验确定，通常，10≤W≤N；σ1Qkr为σijkr中i＝1,j＝Q；σ2Qkr为σijkr中i＝2,j＝Q；以此类推，σQQkr为σijkr中a＝Q,b＝Q；[Σkr]为k、r阶耦合模态的中间变量矩阵；先计算k＝r的情况，得到[Σkr]，计算其中，k＝1,2,…,W；将从大到小排序为其中，km＝1,2,…,W，m＝1,2,…,W为排序的序号；式中，为第k阶模态响应的方差；[Σkk]为第k阶模态的中间变量矩阵，即[Σkr]中k＝r的情况；针对前Z项，使得计算k<r的[Σkr]，根据对称原则，[Σrk]＝[Σkr]，得到k>r的[Σkr]；所述，1≤Z≤W；1≤n≤W；为第km阶模态响应的方差，即取k＝km的结果。其它步骤及参数与具体实施方式一至四之一相同。具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述分子多项式的系数θq和分母多项式的系数λs具体为：对关于频率ω的分子多项式Θ(ω)针对频率ω合并多项式，得到ω2q的系数θq，q＝0,1,2,3；对关于复频率的分母复多项式针对复频率合并多项式，得到的系数λs，s＝1～7。其它步骤及参数与具体实施方式一至五之一相同。具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述I1和I0具体公式为：当cij≠0时，I1=000θ3θ2θ1θ0-cijλ5-λ3λ10000-λ6λ4-λ2λ0000cij-λ5λ3-λ10000λ6-λ4λ2-λ0000-cijλ5-λ3λ10000-λ6λ4-λ2λ0,]]>I0=cij·λ6-λ4λ2-λ0000-cijλ5-λ3λ10000-λ6λ4-λ2λ0000cij-λ5λ3λ10000λ6-λ4λ2-λ0000-cijλ5-λ3λ10000-λ6λ4-λ2λ0;]]>当cij＝0时，I1=00θ3θ2θ1θ0-1λ4-λ20000-λ5λ3-λ10001-λ4λ2-λ0000λ5-λ3λ1000-1λ4-λ2λ0]]>I0=λ5-λ3λ1000-1λ4-λ20000-λ5λ3-λ10001-λ4λ2-λ0000λ5-λ3λ1000-1λ4-λ2λ0.]]>其它步骤及参数与具体实施方式一至六之一相同。具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述步骤六中计算模态响应协方差形成模态响应协方差矩阵；具体过程为：计算模态响应协方差形成模态响应协方差矩阵其它步骤及参数与具体实施方式一至七之一相同。具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：所述步骤七中计算位移响应协方差矩阵[Σx]＝[Ψ]T[Σy][Ψ]，进而可求得其任意响应的协方差矩阵；具体过程为：计算位移响应协方差矩阵[Σx]＝[Ψ]T[Σy][Ψ]，进而可求得其任意响应的协方差矩阵[Σs]＝[I]T[Σx][I]。其它步骤及参数与具体实施方式一至八之一相同。采用以下实施例验证本发明的有益效果：实施例一：本实施例一种基于自相关函数的风振响应高效频域估计方法具体是按照以下步骤制备的：取一实际工程扇形平面悬挑网架算例进行说明并验证。步骤一～三：将风洞试验测得的风压时程进行统计建模分析，如图2a、图2b、图2c、所示，得到各点的平均、脉动风压、风荷载频率以及相干指数，并利用风压谱检验模型的准确性。步骤四：对有限元模型进行质量、刚度、阻尼等矩阵的提取，并进行模态分析，结果如图3a、图3b、图3c、图3d所示。步骤五：计算积分，形成模态响应标准差，如图4a所示，可以看出前2×2的耦合情况。步骤六：计算模态响应协方差矩阵，结果如图4b所示。步骤七：计算结构位移响应标准差，如图5所示。为验证该发明的有效性，本文还将计算结果与传统方法进行了对比，发现传统风振分析方法计算效率较低，而采用不考虑模态耦合的传统频域计算方法，虽然效率较高，但计算误差较大(达16.7％)。本发明提出的方法在保证计算精度的同时，计算效率大大提升。说明该方法实现了大型复杂结构风振响应的高效计算，具有实用价值。本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。当前第1页1 2 3