24小时热门版块排行榜

返回列表

gailisen

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 93.6
帖子: 11
在线: 10.3小时
虫号: 2483886
注册: 2013-05-27
专业: 机械动力学

[求助] 采用runge-kutta法求解常微分方程组时，时间步长怎么设定，积分时间最长为多少

我编制了一个验证弹性碰撞碰撞开始时和碰撞结束时的速度等值反向的FORTRAN程序。
采用的runge-kutta算法，使用的NETLIB上的程序包rksuite。我设置得误差限是1e-6，接触刚度为1E10数量级的，我设置得步长是1e-6s，积分总时间为0~0.1s。
我通过计算发现弹性碰撞开始到弹性碰撞结束发生在1E-6s的数量级上。我把这个过程找到了。并验证了程序的正确性。

但是，我们老师给我提出了一个问题，如果想要用此程序计算第100s，第1000s的结果，我的程序能不能实现。
我通过调试发现：最多能够计算到0.2s左右。我试过了增大、减小积分步长，增大、减小误差限，但是效果都不明显，都只能计算到0.2s左右。
而且，我发现一般参考书上的解常微分方程的积分时间最多也只是1s左右。

我想问问各位前辈，该怎样修改程序才能够计算到任意时间，比如1000s的结果？

希望各位前辈不吝赐教。

我采用的runge-kutta求解器的下载地址为：
http://www.netlib.org/ode/rksuite/rksuite.f
下面是我编写的微分方程的主程序部分：
PROGRAM MAIN
   IMPLICIT NONE
   EXTERNAL F, SETUP, STAT, UT
   INTEGER, PARAMETER :: NEQ = 4, LENWRK = 128, Nr = 1, DOF = 2,
   1                      nDOF = 4, METHOD = 2
   INTEGER(8) K
   INTEGER I, NOUT, STPCST, STPSOK, TOTF, UFLAG
   LOGICAL ERRASS, MESAGE
   DOUBLE PRECISION T,Y(NEQ),YP(NEQ), TSTART, YSTART(NEQ), TEND, TOL,
   1                THRES(NEQ), HSTART, WORK(LENWRK), TWANT,
   2                HNEXT, TINC, TLAST, WASTE, YMAX(NEQ),
   3                Rr, Mr, Jr, L, Ke, Ce, de, mu, Ri, Mi, Mt0, w,
   4                Fx, Fy, Ji, Rou, PI, MCHPES, DWARF

   INTRINSIC COS, SIN

! 各变量赋初始值
   TSTART = 0.0
   YSTART = (/0.095,0.0,0.2,0.0/)
   TLAST = 10
   TEND = 15
! 圆周率
   PI = 3.1415926
! 滚动体参数
! 滚动体端面圆半径(m)
   Rr = 0.0135
! 滚动体质量(kg)
   Mr = 0.213
! 滚动体转动惯量
   Jr = 0.5*Mr*Rr*Rr
! 接触参数
! Hertz线接触系数(N/m^2)
   K = 8.06D10
! 滚动体与滚道接触长度(m)
   L = 0.048
! 接触刚度
   Ke = K*L**(8/9)

! 接触摩擦系数
   mu = 0.0
! 内圈参数
! 内圈滚道半径(m)
   Ri = 0.08
! 外圈参数
! 外圈滚道半径(m)
   Rou = 0.11

! 初始化输出结果
   OPEN(8,FILE='RESULT_3.TXT')
   WRITE(8,'(A,I10)') 'Template 1a with METHOD =', METHOD
   WRITE (8,'(/A/)') ' t          y '
   WRITE (8,'(1X,F11.6,4(3X,F11.6))') TSTART, YSTART

! 设置误差限
   TOL = 1D-4
   THRES = 1D-12

