Methods and systems for modeling, designing, and conducting drilling operations that consider vibrations

[0001] This application is the National Stage of International Application No. PCT/US09/59040, filed 30 Sep. 2009, which claims the benefit of U.S. Provisional Application Nos. 61/117,015; 61/117,016; 61/117,021; each of which was filed on 21 Nov. 2008, and each of which is incorporated herein by reference in their entirety for all purposes.

[0002] The present disclosure provides methods and systems for modeling, designing, and conducting drilling operations that consider vibrations, which may be experienced by a drilling system. In particular, the present disclosure provides systems and methods for modeling bottom hole assembly (BHA) vibration performance during drilling to enable improved design and operation for enhanced drilling rate of penetration, to reduce downhole equipment failure, to extend current tool durability and/or to enhance overall drilling performance. BHA modeling may be used to enhance hydrocarbon recovery by drilling wells more efficiently.

[0003] This section is intended to introduce various aspects of related technology, which may be associated with exemplary embodiments of the present techniques. This discussion is believed to be helpful in providing information to facilitate a better understanding of particular aspects of the present techniques. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.

[0004] The production of hydrocarbons, such as oil and gas, has been performed for many years. To produce these hydrocarbons, one or more wells are typically drilled into subterranean locations, which are generally referred to as subsurface formations or basins. The wells are formed to provide fluid flow paths from the subterranean locations to the surface. The drilling operations typically include the use of a drilling rig coupled to a drillstring and bottom hole assembly (BHA), which may include a drill bit or other rock cutting devices, drill collars, stabilizers, measurement while drilling (MWD) equipment, rotary steerable systems (RSS), hole opening and hole reaming tools, bi-center bits, roller reamers, shock subs, float subs, bit subs, heavy-weight drill pipe, mud motors, and other components known to those skilled in the art. Once drilling operations are complete, the produced fluids, such as hydrocarbons, are processed and/or transported to delivery locations. As is well understood, drilling operations for the preparation of production wells, injection wells, and other wells are very similar. The present methods and systems may be used in cooperation with providing wells for hydrocarbon production, for injection operations, or for other purposes.

[0005] During the drilling operations, various limiters may hinder the rate of penetration (ROP). For instance, vibrations during drilling operations have been identified as one factor that limits the ROP. These vibrations may include lateral, axial and torsional vibrations, which may be present in a coupled or an uncoupled form. Axial vibrations occur as a result of bit/rock interactions and longitudinal drillstring dynamics; this mode may propagate to surface or may be dampened out by contact with the wellbore. Torsional vibrations may involve fluctuations in the torque at the bit and subsequent propagation uphole as a disturbance in the rotary motion of the drillstring. BHA lateral vibrations involve beam bending dynamics in the stiff pipe near the bit and do not usually propagate directly to the surface. However, lateral vibrations may couple to the axial and torsional vibrations and be experienced at the surface. Some authors have identified lateral vibrations as the most destructive vibrational mode to drilling equipment. The identification of the different types and amplitudes of the vibrations may be provided from downhole sensors in MWD equipment to provide either surface readout of downhole vibrations or stored data that can be downloaded at the surface after the “bitrun” or drilling interval is complete.

[0006] As drilling operations are expensive, processes for optimizing drilling operations based on the removal or reduction of system inefficiencies, or founder limiters, such as vibrations, may be beneficial. The downhole failure of a BHA or BHA component may be expensive and significantly increase the costs of drilling a well. The costs of BHA failures may include replacement equipment and additional time for a round-trip of the drillstring in the event of a washout (e.g., loss of drillstem pressure) with no parting of the drillstring. Further compounding these costs, sections of the wellbore may be damaged, which may result in sidetracks around the damaged sections of the wellbore. While many factors affect the durability of a BHA, vibrations have been identified as a factor that impacts equipment durability.

[0007] Accordingly, design tools (e.g., software applications and modeling programs) may be utilized to examine the drill string and BHA configurations and proposed drilling operations before implementation in a drilling operation. For example, vibrational tendencies may be identified along with drilling conditions, configuration designs, materials, and other operational variables that may affect the vibrational tendencies of the drill string and/or BHA during drilling operations. For example, modeling programs may represent the static force interactions in a BHA as a function of stabilizer placement. Although there have been numerous attempts to model BHA dynamics, there is a need for model-based design tools to simulate BHA designs for evaluating vibration effects as described herein.

[0008] In the numerous references cited in this application, there are both time and frequency-domain models of drilling assemblies. Because of the interest in direct force calculations for bit design and the rapid increase in computational capability, recent activity has focused on the use of direct time domain simulations and the finite element methods, including both two-dimensional and three-dimensional approaches. However, these simulations still require considerable calculation time, and therefore the number of cases that can be practically considered is limited. The finite element method has also been used for frequency-domain models, in which the basic approach is to consider the eigenvalue problem and solve for the critical frequencies and mode shapes. Only a couple of references have used the forced-frequency response approach, and these authors chose different model formulations than those discussed herein, including a different selection of boundary conditions. One reference used a similar condition at the bit in a finite element model, but a different boundary condition was specified at the top of the bottom hole assembly. This reference did not proceed further to develop the design procedures and methods disclosed herein.

[0009] Further, as part of a modeling system developed by ExxonMobil, a vibration performance index was utilized to provide guidance on individual BHA designs. A steady-state forced-frequency response dynamic model was developed to analyze a single BHA in batch mode from a command line interface, using output text files for graphical post-processing using an external software tool, such as Microsoft Excel™. This method was difficult to use, and the limitations of the interface impeded its application. The model has been utilized in some commercial applications within the United States since 1992 to place stabilizers to reduce the predicted vibration levels, both in an overall sense and specifically within designed rotary speed ranges. This model provided an End-Point Curvature index for a single BHA configuration. The End-Point Curvature index was limited to looking at performance from the perspective of a single point at the top of the BHA model. Moreover, the operational limitations of this prior model limited its application to individual BHA configurations for the determination of stabilizer placement. It was not capable of considering multiple BHA configurations conveniently or of conveniently varying a plurality of parameters for optimizing one or more factors other than the stabilizer location.

[0010] Other related material may be found in the following: G. Heisig et al., “Lateral Drillstring Vibrations in Extended-Reach Wells”, SPE 59235, 2000; P. C. Kriesels et al., “Cost Savings through an Integrated Approach to Drillstring Vibration Control”, SPE/IADC 57555, 1999; D. Dashevskiy et al., “Application of Neural Networks for Predictive Control in Drilling Dynamics”, SPE 56442, 1999; A. S. Yigit et al., “Mode Localization May Explain Some of BHA Failures”, SPE 39267, 1997; M. W. Dykstra et al., “Drillstring Component Mass Imbalance: A Major Source of Downhole Vibrations”, SPE 29350, 1996; J. W. Nicholson, “An Integrated Approach to Drilling Dynamics Planning, Identification, and Control”, SPE/IADC 27537, 1994; P. D. Spanos and M. L. Payne, “Advances in Dynamic Bottomhole Assembly Modeling and Dynamic Response Determination”, SPE/IADC 23905, 1992; M. C. Apostal et al., “A Study to Determine the Effect of Damping on Finite-Element-Based, Forced Frequency-Response Models for Bottomhole Assembly Vibration Analysis”, SPE 20458, 1990; F. Clayer et al., “The Effect of Surface and Downhole Boundary Conditions on the Vibration of Drillstrings”, SPE 20447, 1990; D. Dareing, “Drill Collar Length is a Major Factor in Vibration Control”, SPE 11228, 1984; A. A. Besaisow, et al., “Development of a Surface Drillstring Vibration Measurement System”, SPE 14327, 1985; M. L. Payne, “Drilling Bottom-Hole Assembly Dynamics”, Ph.D. Thesis, Rice University, May 1992; A. Besaisow and M. Payne, “A Study of Excitation Mechanisms and Resonances Inducing Bottomhole-Assembly Vibrations”, SPE 15560, 1988; and U.S. Pat. No. 6,785,641.

[0011] The prior art does not provide tools to support a design process as disclosed herein (i.e. a direct characterization of the drilling vibration behavior for myriad combinations of rotary speed and weight on bit), and there are no references to design indices or figures of merit to facilitate comparison of the behaviors of different assembly designs. Accordingly, there is a need for such software tools and design metrics to design improved bottom hole assembly configurations and drilling operations to reduce drilling vibrations.

[0012] The technologies of the present disclosure are directed to methods and systems for representing vibrational performance of drilling equipment. In some implementations, the methods consist of: a) constructing at least one surrogate representing at least a portion of a bottom hole assembly; b) associating at least two virtual sensors with each of the at least one surrogates such that the at least two sensors are spaced longitudinally from each other along each bottom hole assembly; c) utilizing at least one frequency-domain model to calculate at least one state of the at least two virtual sensors during one or more simulated drilling operations for each of the at least one surrogates; d) calculating a transmissibility index between the at least two virtual sensors for each of the at least one surrogates, wherein the transmissibility index is based at least in part on at least one of the calculated states; and e) using the calculated transmissibility index for each of the at least one surrogates to determine the transmissibility of vibrations within the bottom hole assembly.

[0013] Each of the steps outlined above may be carried out with various adjustments and/or specifics within the scope of the present disclosure. For example, the calculated at least one state may comprise at least one of displacement, tilt angle, bending moment, and shear force. One or more of these calculated states may be used to calculate accelerations of the at least two virtual sensors, which, in some implementations, may be used to calculate the transmissibility index, such as by ratio. Similarly, the calculated transmissibility index may be a ratio between any one or more of the calculated states or derivations therefrom.

[0014] Depending on how the transmissibility index is calculated, its numerical value may have different meanings. In implementations where the transmissibility index is calculated as a ratio, a transmissibility index greater than one may predict that vibrations would increase between a first virtual sensor and a second virtual sensor. Similarly, a transmissibility index less than 1 may predict that vibrations would decrease between a first virtual sensor and a second virtual sensor.

[0015] In some implementations, at least one of the virtual sensors may be associated with a bit of the at least one bottom hole assembly surrogate and the transmissibility index may be calculated for a plurality of points along the surrogate. In such implementations, the calculated transmissibility indices may produce a plot wherein peaks of the transmissibility plot indicate locations of local peak vibration in the surrogate bottom hole assembly.

[0016] In some implementations, the methods may further include: f) drilling at least a portion of a well with a bottom hole assembly at least substantially embodying a surrogate used to calculate a transmissibility index while measuring acceleration at least at two sensors disposed along the embodied bottom hole assembly; g) calculating a measured transmissibility index using the measured accelerations; and h) comparing the measured transmissibility index with the transmissibility index of the surrogate. Moreover, some implementations may include updating the at least one surrogate to represent a different bottom hole assembly configuration and repeating steps (b)-(e) from above. Additionally or alternatively, the methods may include modifying drilling operations on the well based at least in part on the measured transmissibility index and the surrogate transmissibility index. Still further, some implementations may include updating one or more of the at least one surrogate, the at least two virtual sensors, the at least one frequency-domain model, and the transmissibility index calculations based at least in part on the comparison of the measured transmissibility index and the transmissibility index of the at least one surrogate.

[0017] Additionally or alternatively, the present disclosure provides methods of drilling a well for use in the production of hydrocarbons. For example, a suitable method may include: a) constructing at least one surrogate representing at least a portion of a bottom hole assembly, wherein the at least one surrogate includes at least two virtual sensors; b) calculating a transmissibility index between the at least two virtual sensors for each of the at least one surrogates; c) selecting an optimized bottom hole assembly configuration for a drilling operation based at least in part on the calculated transmissibility index; and d) drilling a well with drilling equipment incorporating a bottom hole assembly at least substantially embodying the selected bottom hole assembly configuration. In some implementations, the step of drilling the well may be conducted according to a drilling plan developed based at least in part on the calculated transmissibility index. Additionally or alternatively, the step of selecting an optimized bottom hole assembly configuration may comprise selecting different bottom hole assembly configurations for different portions of the drilling operation.

[0018] As with all of the implementations described herein, the methods and systems may be implemented and/or utilized in the production of hydrocarbons. For example, the methods may include the step of producing hydrocarbons from a well drilled with drilling equipment incorporating a bottom hole assembly at least substantially embodying a bottom hole assembly surrogate for which a transmissibility index was calculated.

