- 1、本文档共11页,可阅读全部内容。
- 2、原创力文档(book118)网站文档一经付费(服务费),不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
- 3、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。如您付费,意味着您自己接受本站规则且自行承担风险,本站不退款、不进行额外附加服务;查看《如何避免下载的几个坑》。如果您已付费下载过本站文档,您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“版权申诉”(推荐),也可以打举报电话:400-050-0827(电话支持时间:9:00-18:30)。
查看更多
华中科技大学流固耦联悬臂输液管动态特性微分求积法附源代码研究一端固支一端悬臂输液管道的动态特性。管道的流固质量比为,无量纲流速为,最大值为14,悬臂输液管道的无量纲控制方程为为管道的横向位移,为无量纲的轴向位置,范围为0~1,为时间项图1 一端固支一端悬臂输液管模型示意图1、问题分析与方法选择对于一边固支,一边悬臂形式的输液管道,边界条件为当时挠度(横向位移对轴向位移零阶导数)与转角(横向位移对轴向位移一阶导数)均为0,即当时,曲率(横向位移对轴向位移二阶导数)为0,而且曲率的空间变化率(横向位移对轴向位移三阶导数)为0,即最近 30 年,微分求积法(DQM,Differential Quadrature Method)逐渐成为一种求解边值和初值问题的有效方法,逐渐被人们认可,广泛地应用于工程结构的振动分析之中。该方法的本质是把函数在给定网点处的导数值近似地用域上全部网点处函数值的加权和表示对于输液管道模型,传统的作法是采用Galekrin法将控制方程进行离散化,仅提取有限的低阶振型函数作近似处理;DQM则克服了Galekrin法的这种局限性,本质上考虑了所有振型的综合贡献。而且DQM的实施过程避免了Galekrin的大量积分运算,其计算精度和效率较好。所以本文最终选取DQM方法。2、结构离散、计算权系数矩阵微分求积法的核心是将给定网格节点的导数值由由网格域上所有节点的函数值的加权和近似地表示表示。为此,首先需要对结构进行网格离散,对于一维的输液管结构,可以划分成轴向分布的一系列网点。假定待求函数为,自变量的范围为,网点的划分总个数为N,如下图所示,各网点的为待求函数。需要指出的是,在实用传统的DQM求解超过二阶的微分方程时,一般在边界点需引入微小量(),以便于添加系统的边界条件。边界点之外的网点称为内部网点,可以采取不同的划分方式:均匀划分与非均匀划分。,,,(边界网点预设微小量),(内部网点均匀划分),(内部网点非均匀划分)图2输液管网点划分示意图加权系数的计算方法一般由拉格朗日插值多项式推导得到。在区间内,假设函数的n个节点的函数值利用Lagrange插值近似为:其中为拉格朗日插值多项式则一阶导数为将多项式形式转化为权系数矩阵形式,权系数矩阵有如下规则:通过矩阵连乘得到更高阶的权系数矩阵3、将DQM格式应用于输液管的运动方程和边界条件,并组装矩阵,计算求解at at 可得将上述关系合并,写成矩阵形式为了方便求解方程,一般将边界网点(包括网点1,2,N-1,N,用表示)放在一起,剩余为内部场点(用表示),并相应地对矩阵内部元素的位置进行调整,所以其中为4X4矩阵:其中为4X(N-4)矩阵:由第一行的关系,可以得到等式由此可以得到边界点位移到内部点位移的转移矩阵T将矩阵T带入到针对内部网点的式子,消除边界网点即有对于,,从而得到假设解的形式,代入上面的动力学方程组可得:有非零解的条件是:对于这种广义特征值问题,可以令,从而得到利用MATLAB得到矩阵特征值。由此得到复数形式的运动频率。再据此分析悬臂输液管的运动特性3、分析所得结果,绘制运动响应图形(略)分析前三阶结果,绘制分叉图与Argand图图3悬臂输液管的无量纲频率的实部、虚部随无量纲流速的变化曲线图4 悬臂输液管的Argand图时,一阶模态的实频率变为0,反映的是失稳。(First Mode Divergency)时,输液管的运动不稳定,可能同时存在第二阶或第三阶运动时,输液管的运动不稳定,可能同时存在第一阶或第二阶运动时,一阶模态的实频率变为0,反映的是失稳。时,一阶模态的实频率变为0,可能同时存在第二阶或第三阶运动附:程序源代码clc;clear;%??·?í???%1üμà2?á?beta=0.2;%á÷1ì?êá?±è?a0.2umax=14;uu=0:0.1:umax;%á÷?ù·??§ll=length(uu);omega1=1:ll;omega2=1:ll;omega3=1:ll;for mm=1:ll u=uu(mm);%u=2;%á÷?ù%μè??????·??a30·YN=35;delta=1e-4;%?ù?è??·?% xx=1:N;% for ii=3:N-2% xx(ii)=(ii-2)/(N-3);% end% xx(1)=0;% xx(2)=delta;% xx(N-1)=1-delta;% xx(N)=1;%·??ù?è??·?xx=1:N;for ii=3:N-2 xx(ii)=0.5*(1-cos(((ii-2)/(N-3))*pi));endxx(1)=0;xx(2)=delta;xx(N-1)=1-delta;xx(N)=1;%?óμ??óè¨A???óμ?±àD′AA1=zeros(N,N);AA2=ze
文档评论(0)