! 调用子程序SETUP
   MESAGE = .TRUE.
   ERRASS = .FALSE.
   HSTART = 0.0D0
   CALL SETUP(NEQ,TSTART,YSTART,TEND,TOL,THRES,METHOD,'Usual
   1          Task',ERRASS,HSTART,WORK,LENWRK,MESAGE)

! 计算NOUT等分点处的解
   NOUT = 1D8
   TINC = (TLAST-TSTART)/NOUT

   DO 40 I =1, NOUT
      TWANT = TLAST + (I-NOUT)*TINC
      CALL UT(F,TWANT,T,Y,YP,YMAX,WORK,UFLAG)

      IF (UFLAG.GT.2) GO TO 60

! Success. T = TWANT. Output computed and true solution components.
   WRITE (8,'(1X,F11.6,4(3X,F11.6))') T, Y
40 CONTINUE
! The integration is complete or has failed in a way reported in a
! message to the standard output channel.
60 CONTINUE
! YMAX(L) is the largest magnitude computed for the solution component
! Y(L) in the course of the integration from TSTART to the last T.  It
! is used to decide whether THRES(L) is reasonable and to select a new
! THRES(L) if it is not.

   WRITE (8,'(A/)') '          YMAX(I) '
   DO 80 I = 1, NEQ
      WRITE (8,'(13X,1PE8.2)') YMAX(I)
80 CONTINUE

! 调用子程序STAT
   CALL STAT(TOTF,STPCST,WASTE,STPSOK,HNEXT)
   WRITE (8,'(/A,I10)')
   1    ' The cost of the integration in evaluations of F is', TOTF

   STOP
   CLOSE(8)
   END PROGRAM MAIN

! *******************************************************************************

!微分方程组

   SUBROUTINE F(T,Y,YP)
   IMPLICIT NONE
   INTEGER(8) K
   DOUBLE PRECISION T, Y(4), YP(4), nh(3), Fni(3), v(3), vt(3),
   1    Ffi(3), nt(3), Fno(3), Ffo(3), PI, Rr, Mr, Jr, L, Ke, mu,Ri,
   2    Rou, dr, delta_i, delta_o, Ro, Fn_i, Ff_i, d_vt, Fn_o, Ff_o


! 圆周率
   PI = 3.1415926
! 滚动体参数
! 滚动体端面圆半径(m)
   Rr = 0.0135
! 滚动体质量(kg)
   Mr = 0.213
! 滚动体转动惯量
   Jr = 0.5*Mr*Rr*Rr
! 接触参数
! Hertz线接触系数(N/m^2)
   K = 8.06E10
! 滚动体与滚道接触长度(m)
   L = 0.048
! 接触刚度
   Ke = K*L**(8/9)

! 接触摩擦系数
   mu = 0.0
! 内圈参数
! 内圈滚道半径(m)
   Ri = 0.08
! 外圈参数
! 外圈滚道半径(m)
   Rou = 0.11

! 滚动体几何中心到坐标原点距离
   dr = Y(1)**2 + Y(2)**2; dr = sqrt(dr)
! 接触法线方向
   nh =(/Y(1),Y(2),0/) ; nh = nh/dr

! 滚动与内圈滚道之间嵌入量
   delta_i  = 0
! 滚动与外圈滚道之间嵌入量
   delta_o = 0
   IF ((dr-Rr).LT.Ri) THEN
      delta_i = Ri - (dr-Rr)
   ENDIF
   IF ( dr+Rr .GT. Rou ) THEN
      delta_o = (dr+Rr) - Rou
   ENDIF
! Fn: 法向接触力，Ff: 切向摩擦力
   Fn_i = 0 ; Ff_i = 0
! 滚动体与内圈滚道之间作用力大小
   IF (  delta_i .GT. 0 ) THEN
      Fn_i = K*delta_i ; Ff_i = mu*Fn_i
   ENDIF
! 滚动体与内圈滚道之间作用力（定义在全局惯性坐标系中）
   Fni = Fn_i*nh
! 滚动体相对于接触面（内圈）的相对运动速度
   v  = (/Y(3),Y(4),0/)
   vt =  v - nh*dot_product(v,nh)
   d_vt = sqrt(vt(1)**2 + vt(2)**2 + vt(3)**2)
   IF  ( d_vt .LT. 1e-12 ) THEN
      Ffi = (/0,0,0/)
   ELSE
      nt = -vt/d_vt
      Ffi = Ff_i *nt
   ENDIF

! 滚动体与外圈滚道之间作用力大小
   Fn_o = 0 ; Ff_o = 0
   IF (  delta_o .GT. 0 ) THEN
      Fn_o = K*delta_o
      Ff_o = mu*Fn_o
   ENDIF
! 滚动体与外圈滚道之间作用力（定义在全局惯性坐标系中）
   Fno = Fn_o*(-nh)
! 滚动体相对于接触面（外圈）的相对运动速度
   v  = (/Y(3),Y(4),0/)
   vt =  v - nh*dot_product(v,-nh)
   d_vt = sqrt(vt(1)**2 + vt(2)**2 + vt(3)**2)
   IF ( d_vt .LT. 1e-12 ) THEN
      Ffo = (/0,0,0/)
   ELSE
      nt = -vt/d_vt
      Ffo = Ff_o *nt
   ENDIF
! 运动微分方程
   YP(1) = Y(3)
   YP(2) = Y(4)
   YP(3) = ( Fni(1) + Ffi(1) + Fno(1) + Ffo(1) )/Mr
   YP(4) = ( Fni(2) + Ffi(2) + Fno(2) + Ffo(2) )/Mr
   RETURN

   END SUBROUTINE F

回复此楼