[0019] Any one or more of the methods described above may include one or more steps adapted to be performed by a computer-based modeling system. Accordingly, the present disclosure is further directed to modeling systems. An exemplary modeling system may include a processor; a memory coupled to the processor; and a set of computer readable instructions accessible by the processor. The set of computer readable instructions may be configured to: a) construct at least one surrogate representing at least a portion of a bottom hole assembly, wherein the at least one surrogate includes at least two virtual sensors; b) calculate a transmissibility index between the at least two virtual sensors for each of the at least one surrogates; and c) output the transmissibility index for use in selecting an optimized bottom hole assembly configuration for a drilling operation based at least in part on the calculated transmissibility index. An exemplary system may include instructions adapted to calculate the transmissibility index utilizing at least one frequency-domain model to calculate at least one state of the at least two virtual sensors during one or more simulated drilling operations for each of the at least one surrogates. Additionally or alternatively, systems within the present disclosure may provide the output as a graphical representation of the transmissibility index of a bottom hole assembly configuration at one or more points along the bottom hole assembly configuration.

[0020] The technologies of the present disclosure are further directed to methods and systems for representing vibrational performance of drilling equipment. In some implementations, the methods consist of: a) constructing at least one surrogate representing at least a portion of a bottom hole assembly disposed in a well; b) utilizing a frequency-domain model to calculate a sideforce at least at one contact point between the bottom hole assembly and the well, wherein the sideforce is calculated as a function of rotational speed for each surrogate; c) determining at least one sideforce slope index as a function of rotational speed for the at least one contact point; and d) displaying the calculated sideforce slope index as a function of rotational speed.

[0021] Each of the steps outline above may be carried out with various adjustments and/or specifics within the scope of the present disclosure. For example, the frequency-domain model used to calculate the sideforce may consider the coefficient of friction to be non-constant over the rotational speeds considered for at least one of the contact forces. Additionally or alternatively, the at least one sideforce slope indices may be determined graphically and/or numerically. In some implementations, the determined sideforce slope index may be a combined index representative of a plurality of contact points between the bottom hole assembly and the well.

[0022] Depending on how the sideforce slope index is calculated, its numerical value may have different meanings. In some implementations, a non-zero sideforce slope index may indicate a greater potential for vibration in that region of the bottom hole assembly. In some implementations, the absolute value of the sideforce slope index may be plotted as function of rotational speed to determine a quantified potential for vibration, which may be used to identify one or more contact points having greater potential for vibration.

[0023] Additionally or alternatively, the present disclosure provides methods of drilling a well for use in the production of hydrocarbons. For example, a suitable method may include: a) constructing at least one surrogate representing at least a portion of a bottom hole assembly disposed in a well; b) determining at least one sideforce slope index as a function of rotational speed for at least one contact point between the bottom hole assembly and the well; c) selecting an optimized bottom hole assembly configuration for a drilling operation based at least in part on the determined at least one sideforce slope index; and d) drilling a well with drilling equipment incorporating a bottom hole assembly at least substantially embodying the selected bottom hole assembly configuration. In some implementations, the step of drilling the well may be conducted according to a drilling plan developed based at least in part on the determined at least one sideforce slope index. Additionally or alternatively, the step of selecting an optimized bottom hole assembly configuration may comprise selecting different bottom hole assembly configurations for different portions of the drilling operation.

[0024] As with all of the implementations described herein, the methods and systems may be implemented and/or utilized in the production of hydrocarbons. For example, the methods may include the step of producing hydrocarbons from a well drilled with drilling equipment incorporating a bottom hole assembly at least substantially embodying a bottom hole assembly surrogate for which a sideforce slope index was calculated.

[0025] Any one or more of the methods described above may include one or more steps adapted to be performed by a computer-based modeling system. Accordingly, the present disclosure is further directed to modeling systems. An exemplary modeling system may include a processor; a memory coupled to the processor; and a set of computer readable instructions accessible by the processor. The set of computer readable instructions may be configured to: An exemplary system may include instructions adapted to: a) construct at least one surrogate representing at least a portion of a bottom hole assembly disposed in a well; b) determining at least one sideforce slope index as a function of rotational speed for at least one contact point between the bottom hole assembly and the well; and d) output the at least one sideforce slope index for use in selecting an optimized bottom hole assembly configuration for a drilling operation based at least in part on the determined at least one sideforce slope index. In some implementations, the step of determining a sideforce slope index utilizes at least one frequency-domain model to calculate a sideforce at least at one contact point. Additionally or alternatively, systems within the present disclosure may provide the output as a graphical representation of the sideforce slope index of a bottom hole assembly configuration at one or more points along the bottom hole assembly configuration.

[0026] Additionally or alternatively, the technologies of the present disclosure are directed to methods of modeling drilling equipment to represent vibrational performance of the drilling equipment. In some implementations, the method includes a) identifying two or more weighted fundamental excitation modes for a drilling bottom hole assembly; b) constructing at least one surrogate representing at least a portion of a bottom hole assembly; c) utilizing a frequency-domain model to simulate a response of the at least one surrogate to excitations corresponding with the identified fundamental excitation modes; d) determining one or more performance indices for the simulated surrogate; and e) utilizing the one or more performance indices in selecting at least one of one or more bottom hole assembly configurations and one or more drilling plans for use in drilling operations. In some implementations, each fundamental excitation mode may be weighted relative to at least one other fundamental excitation mode. Additionally or alternatively, the excitation modes may be related to at least one vibration-related drilling parameter.

[0027] One or more of the determined performance indices may be based at least in part on the simulated response of the surrogate at least at two fundamental excitation modes and on the relative weight of the at least two fundamental excitation modes. The one or more performance indices may be selected from at least one of an end point curvature index, a BHA strain energy index, an average transmitted strain energy index, a transmitted strain energy index, a root-mean-square BHA sideforce index, a root-mean-square BHA torque index, a total BHA sideforce index, a total BHA torque index, a sideforce slope index, a transmissibility index, and any mathematical combination thereof. Other suitable performance indices may be identified.

[0028] In some implementations of the present methods, the methods may further include drilling a well using at least one of a) the selected one or more bottom hole assembly configurations and b) the selected one or more drilling plans.

[0029] The two or more fundamental excitation modes may be identified in a variety of suitable manners. For example, the fundamental excitation modes may be identified from field data using a method including: a) obtaining field-data dynamic measurements of at least one dynamic state of a drilling bottom hole assembly, wherein each of the measurements is associated with at least one node in the bottom hole assembly; processing the field-data measurements to obtain one or more windows having frequency-domain spectra of at least one of the measured dynamic states; and c) identifying two or more fundamental excitation modes in the one or more windows. The fundamental excitation modes may correspond to regions of the frequency-domain spectra having spectral peaks or accumulations. Additionally, each of the two or more fundamental excitation modes is weighted relative to at least one other fundamental excitation mode.

[0030] Continuing with the exemplary field data-based method, the at least one dynamic state may be selected from one or more of rotary speed, displacement, velocity, acceleration, bending strain, bending moment, tilt angle, and force. The field-data may be collected using one or more near-bit sensors. In some implementations, the field-data measurements may be processed using one or more Fourier transforms to provide frequency-domain spectra. Additionally or alternatively, in some implementations, the one or more windows each may present measured data for an interval in a drilling history, wherein the interval is for at least one of a period of time, a depth range, and a rotary speed applied during the drilling. For example, the one or more windows may present intervals of nearly constant rotary speed and the one or more identified fundamental excitation modes may be associated with one or more multiples of the rotary speed having spectral peaks. The field data-based methods may further include drilling a well using at least one of a) the selected one or more bottom hole assembly configurations and b) the selected one or more drilling plans.

[0031] In some implementations, the fundamental excitation modes may be identified from both simulated data and field data. An exemplary method may include: a) obtaining measurements of at least one parameter of a drilling bottom hole assembly indicative of vibrational performance, wherein the measurements relate to one or more nodes on the drilling bottom hole assembly; b) constructing a surrogate representing at least a portion of the drilling bottom hole assembly; c) utilizing a frequency-domain model to simulate a response of the surrogate to dynamic excitations at one or more reference nodes corresponding to the nodes on the drilling bottom hole assembly, wherein a response is simulated for each of at least two excitation modes; d) determining a vibrational performance index for each of the at least two excitation modes based at least in part on the response of the surrogate to the dynamic excitations; e) comparing the at least two determined vibrational performance indices with the obtained measurements to determine the relative contribution of each excitation mode to the measured vibration performance; and f) weighting each of the excitation modes according to the respective relative contributions to determine at least two fundamental excitation modes, which are weighted relative to each other.

[0032] Continuing with the example utilizing both field and simulated data, the at least one measured parameter may be selected from one or more of rate of penetration, mechanical specific energy, measured downhole acceleration, measured downhole velocity, bending moment, bending strain, shock count, and stick-slip vibrations. Such parameters may be collected in any suitable manner using a variety of equipment and methods readily available. In some implementations, the dynamic excitations of the surrogate may be applied by perturbing at least one model state selected from displacement, tilt angle, moment, and force. Additionally or alternatively, in some implementations, the at least two determined vibrational performance indices may be summed with multiplicative non-negative coefficients to obtain a combined surrogate performance index for comparison with the obtained measurements. The surrogate vibrational performance index may be compared with the obtained measurements while varying the non-negative coefficients for each performance index until differences between the combined performance index and the obtained measurements are at least substantially minimized. When those differences are minimized, excitation coefficients are established corresponding to at least two weighted fundamental excitation modes. As with the other methods described herein, the methods utilizing both field and simulated data may further include drilling a well using at least one of a) the selected one or more bottom hole assembly configurations and b) the selected one or more drilling plans.

[0033] The methods described herein may be implemented and/or utilized in the production of hydrocarbons. For example, the methods may include the step of producing hydrocarbons from a well drilled using at least one of a) the selected one or more bottom hole assembly configurations and b) the selected one or more drilling plans.

[0064] In the following detailed description section, the specific embodiments of the present techniques are described in connection with preferred embodiments. However, to the extent that the following description is specific to a particular embodiment or a particular use of the present techniques, this is intended to be for exemplary purposes only and simply provides a concise description of the exemplary embodiments. Moreover, to the extent that a particular feature or aspect of the present systems and methods are described in connection with a particular embodiment or implementation, such features and/or aspects may similarly be included or used in connection with other embodiments or implementations described herein or otherwise within the scope of the invention claimed in this or related applications. Accordingly, the invention is not limited to the specific embodiments described below, but rather, it includes all alternatives, modifications, and equivalents falling within the true scope of the appended claims.

[0065] The present disclosure is directed to methods and systems for modeling, designing, and utilizing bottom hole assemblies to evaluate, analyze, design, and assist in the drilling of wells and in the production of hydrocarbons from subsurface formations. Under the present techniques, a modeling system may include software or modeling programs that characterize the vibration performance of one or more candidate BHA's graphically in what is referred to as “design mode.” In some implementations, the vibration performance of two or more candidate BHA's may be displayed graphically and simultaneously to facilitate comparison of the candidate BHA's. The BHA used in a drilling system may be selected based on one or more relative vibration performance indices for different BHA surrogates. These indices may include point indices, such as an end-point curvature index, and interval indices, such as a BHA strain energy index, an average transmitted strain energy index, a transmitted strain energy index, a root-mean-square (RMS) BHA sideforce index, an RMS BHA torque index, a total BHA sideforce index, a total BHA torque index, a transmissibility index, and a sideforce slope index, which are discussed further below, in addition to specific static design objectives for the respective assembly.

[0066] Further, the present disclosure provides methods and systems that utilize a “log mode” display to compare predicted vibration characteristics with measured data under specific operating conditions. The same indices used in the design mode may be presented in a log mode display to compare measured drilling data with the indices to assist in assessing the BHA vibration performance and to gain an understanding of how to evaluate the different vibration performance metrics by comparison with field performance data (e.g., measured data). For example, and as will be better understood from the description herein, one or more of the data sets from the design mode, including the vibration performance indices, may be compared against measured data and/or data derived from measured data. The comparison may reveal helpful information such as the components of the BHA most likely contributing to the vibrations, the drilling conditions that will avoid vibrations, relative contributions of particular indices, excitation modes, and/or vibrational modes to the actual performance, and other information to aid in improving the modeling process, the BHA design process, and/or the development of drilling operational plans. Additionally or alternatively, this same data may be plotted in a format similar to that used for the vibration performance indices, with rotary speed and/or bit weight on the independent axes, showing the relationships of the measured data to the vibration performance indices. Since this data is normally obtained in a Drilling Vibrations Data Test, this plot is referred to as the “DVDT” display.

