Drillstring Vibrations: Classic Ideas and Modern Approaches - Wave Propagation in Drilling, Well Logging and Reservoir Applications - page 164

surface boundary condition in Equation 4.2.91, the rock-bit interaction model in
Equations 4.2.92 and 4.2.93, and the displacement source model in Equations
4.2.96 and 4.2.97. Once the tridiagonal matrix solver is called and executed, the
solution for the axial displacement U i,n is available for post-processing as
discussed. The solution at the bit will be used to determine appropriate torsional
bit boundary conditions, discussed later; the complete axial solution is used to
evaluate variable coefficients in the bending equations also described later. The
reader should experiment with our Fortran. For example, by removing the rock-
bit and displacement source logic, and replacing it with B(1) = 1.0, C(1) = 0.0,
and W(1) = UZERO*SIN(WMEGA*T) , deleting the gravity term in the difference
equation and initializing with zero displacements, a transient implementation of
the Dareing and Livesay (1968) model is obtained.
.
.
C INITIALIZE AXIAL, ANGULAR AND BENDING DISPLACEMENTS
DO 130 I=1,IMAX
X = XS*(I-1)/(IMAX-1)
UNM1(I) = RHO*G*(X**2)/(2.*ELAST) +WZERO*X/AE
UNM2(I) = UNM1(I)
TNM1(I) = 0.
TNM2(I) = 0.
VNM1(I) = 0.
VNM2(I) = 0.
WNM1(I) = 0.
WNM2(I) = 0.
130 CONTINUE
C BEGIN TIMEWISE INTEGRATION
T = 0.
DO 900 N=1,NMAX
T = N*DT
C AXIAL VIBRATIONS
DO 200 I=2,IMAXM1
A(I) = 1.
C(I) = 1.
B(I) = -2.-RHO*DX*DX/(ELAST*DT*DT)-GAMA*DX*DX/(2.*ELAST*DT)
W(I) = -RHO*DX*DX*(2.*UNM1(I)-UNM2(I))/(ELAST*DT*DT)
1 -GAMA*DX*DX*UNM2(I)/(2.*ELAST*DT)+RHO*G*(DX**2)/ELAST
200 CONTINUE
A(IMAX) = -AE/DX
B(IMAX) = TBMASS/(DT**2) +BETA/(2.*DT) +AE/DX +SPRING
C(IMAX) = 99.
W(IMAX) = 2.*TBMASS*UNM1(IMAX)/(DT**2)
1 +(BETA/(2.*DT) -TBMASS/(DT**2))*UNM2(IMAX)-TBMASS*G
A(1) = 99.
VEL = (UNM1(1)-UNM2(1))/DT
IF(VEL.LE.0.) GO TO 205
B(1) = -1.
C(1) = +1.
W(1) = 0.
GO TO 215
205 B(1) = -RBALPH/DX +RBBETA/DT +RBLAMB