气动弹性时域分析正在进行中

FAI_gavin

问题描述：利用时域分析方法分析二元翼型的颤振特性，求出当量颤振速度file:///C:/DOCUME~1/gavin/LOCALS~1/Temp/msohtml1/01/clip_image002.gif和当量颤振频率file:///C:/DOCUME~1/gavin/LOCALS~1/Temp/msohtml1/01/clip_image004.gif。
其实就是通过无量纲化之后求解一个二阶微分方程组而已，只不过是带参变量的微分方程组。

即：file:///C:/DOCUME~1/gavin/LOCALS~1/Temp/msohtml1/01/clip_image002.gif






其中系数只有file:///C:/DOCUME~1/gavin/LOCALS~1/Temp/msohtml1/01/clip_image002.gif作为参变量，也是最终要求的结果。如何求解呢？通过判断此方程关于某个file:///C:/DOCUME~1/gavin/LOCALS~1/Temp/msohtml1/01/clip_image002.gif是否收敛。


具体思路：通过ode45来求解。




遇到的麻烦：无论我如何给定file:///C:/DOCUME~1/gavin/LOCALS~1/Temp/msohtml1/01/clip_image002.gif，其都是收敛的。（与 提示不相符合）


正在解决中。。。

FAI_gavin

气动弹性大作业.doc (97.5 KB, 下载次数: 38)

2008-6-21 01:15 上传

点击文件名下载附件

还是给个附件吧

FAI_gavin

因为是初学matlab，进展比较小，先公布如下，不过问题依旧是无法确定关于参变量Vf的发散工况。现把程序写于下，望多加指正：


现把转化为一阶：
function dz=flutter(t,z)
%定义控制方程系数
syms h12 alpha12 h11 alpha11 h10 alpha10 h22 alpha22 h21 alpha21 h20 alpha20 real
%定义已知量
a=-0.4;xa=0.05;ra=0.25^(1/2);omega=sqrt(0.6);u=20;Vf=1.1;
%用已知量表示方程系数
h12=1;alpha12=xa;h11=2*Vf/(sqrt(u));alpha11=(1-2*a)*Vf/(sqrt(u));h10=omega^2;alpha10=0;
h22=xa;alpha22=ra^2;h21=(-2*a-1)*Vf/(sqrt(u));alpha21=2*a^2*Vf/(sqrt(u));h20=0;alpha20=ra^2;
dz(1)=z(2);
dz(2)=((alpha22*h11-alpha12*h21)*z(2)+(alpha11*alpha22-alpha21*alpha12)*z(4)+(h10*alpha22-h20*alpha12)*z(1)+(alpha10*alpha22-alpha20*alpha12)*z(3))/(h22*alpha12-h12*alpha22);
dz(3)=z(4);
dz(4)=((h11*h22-h12*h21)*z(2)+(alpha11*h22-alpha21*h12)*z(4)+(h10*h22-h12*h20)*z(1)+(alpha10*h22-alpha20*h12)*z(3))/(h12*alpha22-h22*alpha12);
dz=[dz(1);dz(2);dz(3);dz(4)];
再引入4阶龙格库塔
function [k,T,Z,P]=RK4(dzdt,a,b,CT,h)
n=fix((b-a)/h);
T=zeros(n+1,1);
Z=zeros(n+1,length(CT));
T=a:h:b
Z(1,:)=CT;
for k=1:n
    k1=feval(dzdt,T(k),Z(k,:));
    t2=Z(k)+h/2;
    z2=Z(k,:)'+k1*h/2;
    k2=feval(dzdt,t2,z2);
    k3=feval(dzdt,t2,Z(k,:)'+k2*h/2);
    k4=feval(dzdt,T(k)+h,Z(k,:)'+k3*h);
    Z(k+1,:)=Z(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6;
    k=k+1;
end
T=T(1:n+1);
Z(1:n+1,:);
k=1:n+1;
P=[k',T',Z];
最后就是plot大致判断
CT=[0;1;0;1];h=0.005;[k,X,Y,P]=RK4(@flutter,0,20,CT,h); plot(X,Y(:,1),'g-')

sigma665

具体有什么错误
请贴错误提示

FAI_gavin

整个程序不出错，只是结果感觉不合理。
总结了两点：（1）定义太多符号变量好像会影响运算速度
            （2）在微分方程中含有参变量的问题，能不能实现程序的自动化。
现在我还是没找着参变量Vf的临界值，因为结果都是收敛的。
谢谢sigma665的关心。

sogooda

如果龙格库塔法是自己编程的，那只要多加一个传入参数就可以实现自动化了吧。
如果使用ode45的话，请参阅帮助或版面搜索相关讨论。

ChaChing

LZ问题解决没?