[0067] Turning now to the drawings, and referring initially to FIG. 1, an exemplary flow chart 100 of a process of modeling and operating a drilling system in accordance with certain aspects of the present techniques is described. In this process, candidate BHA configurations are represented by surrogates that can be utilized in modeling programs. The modeling programs of the present disclosure provide graphical and/or numerical representations of the how the BHA configuration would operate during implementations under one or more operating conditions. The graphical and/or numerical representations may be presented in the form of one or more indices, which may be evaluated on an absolute or comparative basis to identify a preferred BHA for given operating conditions and/or a preferred set of operating conditions for a given BHA.

[0068] The flow chart begins at block 102. At block 104, data may be obtained for use in the methods of the present disclosure. The data may include well operating parameters (e.g., weight on bit (WOB) range, rotary speed range (e.g., rotations per minute (RPM)), nominal borehole diameter, hole enlargement, hole angle, drilling fluid density, depth, and the like). Some model-related parameters may also be obtained, such as the vibrational excitation modes to be modeled (specified as integer and/or non-integer multiples of the rotary speed and/or specific vibration frequencies), element length, boundary conditions, and number of “end-length” elements and the end-length increment value. Then, one or more BHA surrogates may be constructed, as shown in block 106. The construction of the BHA surrogates includes identifying BHA design parameters (e.g., drill collar dimensions and mechanical properties, stabilizer dimensions and locations in the BHA, drill pipe dimensions, length, and the like). As will be described more thoroughly below, the BHA surrogate may be constructed in a variety of suitable manners provided that the surrogate can be modeled using frequency-domain models.

[0069] In block 108, the operation of the BHA surrogate is modeled using one or more frequency-domain models. The modeling of the BHA surrogates may include consideration of the static solutions and the dynamic solutions. The modeling may include two dimensional models and/or three dimensional models, both of which are described in better detail below. The frequency-domain models provide various data about the operation of the BHA surrogate, which can be used to generate at least one vibration performance index. FIG. 1 illustrates at block 110 the step of determining at least one vibration performance index for a BHA surrogate. Examples of illustrative vibration performance indices are provided below together with examples of possible uses and interpretations of such indices. At least one index is then displayed or otherwise presented to a user or an operator, which is represented by block 112 in FIG. 1. The display or presentation of the vibration performance index may communicate the index to the user in any suitable manner and in any suitable format. For example, the vibration performance index may be presented in numerical and/or graphical formats. Additionally, the index may be presented on a computer display, on a printed page, transmitted to a remote location for presentation, stored for later retrieval, etc. With experience, a BHA design engineer may appreciate the design tradeoffs and, by comparing vibration performance index results for different designs, may develop BHA designs with improved operating performance and/or identify better operating parameters. An example of the design iteration process is described further below.

[0070] FIG. 1 further illustrates that following the determination and display of a vibration performance index, various optional steps may be included in the methods within the scope of the present disclosure. Once modeled, one of the BHA configurations represented by a surrogate may be selected, as shown in block 114. The selection may be based on a comparison of multiple BHA surrogates. That is, the modeling of the BHA surrogates may include different displays of the calculated state vectors (e.g., displacement, tilt, bending moment, lateral shear force of the beam, and BHA/wellbore contact forces and torques) as a function of the operating parameters (e.g., RPM, WOB, etc.), distance to the bit, and BHA configuration. The displayed results or solutions, including the vibration performance indices, may include detailed 3-dimensional state vector plots intended to illustrate the vibrational tendencies of alternative BHA configurations. The selection of a BHA configuration may include selecting a preferred BHA configuration in addition to identifying a preferred operating range for the preferred configuration. The selection may be based on the relative and/or absolute performance of the BHA configurations, which may be evaluated using a variety of indices, including end-point curvature index, BHA strain energy index, average transmitted strain energy index, transmitted strain energy index, RMS BHA sideforce index, RMS BHA torque index, total BHA sideforce index, total BHA torque index, transmissibility index, sideforce slope index, and any mathematical combination thereof. In some implementations, the selection of a BHA configuration may include the selection of a configuration that had been represented by one or more of the BHA surrogates. Additionally or alternatively, the selected BHA configuration may incorporate features or aspects from two or more of the BHA surrogates.

[0071] Continuing with the schematic flow chart 100 of FIG. 1, the methods of the present disclosure may optionally include drilling a well with a bottom hole assembly embodying the selected BHA configuration, such as represented by block 116. The drilling of the well may include forming the well to access a subsurface formation with the drilling equipment.

[0072] In some implementations, measured data may then be compared with calculated data and/or determined vibration performance indices for the selected BHA configuration, as shown in block 118. That is, as the drilling operations are being performed or at some time period following the drilling operations, sensors may be used to collect measured data associated with the operation of the drilling equipment. For example, the measured data may include but is not limited to RPM, WOB, axial, lateral, and stick/slip vibration measurements, drilling performance as determined by the Mechanical Specific Energy (MSE), or other appropriate derived quantities. Downhole data may be either transmitted to the surface in real-time or it may be stored in the downhole equipment and received when the equipment returns to the surface. The measured data and/or data derived from the measured data may be compared with calculated data and/or vibration performance indices from the modeling system for the selected BHA configuration.

[0073] The comparison of the measured data (or data derived from the measured data) with the model data and vibration performance indices can be used in a variety of manners, some examples of which are described in more detail herein. An illustrative and non-exhaustive list of such uses includes 1) updating the surrogate to better represent the BHA configuration; 2) updating the frequency-domain model to better simulate the response of the BHA during drilling operations under a variety of conditions; 3) updating the calculations and/or parameters used to determine one or more vibration performance indices; 4) updating the drilling operations plans for a selected bottom hole assembly configuration, such as represented by box 120 in FIG. 1; and 5) using measured vibration data to determine the model input excitation, simulating the response of the surrogates with this input, and comparing the model results with other measured data that is considered to be the system output response. The feedback process facilitates modeling validation and verification. It also helps to determine which of the vibration performance indices warrant greater weighting in the BHA configuration selection process, thus providing learning aids to advance the development of the BHA configuration selection process. Additionally or alternatively, the comparison between the model results and the measurements may enable the vibration performance indices to more accurately predict or indicate the vibrational tendencies of a BHA surrogate, such as by allowing one or more input parameters of a vibration performance index to be further refined or tuned. One example of such vibration performance index improvements includes weighting the various vibrational excitation modes to more accurately consider the modes that are most relevant.

[0074] Once the wellbore is formed, hydrocarbons may be produced from the well, as shown in block 122. The production of hydrocarbons may include completing the well with a well completion, coupling tubing between the well completion and surface facilities, and/or other known methods for extracting hydrocarbons from a wellbore. The process ends at block 124.

[0075] Beneficially, the present techniques may be utilized to design, construct, and/or utilize equipment that can reduce the impact of limiters that may hinder drilling operations. In some implementations, two or more BHA configurations may be compared simultaneously with concurrent calculation and display of model results for two or more surrogates. With this comparison, the merits of alternative BHA configurations can be evaluated. Further, in implementations where the calculated model data and the measured data are associated with the selected BHA configuration, other limiters that may be present during the drilling of the wellbore may be identified and addressed in a timely manner to further enhance drilling operations. For example, if the primary limiter appears to be torsional stick/slip vibrations and the sources of torque in the BHA due to contact forces have been minimized, another possible mitigator is to choose a less aggressive bit that generates less torque for a given applied weight on bit. An example of the modeling of two or more BHA configuration surrogates is described in greater detail below in FIG. 2.

[0076] FIG. 2 is an exemplary flow chart 200 of the modeling of two or more BHA surrogates in accordance with certain aspects of the present techniques. For exemplary purposes, in this flow chart, the modeling of the two or more BHA surrogates is described as being performed by a modeling system. The modeling system may include a computer system that operates a modeling program. The modeling program may include computer readable instructions or code that compares two or more BHA surrogates, which is discussed further below. While FIG. 2 is directed to the comparison of two or more BHA surrogates, the present methods and systems are useful in modeling a single BHA surrogate to identify operational and/or design parameters that can be modified to improve performance by reducing vibrations.

[0077] The flow chart 200 begins at block 202. To begin, the BHA layout and operating parameters are obtained for use in the modeling operations introduced above. At block 204, operating parameters may be obtained. The operating parameters, such as the anticipated ranges of WOB, RPM and wellbore inclination, may be obtained from a user entering the operating parameters into the modeling system or accessing a file having the operating parameters. For the static model, the condition of the BHA model end-point (e.g., end away from the drill bit) can be set to either a centered condition (e.g., the pipe is centered in the wellbore) or an offset condition (e.g., the pipe is laying on the low side of the wellbore).

[0078] The BHA design parameters are then obtained, as shown in block 206. The BHA design parameters may include available drill collar dimensions and mechanical properties, dimensions of available stabilizers, drill pipe dimensions, length, and the like. For example, if the drilling equipment is a section of tubing or pipe, the BHA design parameters may include the inner diameter (ID), outer diameter (OD), length and bending moment of inertia of the pipe, and the pipe material properties. Also, the modeling system may model drilling equipment made of steel, non-magnetic material, Monel, aluminum, titanium, etc. If the drilling equipment is a stabilizer or under-reamer, the BHA design parameters may include blade OD, blade length, and/or distance to the blades from the ends.

[0079] At block 208, the initial BHA surrogates are obtained. Obtaining of the BHA surrogates may include accessing a stored version of a previously modeled or utilized BHA configuration or BHA surrogate, interacting with the modeling system to specify or create a BHA surrogate from the BHA design parameters, or entering a proposed BHA configuration into the model that was provided by the drilling engineer or drilling service provider. The BHA surrogates specify the positioning of the equipment and types of equipment in the BHA, usually determined as the distance to the bit of each component.

[0080] Once the different BHA surrogates are obtained and/or constructed, the results for the selected BHA surrogates are calculated/modeled, as shown in block 210. The calculations may include calculation of the static states to determine force and tilt angle at the bit and static stabilizer contact forces, calculation of dynamic vibration performance indices, calculation of dynamic state values for specific excitation modes as a function of rotary speed, weight on bit, and distance to bit, and the like. More specifically, the calculations may include the dynamic lateral bending (e.g., flexural mode) and eccentric whirl dynamic response as perturbations about a static equilibrium, which may be calculated using the State Transfer Matrix method described below or other suitable method. This flexural or dynamic lateral bending mode may be referred to as “whirl.” The static responses may include the state vector response (e.g., displacement, tilt, bending moment, shear force, and contact forces or torques) as a function of distance from the bit, WOB, fluid density, and wellbore inclination (e.g., angle or tilt angle). For the dynamic response values, the state variables may be calculated as a function of distance from the drill bit, WOB, RPM, excitation mode, and end-lengths. For the lateral bending and eccentric whirl, the model states (e.g., displacement, tilt, bending moment, shear force, and contact forces or torques) may be calculated and displayed as functions of distance from the bit for specified WOB, RPM, excitation mode, and end-length.

[0081] As used herein, the “excitation mode” is the integer and/or non-integer multiple of the rotary speed or specific excitation frequency at which the system is being excited (for example, it is well known that a roller cone bit provides a three times multiple axial excitation, which may couple to the lateral mode). The “end-length” is the length of pipe added to the top of the BHA, often in the heavy-weight drillpipe, to evaluate the vibrational energy being transmitted uphole. Because the response may be sensitive to the location of the last nodal point, one computational approach is to evaluate a number of such possible locations for this nodal point for the purpose of computing the response. Then these different results may be averaged (by root-mean-square (RMS) or another averaging method) to obtain the overall system response for the parametric set of the various excitation modes and end-lengths for each RPM and WOB. Additionally or alternatively, the “worst case” maximum value may also be presented, which is described further below.

[0082] Once the results are calculated, the results are displayed as shown in block 210. When the present methods are implemented for direct comparison of two or more BHA surrogates, the results may be displayed simultaneously on one or more display screens and/or windows or may be displayed in a common window. As described above, the results may similarly be transmitted to remote locations for display or stored for later retrieval. The display may be on a screen or other audiovisual medium or may be printed. Additionally, the display may include graphical and/or numerical representations of the results.

[0083] Continuing with the flow chart of FIG. 2, the results are verified, as shown in block 212. The calculation result verification process may include determining by examination that, for example, there were no numerical problems encountered in the simulation and that all excitation modes were adequately simulated throughout the requested range of rotary speeds, bit weights, and end-lengths. In some implementations, the calculation result verification process may include discarding and/or discounting numerically divergent results in calculating one or more vibration performance indices. Other methods of verifying the results may be implemented.

