关于newmark法问题的请教！

spring_zhao

subroutine newmark(count,WMM,WDM,WKM,WEF,U,DiffU,DiffDiffU)
  use constant
  use crank
  use imsl
  implicit none
  integer i,j !下标
  real(8) WMM(61,61), WDM(61,61), WKM(61,61),WEF(61)                      ！输入
  real(8) U(0:400,61),DiffU(0:400,61),DiffDiffU(0:400,61)                          ！输出位移、速度、加速度
  real(8) tempU(1:11,61),tempDiffU(1:11,61),tempDiffDiffU(1:11,61)          ！中间变量：位移、速度、加速度
  real(8) temp1(61),temp2(61),TempWMM(61,61),TempWDM(61,61)        ！中间变量
  real(8) K(61,61),P(61)                                                                        ！中间变量
  integer count                                                                                    ！代表所求时刻0.1、0.2、0.3**********
  real(8)::a=5,b=0.25,dt=0.01
  real(8) A0,A1,A2,A3,A4,A5,A6,A7,crankrad,temp                                   ！参数
  A0=1/(b*dt**2)
  A1=a/(b*dt)
  A2=1/(b*dt)
  A3=1/(2*b)-1
  A4=a/b-1
  A5=dt*a/(2*b)-dt
  A6=(1-a)*dt
  A7=a*dt
  crankrad=crankspeed*2*pi/60
  temp=crankrad*180/(9*pi)
  WDM=WDM*temp                                                                             ！进行自变量的转换，与方法本身无关
  WMM=WMM*temp**2                                                                       ！进行自变量的转换，与方法本身无关
  tempU(1,:)=U(count-1,:)                                                                    ！设定初始值
  tempDiffU(1,:)=DiffU(count-1,:)
  tempDiffDiffU(1,:)=DiffDiffU(count-1,:)
  TempWMM=A0*WMM
  TempWDM=A1*WDM
  k=TempWMM+TempWDM+WKM
  do i=2,11
    temp1=A0*tempU(i-1,:)+A2*tempDiffU(i-1,:)+A3*tempDiffDiffU(i-1,:)
    temp2=A1*tempU(i-1,:)+A4*tempDiffU(i-1,:)+A5*tempDiffDiffU(i-1,:)
    temp1=WMM.x.temp1
    temp2=WDM.x.temp2
    p=WEF+temp1+temp2
    tempU(i,:)=k.ix.p
    tempDiffDiffU(i,:)=A0*(tempU(i,:)-tempU(i-1,:))-A2*tempDiffU(i-1,:)-A3*tempDiffDiffU(i-1,:)
    tempDiffU(i,:)=tempU(i-1,:)+A6*tempDiffDiffU(i-1,:)+A7*tempDiffDiffU(i,:)
  enddo
  U(count,:)=tempU(11,:)
  DiffU(count,:)=tempDiffU(11,:)
  DiffDiffU(count,:)=tempDiffDiffU(11,:)  
  do i=1,11
    write(8,"(61f60.8)")tempU(i,:)
  enddo
  return
end subroutine newmark

spring_zhao

我这个是在fortran中编的。
问题是：不收敛
说明：
（1）我无论怎么调我的时间步长，都没有用；反而是我时间步长越小，发散的越快
（2）我无论怎么改变参数a、b都不起效果
（3）我初始条件在第一次调用子函数时是自己任意给的，但也试了零和非零的没有效果。后面再调用时就是前一次计算的结果
最后恳请前辈们能帮我分析一下，到底是哪里出了问题？
万分感激！

wanyeqing2003

看看建立的力学模型是否有问题。

[ 本帖最后由 无水1324 于 2007-11-27 13:54 编辑 ]

spring_zhao

刚度矩阵、质量矩阵、阻尼矩阵我都检查了，没检查出毛病。
关键是程序中的temp1，temp2变化很大，刚开始只有0.000001这样子，后来增加到10000000数量级，后面还有更大的，不可思议级了。
我这个程序编的对吗？

