滑动轴承油膜力 计算出现cl^2-x^2-y^2小于零 哪位高人指点

chao

各位好！我刚接触转子动力学，建立了关于一个jeffcott转子的动力学方程组，其中非线性油膜力
采用：
px=-mu*pi*R*l^3*[(w*y+2*x)/(2*(cl^2-x^2-y^2))^(3/2)+3*x*(x*dx+y*dy)/(cl^2-x^2-y^2)^(5/2)]
py=-mu*pi*R*l^3*[(2*dy-w*x)/(2*(cl^2-x^2-y^2))^(3/2)+3*y*(x*dx+y*dy)/(cl^2-x^2-y^2)^(5/2)]
式中，x,y 为轴中心位移，cl为轴承间隙，前面的系数mu*pi*R*l可以不管。
现在是我用龙哥库塔法求解时，有时候cl^2-x^2-y^2出现小于零的情况，从而得到复数的位移解。这是不合理的，如果此项小于零，则说明轴接触到轴承，并且跑到轴承外了。
我的x，y初始值取的是1e-7，cl的初始值是0.01，保证cl^2-x^2-y^2初始值是合理的，但是这种情况还会出现。求高人指点。

TNC

参考http://forum.vibunion.com/forum- ... 744-fromuid-20.html

chao

回复 2 # TNC 的帖子

谢谢！我看了那个帖子，上面讲到：
“带有开方的这种模型在数值分析中确实很头疼，一般解决办法主要有三种途径
很容易出现DOMAIN error这类错误
1. 采取足够小的步长，不过计算代价呈几何指数增加
2. 给出足够好的初值，初值问题本身就是个难题
3. 人为判断sqrt中是否出现了负值，然后重新赋值，这种办法仅对部分问题有效”
其中只有第一种方法具有实际操作性，但是非常费时间。第二种和第三种说着容易做着难。
我想除了步长之外，是不是和动力学方程中的参数选取也有很大关系呢? m,k,c什么的。这样的话参数太多了，不知道改固定哪个，改变哪个？我是第一次做这方面的，还请多多指点。

chao

最新经验。此问题已解决。
1、首先保证初始值满足cl^2-x^2-y^2大于零；
2、设置较高的龙格库塔法迭代精度，
     此时，如出现cl^2-x^2-y^2小于零，则输出该项，看其是10的几次方，然后修改龙格库塔迭代精度高于此精度1到2个数量级。
这是我的解决办法，不知道是否适合其它朋友，如果遇到此种情况，不妨一试。多谢TNC的回答。谢谢！

yejet

chao 发表于 2010-12-1 06:55

                               
登录/注册后可看大图

最新经验。此问题已解决。
1、首先保证初始值满足cl^2-x^2-y^2大于零；
2、设置较高的龙格库塔法迭代精度 ...

2、设置较高的龙格库塔法迭代精度，
     此时，如出现cl^2-x^2-y^2小于零，则输出该项，看其是10的几次方，然后修改龙格库塔迭代精度高于此精度1到2个数量级。

这种处理方式代价视乎不会太小，当然对于解决问题应该问题不大
至于3楼说的步长问题，可以采用自适应步长，效果应该不会比上述处理方法差
而对于初值问题，可以考虑在同伦算法上是否能够做一些文章

chao

回复 5 # yejet 的帖子

莫非是高人？多谢回答。那么自适应步长是怎么设置的？可有代码供参考，多谢。还有，初值问题能否详细说说，你现在的说的有些笼统，什么是同伦算法？如何做一些文章？

yejet

chao 发表于 2010-12-20 11:19

                               
登录/注册后可看大图

回复 5 # yejet 的帖子

莫非是高人？多谢回答。那么自适应步长是怎么设置的？可有代码供参考，多谢。还有 ...

自适应步长简单的来说就是根据计算的收敛情况在程序中自动的改变计算步长
假设原来的步长位移1，当计算过程中发现不收敛的情况，则自动将步长减半，即变为1/2
这样即可以兼顾计算效率，有可以保证计算的收敛性

yejet

chao 发表于 2010-12-20 11:19

                               
登录/注册后可看大图

回复 5 # yejet 的帖子

莫非是高人？多谢回答。那么自适应步长是怎么设置的？可有代码供参考，多谢。还有 ...

关于同伦算法你可以找找相关的数据
同伦算法又称为延拓算法，国外现在有现成的通用代码，开源的你可以找找

chao

回复 8 # yejet 的帖子

多谢回答。我再看看。