[0084] At block 214, FIG. 2 illustrates that a determination may be made whether the BHA configurations represented by the surrogates and/or other parameters are to be modified. If the BHA configurations or specific parameters are to be modified, the BHA configurations and/or parameters may be modified in block 216. The modifications may include changing specific aspects in the operating parameters, BHA surrogates, BHA design parameters and/or adding a new BHA surrogate. As a specific example, the WOB, RPM and/or excitation mode may be changed to model another set of operating conditions. The BHA configurations and corresponding surrogates are typically adjusted by altering the distance between points of stabilization, by changing the sizes or number of stabilizers and drill collars, by relocating under-reamers or cross-overs to a different position in the BHA surrogate, and the like. Once the modifications are complete, the results may be recalculated in block 210, and the process may be iterated to further enhance performance.

[0085] However, if the BHA configurations and/or parameters are not to be modified, the results are provided, as shown in block 218. Providing the results may include storing the results in memory, printing a report of the results, and/or displaying the results on a monitor. For example, a side-by-side graphical comparison of selected BHA surrogates and/or preferred operating parameters may be displayed by the modeling system. The results of one or more of the calculated static and dynamic responses for specified WOB, RPM, excitation mode, end-lengths, and vibration indices may be displayed on two-dimensional or three-dimensional plots. Similarly, the results may be displayed as results for a single BHA surrogate, a comparison of results for two or more BHA surrogates, and/or a comparison of modeling results and measured data during actual drilling operations. While FIG. 2 illustrates that the method ends at block 220, additional steps may follow, such as the implementation of drilling operations incorporating the information learned during the methods of FIG. 2.

[0086] Beneficially, the modeling of the BHA surrogates may enhance drilling operations by providing a BHA more suitable to the drilling environment. For example, if one of the BHA surrogates is based on drilling equipment utilized in a certain field, then other surrogates may be modeled and directly compared with the previously utilized BHA surrogate. That is, one of the BHA surrogates may be used as a benchmark for comparing the vibration tendencies of other BHA surrogates. In this manner, the BHA surrogates may be compared, either simultaneously or as additional surrogates are modeled, to determine a BHA surrogate that reduces the effect of limiters, such as vibrations. To the extent that the modeling system is adapted to compare more than two different BHA surrogates, additional proposed BHA surrogates can be compared against each other or against a baseline surrogate. The comparative approach may be found to be more practical in some implementations. The relevant question to answer for the drilling engineer relates to which configuration of BHA components operates with the lowest vibrations over the operating conditions for a particular drilling operation. A preferred approach to address this design question is to model several alternative configurations and then select the one that performs in an optimal manner over the expected operating range or to operate the selected configuration at operating parameters suggested by the present methods. Such approach can be accomplished iteratively or through direct and simultaneous comparison of the several configurations.

Exemplary BHA Surrogates

[0087] As described above, BHA surrogates are representations of actual BHA configurations that can be input into the modeling systems to simulate the operation or response of the represented BHA configuration in a drilling operation. Accordingly, BHA surrogates, as representations of actual equipment, incorporate one or more assumptions and/or simplifications to allow the equipment to be mathematically modeled. As with most mathematical representations of actual equipment, the representation can be constructed in a variety of manners, some of which may be different but equal in application. Similarly, some of the different surrogate construction techniques may result in different surrogates that are more or less appropriate for different uses.

[0088] The present methods encompass the use of any suitable surrogate that can be used in a frequency-domain model of drilling operations to simulate drilling and associated vibrations. Exemplary surrogates include a lumped parameter surrogate and a distributed mass surrogate. In a lumped parameter surrogate, the BHA configuration is represented by point masses connected by massless beam and damper elements. In a distributed mass surrogate, the BHA configuration is represented by a beam having a distributed mass. Depending on the manner in which the BHA surrogate is constructed, the frequency-domain model(s) used to model the operation of the surrogate may vary, such as the selection of a 2D or a 3D frequency-domain model.

[0089] As suggested above, the BHA surrogates may be constructed in a variety of manners and the frequency-domain models may vary within the scope of the present disclosure. Through implementation of the present methods, it may be determined that one type of surrogate and/or one type of frequency-domain model more accurately represents actual drilling operations for a particular BHA configuration, for particular operating conditions, or for particular environments. For example, it may be found that 2D lumped parameter surrogates and associated modeling results correspond sufficiently closely to measured data for a particular BHA configuration or drilling application. As another example, it may be found that 3D distributed mass surrogates and associated frequency-domain modeling results more closely correspond to measured data for a particular type of vibration or for a particular excitation mode. Accordingly, methods within the scope of the present disclosure include methods where different BHA surrogates and different frequency-domain models are used to represent one or more BHA configurations in a single drilling operation. Additionally or alternatively, mathematical combinations of different surrogates and/or frequency-domain models may be used to improve the accuracy of the modeling results as compared to measured data.

Exemplary Lumped Parameter BHA Vibration Models

[0090] As an example, one exemplary implementation of a BHA vibration model is described. However, it should be noted that other BHA models, for example using one or more of the calculation methods discussed above, may also be used to form a comparative vibration performance index in a similar manner. As used herein, “BHA vibration model” refers to the use of a BHA surrogate and associated frequency-domain modeling principles to model or simulate the vibrations of a drilling operation using the BHA configuration represented by the BHA surrogate. These methods may include but are not limited to two-dimensional or three-dimensional finite element modeling methods. For example, calculating the results for one or more BHA configurations may include generating a surrogate or mathematical model for each BHA configuration; calculating the results of the surrogate for specified operating parameters and boundary conditions; identifying the displacements, tilt angle (first spatial derivative of displacement), bending moment (calculated from the second spatial derivative of displacement), and beam shear force (calculated from the third spatial derivative of displacement) from the results of the surrogate simulation; and determining state vectors and matrices from the identified outputs of the surrogate simulation. In more complex models, these state vectors may be assigned at specific reference nodes, for example at the neutral axis of the BHA cross-section, distributed on the cross-section and along the length of the BHA, or at other convenient reference locations. As such, the state vector response data, calculated from the finite element model results, may then be used to calculate vibration performance indices to evaluate BHA configurations and to compare with alternative BHA configurations, as described herein.

[0091] The BHA vibration model described in this section is a lumped parameter model, which is one embodiment of a mathematical model, implemented within the framework of state vectors and transfer function matrices. The state vector represents a complete description of the BHA system response at any given position in the BHA surrogate, which is usually defined relative to the location of the bit. The transfer function matrix relates the value of the state vector at one location with the value of the state vector at some other location. The total system state includes a static solution plus a dynamic perturbation about the static state. The linear nature of the model for small dynamic perturbations facilitates static versus dynamic decomposition of the system. The dynamic model presented in this section is one variety in the class of forced frequency response models, with specific matrices and boundary conditions as described below. Other dynamic models may be developed for BHA vibration models utilizing alternative BHA surrogates and/or alternative operating parameters.

[0092] Transfer function matrices may be multiplied to determine the response across a series of elements in the model. Thus, a single transfer function can be used to describe the dynamic response between any two points. A lumped parameter model yields an approximation to the response of a continuous system. Discrete point masses in the BHA surrogate are connected by massless springs and/or dampers to other BHA surrogate mass elements and, in one variation, to the wellbore at points of contact by springs and, optionally, damper elements. The masses are free to move laterally within the constraints of the applied loads, including gravity.

[0093] Matrix and State Vector Formulation

[0094] For lateral motion of a lumped parameter model in a plane, the state vector includes the lateral and angular deflections, as well as the beam bending moment and shear load. The state vector u is extended by a unity constant to allow the matrix equations to include a constant term in each equation that is represented. The state vector u may then be written as equation (e1) as follows:

[0095] $\begin{array}{cc}u=\left(\begin{array}{c}y\\ \theta \\ M\\ V\\ 1\end{array}\right)& \left(e1\right)\end{array}$
Where y is lateral deflection of the beam from the centerline of the assembly, θ is the angular deflection or first spatial derivative of the displacement, M is the bending moment that is calculated from the second spatial derivative of the displacement, and V is the shear load of the beam that is calculated from the third spatial derivative of the displacement. For a three-dimensional model, the state vector defined by equation (e1) may be augmented by additional states to represent the displacements and derivatives along an orthogonal axis at each node. The interactions between the motions at each node may, in the general case, include coupled terms.

[0096] By linearity, the total response may be decomposed into a static component u^s and a dynamic component u^d (e.g., u=u^s+u^d).

[0097] In the forced-frequency response methods, the system is assumed to oscillate at the frequency ω of the forced input, which is a characteristic of linear systems. Then, time and space separate in the dynamic response and, using superposition, the total displacement of the beam at any axial point x for any time t may be expressed by the equation (e2):
u(x,t)=u^s(x)+u^d(x) sin (ωt)  (e2)

[0098] State vectors u_i (for element index i ranging from 1 to N) may be used to represent the state of each mass element, and the state vector u₀ is used to designate the state at the bit. Transfer function matrices are used to relate the state vector u_i of one mass element to the state u_i−1 of the preceding mass element. If there is no damping in the model, then the state vectors are real-valued. However, damping may be introduced and then the state vectors may be complex-valued, with no loss of generality.

[0099] Because state vectors are used to represent the masses, each mass may be assumed to have an associated spring and/or damper connecting it to the preceding mass in the model. With the notation M_i denoting a mass transfer matrix, and a beam bending element transfer matrix represented by B_i the combined transfer function T_i is shown by the equation (e3) below.
T_i=M_iB_i  (e3)
Numerical subscripts are used to specify each mass-beam element pair. For example, the state vector u₁ may be calculated from the state u_o represented by the equation (e4).
u₁=M₁B₁u₀=T₁u₀, and thus u_i=T_iu_i−1  (e4)
These matrices can be cascaded to proceed up the BHA to successive locations. For example, the state vector u₂ may be represented by the equation (e5).
u₂=T₂u₁=T₂T₁u₀  (e5)
While continuing up to a contact point, the state vector u_N may be represented by equation (e6).
u_N=T_Nu_N−1=T_NT_N−1 . . . T₁u₀  (e6)

[0100] Accordingly, within an interval between contact points, the state u_j at any mass element can be written in terms of any state below that element u_i using a cascaded matrix S_ij times the appropriate state vector by the equation (e7):
u_j=S_iju_i where for i<j, S_ij=T_jT_j−1 . . . T_i+1  (e7)
Consideration of the state vector solution at the contact points will be discussed below.

[0101] Formulation of Mass Matrices

[0102] The mass transfer function matrix for the static problem is derived from the balance of forces acting on a mass element m. Generally, each component of the BHA is subdivided into small elements, and this lumped mass element is subjected to beam shear loads, gravitational loading (assuming inclination angle φ), wellbore contact with a stiffness k, and damping force with coefficient b. The general force balance for the element may be written as equation (e8) using the “dot” and “double dot” notations to represent the first and second time derivatives, or velocity and acceleration, respectively.
mÿ=V_i−V_i−1−mg sin φ−ky−b{dot over (y)}=0  (e8)

[0103] The lumped mass element transfer function matrix under static loading includes the lateral component of gravity (mg sin φ) and either a contact spring force or, alternatively, a constraint applied in the solution process, in which case the value of k is zero. In the static case, the time derivatives are zero, and thus inertial and damping forces are absent. The static mass matrix may be written as the following equation (e9).

[0104] $\begin{array}{cc}{M}_S=\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ k& 0& 0& 1& \left(\mathrm{mg}\sin \varphi \right)\\ 0& 0& 0& 0& 1\end{array}\right)& \left(e9\right)\end{array}$

[0105] In lateral dynamic bending, the forces applied to the mass consist of the beam shear forces, wellbore contact, and damping loads. Again, the wellbore contact may be either the result of a spring force or an applied constraint relation. However, because the dynamic perturbation about the static state is sought (using the principle of linear superposition), the gravitational force is absent from the dynamic mass matrix.

[0106] In the dynamic example, the applied loads may be unbalanced, leading to an acceleration of the mass element. The mass times lateral acceleration equals the force balance of the net shear load, spring contact, and damping forces, resulting in the equation (e10).
mÿ=V_i−V_i−1−ky−b{dot over (y)}  (e10)
Assuming a complex harmonic forced response y^d˜e^iωt, where i represents the imaginary number equal to √{square root over (−1)}, the solution to equation (e10) may be found in equation (e11).
V_i=V_i−1+(k+ibω−mω²)  (e11)

[0107] The lumped mass element transfer function matrix M_B, for the lateral bending mode dynamic perturbation, is then written by the following equation (e12).

[0108] $\begin{array}{cc}{M}_B=\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ \left(k+ib\omega -m{\omega }^2\right)& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right)& \left(e12\right)\end{array}$