这个程序的含义是这样的：
      我要对某个机构进行机构弹性动力学分析，现在呢就是想求出机构运转一周各个杆件的弹性变形过程。
      我把一个周期分成40个等分，在某一个等分里假设构件的刚度矩阵、质量矩阵、阻尼矩阵都是不变的；再对每一等分应用newmark方法求解，初始值即为前一等分的最后的值，这样是否正确。
      还有就是我直接对整个周期应用newmark方法，请问这两种应用有什么区别吗，我最后算出来的结果是不一样的，虽然都是发散的。

[ 本帖最后由 spring_zhao 于 2007-11-27 14:01 编辑 ]

wanyeqing2003

建议先用一个简单的例子调试，待方法本身通过以后，再带入复杂的模型运算。这样便于检查。

spring_zhao

哦，这个建议不错，我这就去试一试

spring_zhao

Program newmark
  implicit none
  integer i
  real(8) j !下标
  real(8) WMM, WDM, WKM,WEF
  real(8) U(0:30),DiffU(0:30),DiffDiffU(0:30)
  real(8) K,P
  real(8)::a=5,b=0.25,dt=0.1
  real(8) A0,A1,A2,A3,A4,A5,A6,A7
  A0=1/(b*dt**2)
  A1=a/(b*dt)
  A2=1/(b*dt)
  A3=1/(2*b)-1
  A4=a/b-1
  A5=dt*a/(2*b)-dt
  A6=(1-a)*dt
  A7=a*dt
  U(0)=0
  DiffU(0)=0
  DiffDiffU(0)=0
  WMM=0.02
  WDM=0.028
  do i=1,30
    j=i*0.10000000000000000
if(i>=0.and.i<=10) then
   WKM=1
else
   WKM=0.1
endif
WEF=5*(1-j/2.0)**2
    k=A0*WMM+A1*WDM+WKM
p=WEF+WMM*(A0*U(i-1)+A2*DiffU(i-1)+A3*DiffDiffU(i-1))+&
       WDM*(A1*U(i-1)+A4*DiffU(i-1)+A5*DiffDiffU(i-1))
U(i)=p/k
DiffDiffU(i)=A0*(U(i)-U(i-1))-A2*DiffU(i-1)-A3*DiffDiffU(i-1)
DiffU(i)=U(i-1)+A6*DiffDiffU(i-1)+A7*DiffDiffU(i)
write(*,"('time=',3f12.5)")j,WEF,U(i)
  enddo
end Program newmark

这是我的简单例子，但是算出来的结果和正确的结果相差很大，应该是不正确的，请问增量的和我给出的这种newmaik法有什么区别吗，我倒是觉得没什么区别

spring_zhao

DiffDiffU(i)=A0*(U(i)-U(i-1))-A2*DiffU(i-1)-A3*DiffDiffU(i-1)
DiffU(i)=U(i-1)+A6*DiffDiffU(i-1)+A7*DiffDiffU(i)
我发现数值算法，这两个应变量的不同写法，会带来不同的结果

spring_zhao

可能是我的算法出了问题，我发现同样的非线性系统，用全量的newmark法计算是发散的，而且怎么调b都没有用，而用增量的newmark法计算时就可以收敛。
到底是什么机理导致他们这样，还有待于研究。

[ 本帖最后由 spring_zhao 于 2007-11-28 11:19 编辑 ]

spring_zhao

终于找到问题了，全量remark法，也是收敛的。
但是有一个问题，全量和增量有什么区别呢，请问前辈们，他们的区别表现在哪呢？

spring_zhao

原来还是算法的问题
谢谢wanyeqing2003

spring_zhao

看看这个就是我得出的结果，表示杆件上某点的振动，不知道是不是收敛了，它的含义说明什么？

spring_zhao

请高手过来看看阿，都没有人和我讨论问题？:'(

weideyong8

回复 10 # spring_zhao 的帖子

你好，我与你遇到了同样的问题， 能说说你是怎么解决的问题吗？是算法的问题吗  还是编程出的问题？谢谢回复！！！！

mayuanzhuo

楼主你好，我在用newmark法解非线性微分方程组，里面恢复力项有一个非线性回滞环，用newmark法时估计需要用到迭代，求楼主指教！