[0109] The mass matrix in the dynamic whirl model involves a constant-magnitude force which resembles the gravitational force in the static mass matrix. It is assumed that each drill collar has a slightly unbalanced mass, generating a centrifugal force proportional to this unbalanced mass times the square of the rotational frequency. For a small value ε which represents the dimensionless off-axis distance of the unbalanced mass, the equation of motion for forced response is given by equation (e13).
mÿ=V_i−V_i−1+εmω²−ky−b{dot over (y)}  (e8)

[0110] The radial displacement does not change with time for this simplified whirl mode example, and thus the acceleration and velocity may be set to zero. This represents a steady rotational motion, not unlike a rotating gravitational load, in contrast to the lateral bending mode in which the displacement oscillates through a zero value. The resulting whirl matrix is represented in equation (e14).

[0111] $\begin{array}{cc}{M}_W=\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ k& 0& 0& 1& \left(\varepsilon m{\omega }^2\right)\\ 0& 0& 0& 0& 1\end{array}\right)& \left(e14\right)\end{array}$

[0112] The value ε may take either positive or negative signs in order to represent the shape of the whirl response being modeled. The first whirl mode is generally represented by alternating signs on successive intervals of BHA components as one proceeds up the borehole.

[0113] The lumped parameter mass m is defined as the mass of the element piece of the respective BHA component. In addition, the mass of the drill collar, pipe, or other BHA component is effectively increased by the drilling fluid contained within the collar and that which is entrained by the BHA element as it vibrates. The technique of “added mass” may be used to approximate this phenomenon. For this purpose, a crude approximation is to increase the dynamic collar mass by 10%, leading to a slight reduction in natural frequency. This is a representative value only, and calibration of model results with field data may indicate alternative values for the “added mass” effect that may be used in the model. Note that it is not appropriate to apply the added mass to the static solution. As noted above, depending on the solution method, the spring constant may be omitted if the solution is to apply a constraint relationship such that the BHA model is not permitted to extend outside the wellbore by more than a very small amount.

[0114] If the constraint model is not used, then the contact stiffness k in the above relations should be included explicitly. In this example, a factor to be considered in the choice of wellbore contact stiffness k when modeling dynamic excitation is that the value of k should be chosen sufficiently high for the mass m such that the natural frequency √{square root over (k/m)} is greater than the maximum excitation frequency ω to be evaluated, so that resonance due to this contact representation is avoided. Thus, for an excitation mode of n times the rotary speed, the contact stiffness k may be greater than m(nω)² (e.g., k>m(nω)²).

[0115] Alternatively, and in the preferred embodiment, compliance at the points of contact between BHA and wellbore may be neglected and a fixed constraint relationship applied in the solution method, with k=0 in the matrices above. This approach is described further below.

[0116] Formulation of Stiffness Matrix

[0117] The Euler-Bernoulli beam bending equation for a uniform beam with constant Young's modulus E, bending moment of inertia I, and axial loading P may be written as the fourth-order partial differential equation (e15).

[0118] $\begin{array}{cc}\mathrm{EI}\frac{{\partial }^{4}y}{\partial x^4}-P\frac{{\partial }^{2}y}{\partial x^2}=0& \left(e15\right)\end{array}$

[0119] The characteristic equation for the general solution is represented by equation (e16)
y=e^βx  (e16)
This equation expresses the lateral displacement as the exponential power of a parameter β times the distance x from a reference point, in which the term β is to be found by replacing this solution in equation (e15) and solving with equations (e17) and (e18) below.

[0120] $\begin{array}{cc}{\beta }^2\left({\beta }^2-\frac{P}{\mathrm{EI}}\right)=0& \left(e17\right)\\ \beta =0,\pm \sqrt{\frac{P}{\mathrm{EI}}}& \left(e18\right)\end{array}$

[0121] Note that β is either real (beam in tension), imaginary (beam in compression), or 0 (no axial loading). The appropriate particular solution is a constant plus linear term in x. Thus, the displacement of an axially loaded beam may be represented by the equation (e19).
y=a+bx+ce^βx+de^−βx  (e19)
where the constants a, b, c, and d are found by satisfying the boundary conditions.

[0122] The remaining components of the state vector are determined by the following equations in the spatial derivatives of lateral displacement with the axial coordinate x (e20).

[0123] $\begin{array}{cc}\theta =\frac{\partial y}{\partial x}M=\mathrm{EI}\frac{{\partial }^{2}y}{\partial x^2}V=-\mathrm{EI}\frac{{\partial }^{3}y}{\partial x^3}& \left(e20\right)\end{array}$

[0124] The resulting beam bending stiffness transfer function matrix B may be represented by the following equation (e21).

[0125] $\begin{array}{cc}B=\left(\begin{array}{ccccc}1& L& \left(\frac{-2+e^{\beta L}+e^{-\beta L}}{2P}\right)& \left(\frac{2\beta L-e^{\beta L}+e^{-\beta L}}{2P\beta }\right)& 0\\ 0& 1& \left(\frac{e^{\beta L}-e^{-\beta L}}{2\beta \mathrm{EI}}\right)& \left(\frac{2-e^{\beta L}-e^{-\beta L}}{2P}\right)& 0\\ 0& 0& \left(\frac{e^{\beta L}+e^{-\beta L}}{2}\right)& \left(\frac{-e^{\beta L}+e^{-\beta L}}{2\beta }\right)& 0\\ 0& 0& \left(\frac{-\beta e^{\beta L}+{\mathrm{\beta e}}^{-\beta L}}{2}\right)& \left(\frac{e^{\beta L}+e^{-\beta L}}{2}\right)& 0\\ 0& 0& 0& 0& 1\end{array}\right)& \left(e21\right)\end{array}$

[0126] Boundary Conditions and System Excitation

[0127] With the mass and beam element transfer functions defined, the boundary conditions and system excitation are determined to generate frequency-domain model predictions. Separate boundary conditions are used to model the static bending, dynamic lateral bending, and eccentric whirl problems.

[0128] In each of these examples of lumped parameter BHA vibration models, the solution proceeds from the bit to the first stabilizer or other contact point, then from the first stabilizer to the second stabilizer or other contact point, and so on, proceeding uphole one solution interval at a time (e.g., from the bit as the starting interval). Finally, the interval from the top contact point to the end point is solved. As suggested, contact points are often provided by stabilizers, but may be provided by other BHA components, such as an under-reamer, or perhaps even by contact of one or more BHA components at intermediate points between specific contact points, such as drill collars resting on the wall between stabilizers. For convenience and brevity herein, the exemplary stabilizer will be used to refer to the variety of BHA components that may provide a contact point. The end point is the upper node of the BHA model, and it may be varied to consider different possible nodal points at the “end-length.” An appropriate lateral displacement for this end point is assumed in the static model, based on the amount of clearance between the pipe and the wellbore.

[0129] In these methods, the states in each solution interval are determined by three conditions at the lower element (bit or bottom stabilizer in the interval), and one condition at the upper element (end point or top stabilizer in the interval). With these four conditions and the overall matrix transfer function from the lower to the upper element, the remaining unknown states at the lower element may be calculated.

[0130] Beginning at the bit, the displacement of the first stabilizer is used to determine the bit state, and thus all states up to the first stabilizer are determined using the appropriate transfer function matrices. By continuity, the displacement, tilt, and moment are now determined at the first stabilizer point of contact. The beam shear load is undetermined, as this state does not have a continuity constraint because there is an unknown side force acting between the stabilizer and the wellbore. The displacement of the next stabilizer is used to provide the fourth condition necessary to obtain the solution over the next interval, and thus the complete state at the stabilizer is determined. The contact force between stabilizer and wellbore may be calculated as the difference between this state value and the prior shear load calculation from the previous BHA section. Using the cascaded matrix formulation in equation (e22)

[0131] $\begin{array}{cc}\left(\begin{array}{c}{y}_j\\ {\theta }_j\\ M_j\\ V_j\\ 1\end{array}\right)=S_{\mathrm{ij}}\left(\begin{array}{c}{y}_i\\ {\theta }_i\\ M_i\\ V_i\\ 1\end{array}\right)\mathrm{with}\mathrm{the}\mathrm{conditions}\left(\begin{array}{c}{V}_i·\mathrm{unknown}\\ y_j=0\end{array}\right)& \left(e22\right)\end{array}$
Then the unknown shear load at the lower stabilizer is calculated using an equation (e23) to obtain zero displacement at the upper position.
0=S₁₁y_i+S₁₂θ_i+S₁₃M_i+S₁₄V_i+S₁₅  (e23)

[0132] The beam shear load is discontinuous across the contact points, and the sideforce at such a node may be calculated as the difference between the value obtained by propagating the states from below, V_i⁻, and the value calculated to satisfy the constraint relation for the next segment, V_i⁺. Therefore, the contact sideforce may be represented by the equation (e24).
F_i=V_i⁺−V_i⁻  (e24)

[0133] For the static example, the tilt and sideforce are unknown at the bit. A trial bit tilt angle is used to generate a response and the state vectors are propagated uphole from one contact point to the next, finally reaching the end-point. The final value for the bit tilt angle and sideforce are determined by iterating until the appropriate end condition is reached at the top of the model, for instance a condition of tangency between the pipe and borehole wall. Alternatively, the solution may start at an uphole point of tangency and proceed downhole to the bit, iterating to convergence on the bit condition by varying the distance to the point of tangency. Other convergence methods may also be selected.

[0134] Additionally or alternatively, the stabilizer configuration or the configuration of another element of the BHA surrogate may suggest an additional constraint. For example, for full-gauge stabilizers, it may be appropriate to further constrain the BHA vibration model to a tilt angle of zero. Such a constraint may be appropriate due to the interactions between the full-gauge stabilizers and the wellbore wall. A BHA vibration model incorporating this additional constraint will result in an increase in the reaction sideforces in the model and a discontinuity in the bending moment will occur to represent the reaction torque due to the fixed tilt angle constraint. Additionally, equation (e22) above will be modified with the additional constraint as equation (e22′) for the node j where the tilt constraint applies.

[0135] $\begin{array}{cc}\left(\begin{array}{c}{y}_j\\ {\theta }_j\\ M_j\\ V_j\\ 1\end{array}\right)=S_{\mathrm{ij}}\left(\begin{array}{c}{y}_i\\ {\theta }_i\\ M_i\\ V_i\\ 1\end{array}\right)\mathrm{with}\mathrm{the}\mathrm{conditions}\left(\begin{array}{c}{y}_j=0\\ {\theta }_j=0\\ M_j\mathrm{unknown}\\ V_i\mathrm{unknown}\end{array}\right)& \left(e{22}^{\prime }\right)\end{array}$

[0136] For the dynamic models (flexural bending, whirl, and twirl), a reference bit excitation sideforce is applied, e.g., V_bit=const. The first stabilizer is assumed to be constrained by the pinned condition (e22) or by the built-in condition (e22′). If applying (e22), two more conditions are specified to uniquely solve the equations from the bit up to the first stabilizer. One choice for the boundary conditions is to assume that for small lateral motion, the tilt and moment at the bit are zero. This set of boundary conditions may be written as shown in equation (e25):
y_stab=θ_bit=M_bit=0 V_bit=const  (e25)

[0137] An alternate set of boundary conditions may be considered by assuming that the tilt angle at the first stabilizer is zero as in (e22′), equivalent to a cantilevered condition. One selection for the remaining constraint is to assume that there is no moment at the bit. This alternate set of boundary conditions may be written as shown in equation (e26):
y_stab=θ_stab=M_bit=0 V_bit=const  (e26)

[0138] As is well understood, the bit may be excited in a variety of manners leading to dynamic vibrations. The bit excitation by way of an applied sideforce described above is one common manner. The present methods may be adapted to provide BHA vibration models for other forms of bit excitation as well. As one exemplary modification to accommodate or consider an alternative form of bit excitation, the bit may be excited by an applied moment at the bit, such as may occur in drilling operations when the bit is penetrating a laminated formation. The BHA vibration model may be run utilizing an applied excitation frequency for an applied moment at the bit at any multiple of the rotary speed. In some implementations, it may be preferred to run the model multiple times utilizing various multiples of the rotary speed and considering averages, maximums, and/or historical/measured data to provide a more robust and/or accurate model. For example, while a 1× multiple of the rotary speed may be the most likely excitation frequency, the model may be run using various multiples, including non-integer multiples such as 1.5×, 1.75×, etc. Additionally and alternatively, a fixed excitation frequency can be applied to the bit to represent certain excitation sources that are constant in frequency and not multiples of the rotary speed. One example is the drilling fluid pressure which has pressure pulses in accordance with the stroke rate of the mud pumps. These pulses may cause lateral motion at the bit due to the time-varying pressure drop through the bit nozzles.

[0139] When the excitation is modeled as an applied moment at the bit, equations (e25) and (e26) above will be changed to correspond to the changed excitation mode. For the boundary conditions identified above in connection with equations (e25) and (e26), the boundary conditions for applied moment at the bit may be written as in equations (e25′) and (e26′), as seen below.
y_stab=y_bit=V_bit=0 M_bit=const  (e25′)
y_stab=θ_stab=V_bit=0 M_bit=const  (e26′)

[0140] Regardless of the form of excitation applied, which may include one or more of those described above and/or other common excitations, the solution marches uphole one stabilizer at a time. The solution, or rather the implementation of the model, terminates at the last node of the BHA surrogate, which is arbitrarily chosen but may be located at different “end-lengths” in the dynamic case. By selecting different end-lengths and RMS-averaging the results, vibration performance indices may be formed that are robust. To guard against strong resonance at an individual nodal point, the maximum result may be examined as well, and conversely, the minimum value may be examined to evaluate possible preferred operating regions. These techniques of RMS-averaging and examining the maximum may be preferred when determining a vibration performance index that is sensitive to the selection of a nodal point location. For example, an end-point curvature index and other point indices may be more sensitive to nodal point locations than the interval indices described herein. Being less sensitive to the end location and condition, the interval indices may be preferred in some implementations. It should further be noted, as indicated above, that BHA contact with the borehole at locations between stabilizers may optionally be treated as a nodal point in this analysis method, and the solution propagation modified accordingly.

Exemplary Distributed Mass BHA Vibration Models

[0141] As introduced above, bottom hole assemblies can be represented by surrogates of a variety of construction methods. The exemplary model described above considered a lumped parameter BHA vibration model in great detail. As an illustration of other BHA surrogates and associated frequency-domain models that may be used in BHA vibration models within the scope of the present disclosure, an exemplary distributed mass BHA vibration model will now be described with reference to the discussions above.

[0142] As with the discussion above, the system state vector for a distributed mass BHA vibration model is written as equation (e27):

[0143] $\begin{array}{cc}u=\left(\begin{array}{c}y\\ \theta \\ M\\ V\\ 1\end{array}\right)=\left(\begin{array}{c}\mathrm{lateral}\mathrm{deflection}\\ \mathrm{angular}\mathrm{deflection}\\ \mathrm{bending}\mathrm{moment}\\ \mathrm{shear}\mathrm{load}\\ \mathrm{unity}\mathrm{constant}\end{array}\right)& \left(e27\right)\end{array}$
By linearity, the total response may be decomposed into static and dynamic components. In the forced-frequency response method, the system is assumed to oscillate at the frequency of the forced input; this is a characteristic of linear systems. Then time and space separate in the dynamic response, and, using superposition, one may write the state vector as a function of time and space as equation (e28):
u(x,t)=u^s(x)+u^d(x)sin(ωt)  (e28)

[0144] While many of the principles from the above lumped parameter discussion are relevant to this distributed mass example, several of the factors and relationships described above depended on the mass of the BHA as represented in the surrogates. As the distributed mass BHA surrogate does not simplify the BHA configuration as point mass and springs and/or dampers, several of the relationships described above are adapted to correspond with the construction of the BHA surrogate, as will be seen below. The relations to describe the beam deflections depend on the beam properties (E and I) and the distributed weight per unit length, W. The axial load, P, is also a factor that translates directly from the lumped parameter model to the distributed mass model. As with the above example, the present example will consider both the static case and the dynamic case (or dynamic perturbations around the static solution).

[0145] Static Case Solution

[0146] Considering first the impact of the distributed mass on the static solution described above for the lumped parameter model, a primary difference between the lumped parameter model and the distributed mass model lies in the transfer matrices used in the models. In the above discussion the BHA surrogate was represented by both a mass transfer matrix and a beam bending element transfer matrix (see equation (e3)). However, in the distributed mass models, the mass is distributed along the length of the beam and the two elements (mass and beam bending) can be considered together in a single transfer function, as seen below. Additionally, the relevant mass effect is the component of gravity orthogonal to the axis of the wellbore. Accordingly, it is necessary to adjust the material weight per unit length by the sine of the inclination angle, φ. Therefore, using the term W=(−ρAg) for density ρ, cross-sectional area A, and gravitational constant g, equation (e15) from above is modified as equation (e29):

[0147] $\begin{array}{cc}\mathrm{EI}\frac{d^{4}y}{d{x}^4}-P\frac{d^{2}y}{d{x}^2}-W\sin \left(\phi \right)=0& \left(e29\right)\end{array}$

[0148] Assuming an exponential solution to the homogeneous equation of the form e^βx (see equation (e16) above) yields a characteristic equation with the solution for β of the form shown above as equation (e18) and repeated here for convenience as equation (e30):

[0149] $\begin{array}{cc}\beta =0,\pm \sqrt{\frac{P}{EI}}& \left(e30\right)\end{array}$

[0150] The term β is either real (beam in tension), imaginary (beam in compression), or zero (no axial loading). The particular solution is the sum of linear and quadratic terms in x plus a constant, and the homogeneous solution includes exponential functions with both possible values for β. Thus, the displacement of an axially loaded beam may be represented by the equation (e31):
y=ax²+bx+c+de^βx+fe^−βx  (e31)

[0151] As before, the derivatives can be identified relative to the system state variables,

[0152] $\begin{array}{cc}\theta =\frac{\partial y}{\partial x}M=\mathrm{EI}\frac{{\partial }^{2}y}{\partial x^2}V=-\mathrm{EI}\frac{{\partial }^{3}y}{\partial x^3}& \left(e32\right)\end{array}$

[0153] The matrix that relates the state vector at x=0 to the state at x=L for the lateral beam bending of a distributed mass with axial loading is then written as equation (e33). Here, we denote the matrix as T to identify the distributed mass matrix with the combined mass and stiffness matrix as shown in (e3). The subsequent matrix operations to obtain the solution then follow as described above, with a simple change in the matrix calculations to reflect the distributed mass model surrogates.

[0154] $\begin{array}{cc}{T}_S=\left(\begin{array}{ccccc}1& L& \left(\frac{\begin{array}{c}-2+e^{\beta L}+\\ e^{-\beta L}\end{array}}{2P}\right)& \left(\frac{\begin{array}{c}2\beta L-e^{\beta L}+\\ e^{-\beta L}\end{array}}{2P\beta }\right)& \left(\frac{\begin{array}{c}-{\left(\beta L\right)}^2-2+\\ e^{\beta L}+e^{-\beta L}\end{array}}{2P{\beta }^2}\right)·\\ & & & & W\sin \left(\phi \right)\\ 0& 1& \left(\frac{e^{\beta L}-e^{-\beta L}}{2\beta \mathrm{EI}}\right)& \left(\frac{2-e^{\beta L}-e^{-\beta L}}{2P}\right)& \left(\frac{\begin{array}{c}-2\beta L+\\ e^{\beta L}-e^{-\beta L}\end{array}}{2P\beta }\right)·\\ & & & & W\sin \left(\phi \right)\\ 0& 0& \left(\frac{e^{\beta L}+e^{-\beta L}}{2}\right)& \left(\frac{-e^{\beta L}+e^{-\beta L}}{2\beta }\right)& \left(\frac{-2+e^{\beta L}+e^{-\beta L}}{2{\beta }^2}\right)·\\ & & & & W\sin \left(\phi \right)\\ 0& 0& \left(\frac{\begin{array}{c}-\beta e^{\beta L}+\\ {\mathrm{\beta e}}^{-\beta L}\end{array}}{2}\right)& \left(\frac{e^{\beta L}+e^{-\beta L}}{2}\right)& \left(\frac{-e^{\beta L}+e^{-\beta L}}{2\beta }\right)·\\ & & & & W\sin \left(\phi \right)\\ 0& 0& 0& 0& 1\end{array}\right)& \left(e33\right)\end{array}$

[0155] For a beam in compression, P is negative, β is imaginary, and we have equations (e34):

[0156] $\begin{array}{cc}\begin{array}{ccc}\beta =i\lambda & \lambda =\pm \sqrt{\frac{-P}{\mathrm{EI}}}& {\lambda }^2=\frac{-P}{\mathrm{EI}}\end{array}\begin{array}{cc}{e}^{\mathrm{i\lambda }L}+e^{-i\lambda L}=2\cos \left(\lambda L\right)& e^{i\lambda L}-e^{-i\lambda L}=2i\sin \left(\lambda L\right)\end{array}& \left(e34\right)\end{array}$
and the beam matrix equation reduces to equation (e35):

[0157] $\begin{array}{cc}{T}_{S,\mathrm{COMP}}=(\begin{array}{ccccc}1& L& \left(\frac{\cos \left(\lambda L\right)-1}{P}\right)& \left(\frac{\begin{array}{c}\lambda L-\\ \sin \left(\lambda L\right)\end{array}}{P\lambda }\right)& \left(\frac{1-\cos \left(\lambda L\right)}{{\lambda }^2}-\frac{L^2}{2}\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 1& \left(\frac{\sin \left(\lambda L\right)}{\lambda \mathrm{EI}}\right)& \left(\frac{1-\cos \left(\lambda L\right)}{P}\right)& \left(\frac{\sin \left(\lambda L\right)-\lambda L}{\lambda }\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 0& \cos \left(\lambda L\right)& \left(\frac{-\sin \left(\lambda L\right)}{\lambda }\right)& \mathrm{EI}·\left(\cos \left(\lambda L\right)-1\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 0& \lambda \sin \left(\lambda L\right)& \cos \left(\lambda L\right)& \lambda ·\mathrm{EI}·\sin \left(\lambda L\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 0& 0& 0& 1\end{array})& \left(\mathrm{e35}\right)\end{array}$

[0158] For a beam in tension, β is real-valued and we have equations (e36):
e^βL+e^−βL=2 cos h(βL) e^βL−e^−βL=2 sin h(βL)  (e36)
and the beam matrix equation reduces to equation (e37):

[0159] $\begin{array}{cc}{T}_{S,\mathrm{TENS}}=(\begin{array}{ccccc}1& L& \left(\frac{\cosh \left(\beta L\right)-}{P}\right)& \left(\frac{\begin{array}{c}\beta L-\\ \sinh \left(\beta L\right)\end{array}}{P\beta }\right)& \left(\frac{\cosh \left(\beta L\right)-1}{{\beta }^2}-\frac{L^2}{2}\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 1& \left(\frac{\sinh \left(\beta L\right)}{\beta \mathrm{EI}}\right)& \left(\frac{\begin{array}{c}1-\\ \cosh \left(\mathrm{\beta L}\right)\end{array}}{P}\right)& \left(\frac{\sinh \left(\beta L\right)-\beta L}{\beta }\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 0& \cosh \left(\beta L\right)& \left(\frac{-\sinh \left(\beta L\right)}{\beta }\right)& \mathrm{EI}·\left(\cosh \left(\beta L\right)-1\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 0& -\beta \sinh \left(\beta L\right)& \cosh \left(\beta L\right)& \left(-\beta \right)·\mathrm{EI}·\sinh \left(\beta L\right)·\frac{W\sin \left(\phi \right)}{P}\\ 0& 0& 0& 0& 1\end{array})& \left(\mathrm{e37}\right)\end{array}$

[0160] For a beam with no axial loading, the differential equation is simplified because the term involving P drops out. The solution is a fourth-order polynomial, and the corresponding matrix result is the following equation (e38):

[0161] $\begin{array}{cc}{T}_{S,\mathrm{ZERO}}=\left(\begin{array}{ccccc}1& L& \frac{L^2}{2\mathrm{EI}}& \left(-1\right)·\frac{L^3}{6\mathrm{EI}}& \frac{W\sin \left(\phi \right)·L^4}{24\mathrm{EI}}\\ 0& 1& \frac{L}{\mathrm{EI}}& \left(-1\right)·\frac{L^2}{2\mathrm{EI}}& \frac{W\sin \left(\phi \right)·L^3}{6\mathrm{EI}}\\ 0& 0& 1& \left(-L\right)& \frac{W\sin \left(\phi \right)·L^2}{2}\\ 0& 0& 0& 1& \left(-L\right)·W\sin \left(\phi \right)\\ 0& 0& 0& 0& 1\end{array}\right)& \left(\mathrm{e38}\right)\end{array}$

[0162] Having identified the beam matrices for the different conditions under which the BHA surrogate may be placed during the simulations, the matrices may be used to calculate the static solutions using methods analogous to those described above.

[0163] Dynamic Bending Solution

[0164] Turning now to consider the dynamic perturbation about the static solution, by separation of variables, the total displacement is the product of a function of space and a function of time, as seen in equation (e39):
u^d(x,t)=y(x)τ(t)  (e39)
The components of the total displacement that are a function of time can be further described by equations (e40):

[0165] $\begin{array}{cc}\begin{array}{cc}\frac{d^2{\tau }_1}{d{t}^2}=-{\omega }_n^2\tau & \tau \left(t\right)=A\cos \left({\omega }_{n}t\right)+B\sin \left({\omega }_{n}t\right)\end{array}& \left(e40\right)\end{array}$

[0166] With reference to the above discussions as a framework, the dynamic equation for an interval with a constant axial load P, weight per unit length W, and gravitational constant g may be written as equation (e41):

[0167] $\begin{array}{cc}\mathrm{EI}\frac{d^{4}y}{d{x}^4}-P\frac{d^{2}y}{d{x}^2}-\frac{W}{g}{\omega }_n^{2}y=0& \left(e41\right)\end{array}$

[0168] Writing y as an exponential function of x, the characteristic polynomial is equation (e42):

[0169] $\begin{array}{cc}{r}^4-\frac{P}{\mathrm{EI}}{r}^2-\frac{W{\omega }_n^2}{g·EI}=0& \left(e42\right)\end{array}$

[0170] This fourth-order equation has two solutions, κ and λ, shown in equations (e43):

[0171] $\begin{array}{cc}{\kappa }^2={\left[\frac{P^2}{4{\left(\mathrm{EI}\right)}^2}+\frac{W{\omega }_n^2}{g·\mathrm{EI}}\right]}^{1/2}+\frac{P}{2\mathrm{EI}}{\lambda }^2={\left[\frac{P^2}{4{\left(\mathrm{EI}\right)}^2}+\frac{W{\omega }_n^2}{g·\mathrm{EI}}\right]}^{1/2}-\frac{P}{2\mathrm{EI}}& \left(\mathrm{e43}\right)\end{array}$
The solution to the equation can be given in general form as equation (e44):
y(x)=c₁ cos h(κx)+c₂ sin h(κx)+c₃ cos(λx)+c₄ sin(λx)  (e44)

[0172] Within an element, the position x=0 can be chosen at one face, and the opposite face is then at x=L. At the origin, the cosine functions have unity value, and the sine functions are zero, which can be presented in normalized state vector form as equation (e45):

[0173] $\begin{array}{cc}\left(\begin{array}{c}{y}_o\\ {\theta }_o\\ M_o\\ V_o\end{array}\right)=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& \kappa & 0& \lambda \\ {\kappa }^2\mathrm{EI}& 0& \left(-{\lambda }^2\right)\mathrm{EI}& 0\\ 0& \left(-{\kappa }^3\right)\mathrm{EI}& 0& {\lambda }^3\mathrm{EI}\end{array}\right]·\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\end{array}\right)& \left(\mathrm{e45}\right)\end{array}$
This matrix is invertible, so the coefficients can be solved for in terms of the state vector at one face of the element, as shown in equation (e46).

[0174] $\begin{array}{cc}\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\end{array}\right)=\left(\frac{1}{{\lambda }^2+{\kappa }^2}\right)\left[\begin{array}{cccc}{\lambda }^2& 0& \frac{1}{\mathrm{EI}}& 0\\ 0& \frac{{\lambda }^2}{\kappa }& 0& \left(\frac{-1}{\kappa ·\mathrm{EI}}\right)\\ {\kappa }^2& 0& \left(\frac{-1}{\mathrm{EI}}\right)& 0\\ 0& \frac{{\kappa }^2}{\lambda }& 0& \frac{1}{\lambda \mathrm{EI}}\end{array}\right]·\left(\begin{array}{c}{y}_o\\ {\theta }_o\\ M_o\\ V_o\end{array}\right)& \left(\mathrm{e46}\right)\end{array}$

[0175] The state vector at location x=L can now be determined, using equation (e47)

[0176] $\begin{array}{cc}\left(\begin{array}{c}{y}_L\\ {\theta }_L\\ M_L\\ V_L\end{array}\right)=\left(\frac{1}{{\lambda }^2+{\kappa }^2}\right)\hspace{1em}[\begin{array}{cccc}\cosh \left(\kappa L\right)& \sinh \left(\kappa L\right)& \cos \left(\lambda L\right)& \sin \left(\mathrm{\lambda L}\right)\\ \mathrm{\kappa sinh}\left(\kappa L\right)& \mathrm{\kappa cosh}\left(\kappa L\right)& -\mathrm{\lambda sin}\left(\lambda L\right)& \mathrm{\lambda cos}\left(\lambda L\right)\\ {\kappa }^2\mathrm{EI}·\cosh \left(\kappa L\right)& {\kappa }^2\mathrm{EI}·\sinh \left(\kappa L\right)& -{\lambda }^2\mathrm{EI}·\cos \left(\lambda L\right)& -{\lambda }^2\mathrm{EI}·\sin \left(\lambda L\right)\\ -{\kappa }^3\mathrm{EI}·\sinh \left(\kappa L\right)& -{\kappa }^3\mathrm{EI}·\cosh \left(\kappa L\right)& -{\lambda }^3\mathrm{EI}·\sin \left(\lambda L\right)& {\lambda }^3\mathrm{EI}·\cos \left(\lambda L\right)\end{array}]·\hspace{1em}\left[\begin{array}{cccc}{\lambda }^2& 0& \frac{1}{\mathrm{EI}}& 0\\ 0& \frac{{\lambda }^2}{\kappa }& 0& \left(\frac{-1}{\kappa ·\mathrm{EI}}\right)\\ {\kappa }^2& 0& \left(\frac{-1}{\mathrm{EI}}\right)& 0\\ 0& \frac{{\kappa }^2}{\lambda }& 0& \frac{1}{\lambda \mathrm{EI}}\end{array}\right]·\left(\begin{array}{c}{y}_o\\ {\theta }_o\\ M_o\\ V_o\end{array}\right)& \left(e47\right)\end{array}$

[0177] The matrices may be multiplied to obtain equations (e48), for which the components of the transfer function matrix T are written individually.

[0178] $\begin{array}{cc}\left(\begin{array}{c}{y}_L\\ {\theta }_L\\ M_L\\ V_L\end{array}\right)=\left(\frac{1}{{\lambda }^2+{\kappa }^2}\right)·T·\left(\begin{array}{c}{y}_o\\ {\theta }_o\\ M_o\\ V_o\end{array}\right)T_{11}={\lambda }^2\cosh \left(\kappa L\right)+{\kappa }^2\cos \left(\lambda L\right)T_{12}=\frac{{\lambda }^2}{\kappa }\sinh \left(\kappa L\right)+\frac{{\kappa }^2}{\lambda }\sin \left(\lambda L\right)T_{13}=\frac{1}{\mathrm{EI}}\left(\cosh \left(\kappa L\right)-\cos \left(\lambda L\right)\right)T_{14}=\frac{1}{\mathrm{EI}}\left(\frac{\sin \left(\lambda L\right)}{\lambda }-\frac{\sinh \left(\kappa L\right)}{\kappa }\right)T_{21}=\kappa \lambda ·\left(\mathrm{\lambda sinh}\left(\kappa L\right)-\kappa \sin \left(\lambda L\right)\right)T_{22}={\lambda }^2\cosh \left(\kappa L\right)+{\kappa }^2\cos \left(\lambda L\right)T_{23}=\frac{1}{\mathrm{EI}}\left(\mathrm{\kappa sinh}\left(\kappa L\right)+\mathrm{\lambda sin}\left(\lambda L\right)\right)T_{24}=\frac{1}{\mathrm{EI}}\left(\cos \left(\lambda L\right)-\cosh \left(\kappa L\right)\right)T_{31}={\lambda }^2{\kappa }^2\mathrm{EI}\left(\cosh \left(\kappa L\right)-\cos \left(\lambda L\right)\right)T_{32}=\lambda \kappa ·\mathrm{EI}\left(\mathrm{\lambda sinh}\left(\kappa L\right)-\mathrm{\kappa sin}\left(\lambda L\right)\right)T_{33}={\kappa }^2\cosh \left(\kappa L\right)+{\lambda }^2\cos \left(\lambda L\right)T_{34}=\left(-\kappa \sinh \left(\kappa L\right)-\mathrm{\lambda sin}\left(\lambda L\right)\right)T_{41}=\left(-{\lambda }^2{\kappa }^2\mathrm{EI}\left(\mathrm{\kappa sinh}\left(\kappa L\right)+\mathrm{\lambda sin}\left(\lambda L\right)\right)\right)T_{42}={\lambda }^2{\kappa }^2\mathrm{EI}\left(\cos \left(\lambda L\right)-\cosh \left(\kappa L\right)\right)T_{43}={\lambda }^3\sin \left(\lambda L\right)-{\kappa }^3\sinh \left(\kappa L\right)T_{44}={\kappa }^2\cosh \left(\kappa L\right)+{\lambda }^2\cos \left(\lambda L\right)& \left(\mathrm{e48}\right)\end{array}$

[0179] In the absence of contact within the element, when all of the states are known at the first location (x=0), the states at a second location can be calculated (x=L). Just as in the above solution for the lumped parameter model, the intermediate states and matrices may be combined so that the calculation comprises a matrix relationship from one contact point to the next. However, the shear load V_o is not known because there will be a dynamic sideforce at the first contact point to match the constraint at the second contact, namely that the displacement is equal to zero for the dynamic perturbation. The four known quantities then facilitate calculation of the unknowns in the same manner as for the lumped parameter model.

[0180] Additionally and alternatively, it should be noted that the variables M and V may be normalized to new variables μ and v, respectively, by dividing by a scale factor EI_o that is characteristic of the BHA. The terms in the equations above, and the corresponding static and dynamic computations, may then be adjusted by the scale factor for the new state vector in the scaled variables, (y θ μ v 1)^T. In all other respects, the solution methods for the continuous mass equations are the same as those for the lumped parameter model.

Three Dimensional Frequency-Domain Models

[0181] The models and methods described above are essentially two-dimensional, considering the lateral dynamic bending vibration in a plane for the “flex” mode and the centrifugal effects in the “twirl” mode. By extending these methods to include both transverse coordinates, and by preserving the frequency-domain approach, advanced models can be developed to provide a three-dimensional representation to more accurately represent these bending and centrifugal vibrations. These revised and improved models would consider the dynamic effects of angular momentum and its effect on the BHA vibrations, including the gravitational effects. The following discussion provides an example of extending the above methods into a three-dimensional frequency-domain model. The teachings of the following examples can be adapted in a variety of ways depending on the configuration of the bottom hole assembly being considered. Additionally or alternatively, certain assumptions or conventions utilized in the exemplary methods below can be adjusted with alternative assumptions and/or conventions without departing from the scope of the present disclosure and claims.

[0182] FIGS. 3A-D provide a schematic view of a conventional bottom hole assembly 300 with drill collars and stabilizers 312. FIG. 3A illustrates a perspective view of the bottom hole assembly 300 as it may be bent during rotation; FIG. 3B illustrates a top-view looking downhole at a cross section of the bottom hole assembly 300. As illustrated in FIG. 3, the x-axis is oriented uphole, the z-axis is in the vertical plane orthogonal to x, and the y-axis forms the third orthogonal direction in a right-handed system. A short section of this bottom hole assembly rotates about the centerline of the wellbore with frequency, Ω, at a distance, r, from the axis, as best seen in FIG. 3B. To consider periodic motion, the distance r will be defined as a function of the rotation angle about the centerline. The pipe, or segment of the bottom hole assembly, is “turning to the right” looking down the wellbore, for the purposes of discussion at angular velocity ω_o, which is in the negative sense.

[0183] The bottom hole assembly section is subjected to an applied axial loading P, shear load V at one end and V+dV at the other end, and bending moments M and M+dM, respectively, as best seen in FIG. 3C. The loads applied to this element at the ends of the section arise from the connection to similar bottom hole assembly elements above and below this bottom hole assembly section. While the representation in FIG. 3C can appear complex, starting with the fundamental physics allows the scenario of FIG. 3C to be understood in its basic elements. For example, the net force Σ{right arrow over (F)} is equal to the rate of change of linear momentum {right arrow over (P)}, and the net torque Σ{right arrow over (M)} is equal to the rate of change of angular momentum {right arrow over (H)}. In equation form, these relationships can be written as equations (e49):

[0184] $\begin{array}{cc}\begin{array}{cc}\frac{d}{dt}\left(\overrightarrow{P}\right)=\sum \overrightarrow{F}& \frac{d}{dt}\left(\overrightarrow{H}\right)=\sum \overrightarrow{M}\end{array}& \left(e49\right)\end{array}$

[0185] Kinematics

[0186] With continuing reference to FIG. 3C, the center of mass of the element is located at position {right arrow over (R)}(t), which can be calculated by equation (e50), using θ=Ωt:
{right arrow over (R)}(t)=r(t)cos(Ωt){right arrow over (j)}+r(t)sin(Ωt){right arrow over (k)}  (e50)
According to the conventions of the present exemplary methods, the pipe section is revolving about the centerline of the well at a rate, Ω, and the pipe is spinning at rotary speed (−ω_o) about its own axis, relative to the uphole positive x-axis. Accordingly, the total angular velocity vector relative to the inertial reference frame may be written as equations (e51) where the angles φ and ψ represent the rotation angles about the y′ and z′ axes respectively:
{right arrow over (ω)}_I=−ω_o·{right arrow over (i)}″+Ω·{right arrow over (i)}
{right arrow over (ω)}_I=ω_o·({right arrow over (i)}+ψ·{right arrow over (j)}−φ·{right arrow over (k)})+Ω·{right arrow over (i)}
{right arrow over (ω)}_I=(−ω_o+Ω)·{right arrow over (i)}−ψω_o{right arrow over (j)}+φω_o·{right arrow over (k)}  (e51)

[0187] For the purposes of the present example, the motion is assumed to be in a plane. While small bending strains may be present, for the present illustration all of the angular velocity is assumed to be directed along the wellbore axis, such as shown in FIG. 3D. Accordingly, the kinematics simplify and we can write the angular velocity vector as equation (e52):
{right arrow over (ω)}_I=(−ω_o+Ω)·{right arrow over (i)}  (e52)

[0188] Since the allowable motions in these exemplary methods are constrained to the y-z plane and the motions can be resolved in a bottom hole assembly section body coordinate system that is rotating about the borehole centerline (but not spinning with the bottom hole assembly), the derivative operator applied to a vector {right arrow over (Q)} may be written as equation (e53):

[0189] $\begin{array}{cc}\frac{d\overrightarrow{Q}}{dt}={\left(\frac{\partial \overrightarrow{Q}}{\partial t}\right)}_{\mathrm{rel}}-\Omega Q_z\overrightarrow{j}+\Omega Q_y\overrightarrow{k}& \left(e53\right)\end{array}$
Additionally, the position vector may be written as equation (e54):

[0190] $\begin{array}{cc}\overrightarrow{R}=\left(\begin{array}{c}0\\ R_y\\ R_z\end{array}\right)=\left(\begin{array}{c}0\\ r\cos \left(\theta \right)\\ r\sin \left(\theta \right)\end{array}\right)& \left(e54\right)\end{array}$

[0191] Then we may write the velocity as equations (e55):

[0192] $\begin{array}{cc}\frac{d\overrightarrow{R}}{dt}=\left(\begin{array}{c}0\\ v_y\\ v_z\end{array}\right)\mathrm{for}\begin{array}{c}{v}_y=\stackrel{.}{r}\cos \left(\theta \right)-r\Omega \sin \left(\theta \right)\\ v_z=\stackrel{.}{r}\sin \left(\theta \right)+r\Omega \cos \left(\theta \right)\end{array}& \left(\mathrm{e55}\right)\end{array}$
And the acceleration may be written as equations (e56):

[0193] $\begin{array}{cc}\frac{d^2\overrightarrow{R}}{d{t}^2}=\left(\begin{array}{c}0\\ a_y\\ a_z\end{array}\right)\mathrm{for}\begin{array}{c}{a}_y=\ddot{r}\cos \left(\theta \right)-2\stackrel{.}{r}\Omega \sin \left(\theta \right)-r{\Omega }^2\cos \left(\theta \right)\\ a_z=\ddot{r}\sin \left(\theta \right)+2\stackrel{.}{r}\Omega \cos \left(\theta \right)-r{\Omega }^2\sin \left(\theta \right)\end{array}& \left(\mathrm{e56}\right)\end{array}$

[0194] For θ=0, velocity and acceleration may be written as equations (e57):

[0195] $\begin{array}{cc}\frac{d\overrightarrow{R}}{dt}=\left(\begin{array}{c}0\\ \stackrel{.}{r}\\ r\Omega \end{array}\right)\frac{d^2\overrightarrow{R}}{d{t}^2}=\left(\begin{array}{c}0\\ \ddot{r}-r{\Omega }^2\\ 2\stackrel{.}{r}\Omega \end{array}\right)& \left(\mathrm{e57}\right)\end{array}$
And for θ=π/2, velocity and acceleration may be written as equations (e58):

[0196] $\begin{array}{cc}\frac{d\overrightarrow{R}}{dt}=\left(\begin{array}{c}0\\ -r\Omega \\ \stackrel{.}{r}\end{array}\right)\frac{d^2\overrightarrow{R}}{d{t}^2}=\left(\begin{array}{c}0\\ -2\stackrel{.}{r}\Omega \\ \ddot{r}-r{\Omega }^2\end{array}\right)& \left(\mathrm{e58}\right)\end{array}$

[0197] Linear Momentum

[0198] Using the equations derived above, the equations for linear momentum are simply written as equations (e59):

[0199] $\begin{array}{cc}\frac{d\overrightarrow{P}}{dt}=\frac{d}{dt}\left(m\frac{d\overrightarrow{R}}{dt}\right)=m\frac{d^2\overrightarrow{R}}{d{t}^2}=m\left(\begin{array}{c}0\\ a_y\\ a_z\end{array}\right)=\sum \overrightarrow{F}& \left(\mathrm{e59}\right)\end{array}$

[0200] Angular Momentum

[0201] The total angular momentum may also be developed using the principles and methods described above. For example, the total angular momentum is the sum of the angular momentum of the center of mass about the borehole axis plus a term used to represent the spinning bottom hole assembly section, which can be written as equation (e60):
{right arrow over (H)}=(mr²Ω−I_xω_o)·{right arrow over (i)}  (e60)

[0202] There are no components in the orthogonal directions, and the moments of inertia are defined along a principal component system of the bottom hole assembly section, so the relation for the derivative results in the equation (e61):

[0203] $\begin{array}{cc}\frac{d}{dt}\left(\overrightarrow{H}\right)=2\mathrm{mr}\stackrel{.}{r}\Omega ·\overrightarrow{i}=\sum \overrightarrow{M}& \left(\mathrm{e61}\right)\end{array}$
Without including the terms involving the inclination angles of the element, there are no additional terms along the y and z directions to be considered in this exemplary method and model. Other models within the scope of the present methods may relax the assumptions regarding the ranges of allowable motions.

[0204] Formulation of the Differential Equations of Motion

[0205] With the foregoing equations and discussion as a backdrop, the differential equations of motion for a representative bottom hole assembly section can be formulated. FIG. 4 provides a schematic illustration of a beam element model 400 of a section of bottom hole assembly with uniform properties (density ρ and cross-section A) in a wellbore inclined at an angle φ. The gravitational force component per unit length in the z-direction is therefore (−ρAg sin(φ)). The element is oriented at an angle θ relative to the wellbore axis. The axial load is applied perpendicular to the cross-section of the element, and the shear load is parallel to the end face of the element. A differential increment in force or moment is assumed at the right-hand end of the beam. A force and moment balance on this element will provide the differential equation of motion for the beam in the z-direction.

[0206] Assuming small angles and neglecting higher order terms, and allowing for a force imbalance in the z-direction, the acceleration of the beam in the z-direction can be written as equation (e62), again using the term W=(−ρAg):
−Pθ−V+W sin(φ)dx+(V+dV)+(P+dP)(θ+dθ)=(ρAdx)a_z  (e62)
Simplifying, the acceleration of the beam may be written as equation (e63):

[0207] $\begin{array}{cc}W\sin \left(\phi \right)+\frac{dV}{dx}+P\frac{d\theta }{dx}=\left(\rho A\right)a_z& \left(\mathrm{e63}\right)\end{array}$
Continuing with the development of the equations of motion in the z-direction, it can be assumed that the moments balance to zero for the modeled element and that neglecting higher order terms is appropriate, the moments in the z-direction may be written as equations (e64):

[0208] $\begin{array}{cc}-M+\left(M+\mathrm{dM}\right)+W\sin \left(\phi \right)\mathrm{dx}\frac{\mathrm{dx}}{2}+\left(V+\mathrm{dV}\right)\mathrm{dx}=0\frac{dM}{dx}+V=0& \left(\mathrm{e64}\right)\end{array}$

[0209] The moment can be related to the deformation of the element and the EI product, as seen in equations (e65) and (e66):

[0210] $\begin{array}{cc}M=\mathrm{EI}\frac{d^{2}z}{d{x}^2}& \left(\mathrm{e65}\right)\\ \mathrm{EI}\frac{d^{4}z}{d{x}^4}+\frac{dV}{dx}=0& \left(\mathrm{e66}\right)\end{array}$

[0211] Combining the moment balance with the force balance in the z-direction, the motion in the z-direction may be written as equation (e67):

[0212] $\begin{array}{cc}\mathrm{EI}\frac{d^{4}z}{d{x}^4}-P\frac{d^{2}z}{d{x}^2}-W\sin \left(\phi \right)+\left(\rho A\right)a_z=0& \left(\mathrm{e67}\right)\end{array}$

[0213] The equations of motion in the y-dimension may be similarly derived. However, because there is no gravitational load in the transverse y-dimension, the corresponding equation is written as equation (e68):

[0214] $\begin{array}{cc}\mathrm{EI}\frac{d^{4}y}{d{x}^4}-P\frac{d^{2}y}{d{x}^2}+\left(\rho A\right)a_y=0& \left(\mathrm{e68}\right)\end{array}$
It should be noted that (e56) provides relationships for a_y and a_z in terms of r, θ, and Ω to be used with (e67) and (e68).

[0215] Solving the Differential Equations of Motion

[0216] The above discussion provides the differential equations of motion for the three-dimensional frequency-domain models. As one example, the exemplary equations above for the z-axis are identified as inhomogeneous differential equations because of the presence of the gravitational term. The inhomogeneous differential equation can be solved by combining the solution for the homogeneous case plus a term for the particular solution to reflect the gravity effect. The above discussion combining the static and dynamic solutions for distributed mass BHA vibration model provides a general solution for the inhomogeneous differential equation that can be represented by equation (e69), which is analogous to equation (e44) from above.
z(x)=ax²+c₁ cos h(κx)+c₂ sin h(κx)+c₃ cos(λx)+c₄ sin(λx)  (e69)
Where, as before, the terms κ and λ are defined by equations (e70) and are the same for both the z-dimension and the y-dimension.

[0217] $\begin{array}{cc}{\kappa }^2={\left[\frac{P^2}{4{\left(\mathrm{EI}\right)}^2}+\frac{W{\omega }_n^2}{g·\mathrm{EI}}\right]}^{1/2}+\frac{P}{2\mathrm{EI}}{\lambda }^2={\left[\frac{P^2}{4{\left(\mathrm{EI}\right)}^2}+\frac{W{\omega }_n^2}{g·\mathrm{EI}}\right]}^{1/2}-\frac{P}{2\mathrm{EI}}& \left(\mathrm{e70}\right)\end{array}$

[0218] As before, the derivatives can be identified relative to the system state variables,

[0219] $\begin{array}{cc}\frac{dz}{dx}={\theta }_z{M}_z=\mathrm{EI}\frac{d^{2}z}{d{x}^2}{V}_z=-\mathrm{EI}\frac{d^{3}z}{d{x}^3}& \left(\mathrm{e71}\right)\end{array}$

[0220] Within an element, the position x=0 can be chosen at one face, and the opposite face is then at x=L. At the origin, the cosine functions have unity value, and the sine functions are zero. In normalized state vector form, the equations for z and its derivatives may be written in matrix form as equation (e72):

[0221] $\begin{array}{cc}\left(\begin{array}{c}{z}_o\\ {\theta }_o\\ M_o\\ V_o\\ 1\end{array}\right)=\left[\begin{array}{ccccc}1& 0& 1& 0& 0\\ 0& \kappa & 0& \lambda & 0\\ {\kappa }^2\mathrm{EI}& 0& \left(-{\lambda }^2\right)\mathrm{EI}& 0& 2a·\mathrm{EI}\\ 0& \left(-{\kappa }^3\right)\mathrm{EI}& 0& {\lambda }^3\mathrm{EI}& 0\\ 0& 0& 0& 0& 1\end{array}\right]·\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ 1\end{array}\right)& \left(\mathrm{e72}\right)\end{array}$
Where the entity a is defined as follows in equation (e73):