HARP软件说明书之CHARM3D

CHARM3D USER’S MANUALThe program CHARM3D is written with standard FROTRAN 77 language and is expected to be machine independent.1. INTRODUCTIONCHARM3D is a finite element program for the coupled dynamic analysis of moored multi-column compliant offshore structures, such as TLPs SPARs, and FPSOs. The program performs coupled dynamic analyses for the platform-mooring (riser) system both in time domain and frequency domain. In the time domain analysis, various nonlinearities, such as the drag force on the mooring lines and columns of the platform, the large (translational) motion of the platform, the free surface effects, and the geometric nonlinearity of the mooring system, are included in the time marching scheme. In the frequency domain analysis, the platform--mooring system is assumed undergoing small motions around the mean position so that a linear analysis can be applied. Meanwhile, the drag forces on the mooring lines and columns of the platform are linearized using a statistical linearization technique (Rodenbush et al, 1986). Except the drag force, all the hydrodynamic forces on the platform, which include first-order and second-order sum- and difference-frequency potential forces, are computed from a hydrodynamics program WINTCOL and the computed results are imported into CHARM3D for the ensuing analysis. Therefore, CHARM3D serves as a “post-processor” of the program WINTCOL. (now a post-processor of WAMIT as well)In CHARM3D, the environmental inputs include regular and irregular waves, current, and steady wind force on the platform. The program does not compute dynamic wind forces. Nevertheless, it can read, if necessary, the dynamic-wind-force time series from a user pre-generated data file and include it in the time marching scheme.In CHARM3D, the elastic rod model is used for modeling the mooring system. This model is ideal for small strain, large displacement structural analysis of slender members such as tether, riser and catenary mooring lines. A single global coordinate system is used in the finite element formulation of the rod model. Therefore, the model is simpler and more efficient than other conventional nonlinear models, such as the updated Lagrangian beam model. Detailed theory and finite element modeling of the rod can be found in Nordgren (1974) and Garrette (1982). In CHARM3D, the model was improved to include the stretch of the element under the axial tension. Special situations such as water surface piercing and sea bottom the mooring lines (risers) are also considered in the program.In CHARM3D, the platform but FPSO is assumed to be a rigid body undergoing small (rotational) motions in waves, wind and current. Under this assumption, time domain hydrodynamic forces on the platform, including second-order sum- and difference-frequency forces, can be generated using the information from frequency-domain computations (Ran & Kim, 1995,1996). The viscous effect on the platform is also includedin the form of Morison’s drag formula if the platform is formed by relatively slender elements.For a turret moored FPSO analysis, the platform undergoes a large-angle rotation about the vertical axis (yaw rotation). Therefore, wave, wind, and current forces are computed based on large rotational motion assumptions. For wind and current forces, the modeling follows the method suggested by Oil Company International Marine Forum (OCIMF). For wave loads, the forces are computed by dividing yaw angle into several intervals. In each interval, linear and second order wave loads and motions are computed.In CHARM3D, the connections between the platform and mooring are modeled by linear translational/rotational springs, and linear translational dampers. By controlling the stiffness coefficients of these springs, various types of connections including hinged and fixed connections can be modeled. The connections between the foundations and the moorings can be similarly modeled in addition to using conventional boundary conditions.In time domain, the program CHARM3D can perform coupled quasi-static analysis without include the mooing line dynamics; and perform uncoupled motion analyses for the platform with user defined linear restoring-stiffness matrix modeling of the mooring system. It can also perform uncoupled motion analyses for the mooring (riser) system with user defined platform motions.The program also has an option that allows users to define the wave components for random waves. This option is useful if the time history of an incident wave is known. Otherwise, the phases of component waves are generated from a random-number generator. For example, the user can use FFT to obtain the wave components (amplitudes and phases) from the wave elevation time history of a model test and input them into CHARM3D so that a direct comparison of motion (or force) time series, instead of statistics, can be made between the simulation and the model test for the same wave.Since the CHARM3D is an independent program which reads the hydrodynamic information of the platform from a data file, it can also be used with other hydrodynamic programs such as WAMIT. In that case, users need to make sure that the hydrodynamic information follows the format of CHARM3D when preparing the input data files. A program named “” is developed to convert the WAMIT output data into CHARM3D format.2. PROGRAM COMPILING AND RUNNINGThe program consists of two files: CHARM3D.FOR and WPARA.BLK. The first file contains all the source codes of the program, and the second file specifies the parameters that control the dimensions of the variables used in the program. These parameters are explained in the file WPARA.BLK and users can easily make changes if necessary. It is recommended that users always check the file WPARA.BLK when they prepare new input data files to ensure that the dimensions of the variables are large enough for the type of analysis they want to do. Users need to re-compile the program CHARM3D.FOR whenever a change is made in the file WPARA.BLK. In others words, if CHARM3D aborts when running on a new model for the first time, the chance is that the WPARA.BLK needs to be changed, followed by a recompilation.As an alternative, the user may compile CHARM3D.FOR into different versions of executables based on different parameters within WPARA.BLK, so that he/she can choose to run an executable for certain cases without recompilation. If executables of CHARM3d are placed in a directory included in the system search path, they can be used in the same way as system commands.The program has been tested under UNIX and VAX operating systems on a work station, and WINDOWS NT/98 on a PC (with Watcom, Lahey and DEC FORTRAN compilers). Depending on the variable size, which is defined by the parameter in WPARA.BLK, the memory size required by the program varies significantly. It is recommended to have over 64 MB of memory size on PC to run the program efficiently. Under Unix environment, CHARM3D.FOR should be compiled with an explicit "-static" option, if it is not the default option. And a moderate, rather than agressive, optimization level may enhance the performance of the executable.3. INPUT DATA FILEThere are at least two input data files to be prepared by users. The first one, CHARM3D.IN, contains the information about mooring system (element definition, material property, platform-mooring coupling, mooring line boundary conditions) and run-time command. The second file, CHARM3D.WV, contains hydrodynamic data from WINTCOL (added mass, radiation damping, first-order forces and second-order sum- and difference-frequency forces, and wave drift damping), wave and current parameters, and information about platforms. Figure 1 shows the architecture of how to run CHARM3D.EXE and the associated input/output files. Different from other floating structures, turret-moored FPSO requires additional input files. These additional input files are also described in figure 1.The input and output data are dimensional and standard unit systems are used (meters, kg, Newton and second for SI system; feet, slug, lb and second for English system). The name of each input data in the input file follows the convention of FORTRAN Language, i.e., the variables starting from I to N are integers, otherwise, the data are real numbers. The data given in {} are needed only for three-dimensional analyses.There are two coordinate systems in the program: global coordinate system and rigid-body (platform) coordinate system. In the 3D analysis, the x and y axis define a horizontal plane and z axis is positive upward. In the 2D analysis, x-z coordinate system is used with z axis positive upward. The origins of both the global and rigid-body coordinate systems are located on the mean water surface with the corresponding axis parallel to each other. All the information related to the mooring system is defined with respect to the global coordinate system. On the other hand, hydrodynamic quantities for the platform are given with respect to the rigid-body system. In the following sections, all the input and output data are given in global coordinate except those specified otherwise.Figure 1. Flowchart of the main program CHARM3D.EXE3.1 CHARM3D.INThe input data in CHARM3D.IN can be divided into four groups. The first group of data contains general information needed in the program and is listed as followsNDIM, NLEG, NMAT, GRAVITY, RHOWMITER, NSTEP, DT, NINTBX0(1), BX0(2), {BX0(3)}IFLAGFPSOBETWV, NHEADHEAD1, HEAD2, HEAD3, HEAD4,…., HEAD NHEADIn this group, NDIM defines the dimension of the analysis. It is 2 for two-dimensional analysis and 3 for three-dimensional analysis. NLEG is the total number of tethers, risers and mooring lines in the analysis. In the following, the tethers, risers, and mooring lines are referred as “legs”. NMAT is the total number of sets of material properties in a mooring system. GRAVITY is the gravitational acceleration. GRAVITY=9.8065 in SI unit system, and GRAVITY=32.17 in English unit system. RHOW is the water density. MITER is the maximum number of iterations specified in the static analysis. In the static analysis, the nonlinear equation is solved by using Newton iteration scheme. If the tolerance that is specified later in the file is not satisfied after MITER iterations, the program will stop and print a massage on screen. NSTEP is the total number of steps and DT is the time interval in the time-marching scheme. The accuracy of the time domain simulation is primarily controlled by the time interval, namely, the value of DT. In general, a stronger nonlinear problem requires smaller time interval. For instance, in the coupled analysis of TLP or SPAR, it is recommended that the DT have a value that allows at least 30 time steps in the smallest period of interest. The NINT is a parameter that controls the output in the time domain simulation. The program will output the results for every NINT steps in the time domain simulation. This allows the user to control the amount of output data without changing the time interval DT. For example, with time interval DT=0.005 seconds, the program will output the results for every 0.005 second interval if NINT=1, and will output the results for every 0.1 second interval if NINT=20. Those two variables are used for the mooring lines or flexible risers with part of the line lying on the sea floor. BX0 defines the position (in global coordinate system) of the origin of the local rigid-body coordinate system. IFLAGFPSO is the integer used to specify whether a turret-moored FPSO is evaluated.IFLAGFPSO= 0: Evaluate a platform except FPSOIFLAGFPSO= 1: Evaluate a FPSO platformBETWV is wave headings to be analyzed.NHEAD is equal to the number of selected angle between wave heading and rigid-body coordinate system.These angles are related to how the hydrodynamic forces are computed using a hydrodynamics program such as WAMIT.HEAD1, HEAD2, HEAD3, HEAD4,…., HEAD NHEAD is the angles (degree) between wave heading and rigid-body coordinate system.HEAD1 must be defined to the smallest angle while HEAD NHEAD is the largest one. A convention about these angles is illustrated in figure 2.Figure 2. Definition of HEAD’s AnglesThe second group of data in CHARM3D.IN contains all the material properties of the legs. In CHARM3D, each leg is described by a finite number of rod elements. The material properties (stretch and bending stiffness, etc.) are constant in each element. However, they may vary from element to element so that a non-uniform leg such as a mixed wire-chain mooring line can be modeled. There will be NMAT sets of material properties with NMAT being the total number of sets as given in the first group. The input format for material properties is as follows:GAE1, GEI1, GRHOL1, GRHOA1, GCI1, GCD1, GAS1, BKTEN1...GAE n, GEI n, GRHOL n, GRHOA n, GCI n, GCD n, GAS n, BKTEN n...GAE NMAT,GEI NMAT,GRHOL NMAT,GRHOA NMAT,GCI NMAT,GCD NMAT,GAS NMAT, BKTEN NMATwhere GAE is the axial stiffness (Young’s modulus cross sectional area of the leg), and GEI is the bending stiffness (Young’s modulus moment of inertia of cross section). GRHOL is the mass per unit length of the element, and GRHOA is the displaced mass per unit length if the element is in water (GRHOA=0 if the element is in air). GCI is the coefficient of inertial force, i.e., the inertia force per unit length at unit acceleration. GCD is the coefficient of drag force, i.e., the drag force per unit length at unit relative velocity. GAS is the cross sectional area of the element. GAS is only used in calculating axial stretch of the element in water. If the actual cross section is used, the stretch is computed using actual tension in the element. If GAS=0, then the stretch is computed using effective tension (effective tension = actual tension + hydrostatic pressure GAS). BKTEN is the linebreaking tension which will be used as a criteria to determine if the line is broken during the simulation.The third group of input data in CHARM3D.IN contains information about legs. This group of data must be prepared for every leg. In other words, one should repeat this group of input data for NLEG times, where NLEG is the total number of legs as defined in group one. The first two parameters in this group areNELEM,IFLAG1NELEM is the total number of elements in a leg, and IFLAG1 is a flag defining two ways of inputing the elements’ geometry . If IFLAG=0, the elements are defined in the following format:R(1)1,R(2)1,{R(3)1},RP(1)1,RP(2)1,{RP(3)1},TZER1...R(1)n,R(2)n,{R(3)n},RP(1)n, RP(2)n, {RP(3)n}, TZER n...R(1)NELEM+1,R(2)NELEM+1,{R(3)NELEM+1},RP(1)NELEM+1,RP(2)NELEM+1,{RP(3)NELEM+1}, TZER NELEM+1GLEN1, IOPTN1, MAT1...GLEN n, IOPTN n, MAT n...GLEN NELEM, IOPTN NELEM, MAT NELEMwhere R n is the nodal coordinate of the n-th node in the leg, and RP n is the unit tangential vector (directional cosine) of the leg at that node. The tangential vector of a node always points to the direction of a higher numbered node nearby. There are a total of NELEM+1 nodes in the leg. The nodal numbering should be in a consecutive manner from one end to the other end of the leg. For a coupled analysis, the first node in a leg should be at the end which is connected to the bottom boundary (sea floor), and the last node should be at the other end of the leg which is connected to the rigid body (platform). For an uncoupled analysis, the numbering of the nodes can start from either ends of the leg. TZER n is the pretension in the leg at n-th node. This input affects axial deformation of the element in the analysis. For example, in a TLP analysis, if the stretch of the tether caused by the net buoyancy force is already included in the initial element length (GLEN), one needs to set the TZER as the actual static pretension so that any axial deformation in the calculation is caused only by dynamic tension. GLEN n is the length of the n-th element that is between node n and n+1. IOTPN is an option for each element. It is 1 if that element touches or is likely to touch the foundation (sea floor), and 2 if that element is at or near free surface and likely to pierce the water surface. Otherwise, it is 0. When IOTPN=1, the program will check whether that element touches the foundation and whether the stiffness matrix of the element needs to be modified to include the contribution of the foundation. When IOPTN=2, the program will check the position of free surface with respect to the element and change the added mass and hydrodynamic force computation accordingly. MAT n is theset number of the material properties for the n-th element. All sets of the material properties are previously defined in the second group.The above input format for the element geometry is suitable for the legs whose initial positions are known, such as the vertical tether of a TLP. However, it is inconvenient to use this format for a catenary mooring line since the initial nodal positions of the lines are not easy to find. In that case, user can use a more convenient format by setting IFLAG=1.If IFLAG1=1, the following format for the definition of element position should be used: PA(1), PA(2), {PA(3)}PB(1), PB(2), {PB(3)}GLEN1, IOPTN1, MAT1,TZER1...GLEN n, IOPTN n, MAT n, TZER n...GLEN NELEM, IOPTN NELEM, MAT NELEM, TZER NELEM+1where the PA is the nodal position of the first node, and PB is the nodal position of the last node. For coupled analysis, the first node in a leg should be at one end which is connected to the bottom boundary (sea floor), and the last node should be at the other end of leg which is connected to the rigid body (platform). The definitions for GLEN, IOPTN, MAT and TZER are the same as those for IFLAG=0.Following the element geometry input, leg boundary conditions and other information for each leg are to be given as follows:NBVPIBVP(1)1,IBVP(2)1,IBVP(3)1,BVP1...IBVP(1)NBVP,IBVP(2)NBVP,IBVP(3)NBVP,BVP NBVPNBVSIBVS(1)1,IBVS(2)1,IBVS(3)1,BVS1...IBVS(1)NBVS,IBVS(2)NBVS,IBVS(3)NBVS,BVS NBVSNEND1GSL(1),GSL(2),{GSL(3)},GSR,GE(1),GE(2),{GE(3),}GX(1),GX(2),{GX(3)}GDLNEND2GSL(1),GSL(2),{GSL(3)},GSR,GE(1),GE(2),{GE(3),}GX(1),GX(2),{GX(3)}GDLNGAPNOD 1, NSPG 1GGSL 1 (1), GGSL 1 (2), GGSL 1 (3)GGAP 1 (1), GGAP 1 (2), GGAP 1 (3)GGX 1 (1), GGX 1 (2), GGX 1 (3)GGCS 1,GGCD 1,GPLOAD 1 …….NOD n , NSPG nGGSL n (1), GGSL n (2), GGSL n (3)GGAP n (1), GGAP n (2), GGAP n (3)GGX n (1), GGX n (2), GGX n (3)GGCS n ,GGCD n ,GPLOAD n …….NOD NGAP , NSPG NGAPGGSL NGAP (1), GGSL NGAP (2), GGSL NGAP (3)GGAP NGAP (1), GGAP NGAP (2), GGAP NGAP (3)GGX NGAP (1), GGX NGAP (2), GGX NGAP (3)GGCS NGAP ,GGCD NGAP ,GPLOAD NGAPNBUOYIBUOY 1, BUOY 1 (1), BUOY 1 (2), BUOY 1 (3)…….IBUOY n , BUOY n (1), BUOY n (2), BUOY n (3)…….IBUOY NBUOY , BUOY NBUOY (1), BUOY NBUOY (2), BUOY NBUOY (3)NSPRINGISPRING 1BGSL 1 (1), BGSL 1 (2), BGSL 1 (3), BGSR 1,BGE 1 (1), BGE 1 (2), BGE 1 (3), BGX 1 (1), BGX 1 (2), BGX 1 (3) BGDL 1……..ISPRING nBGSL n (1), BGSL n (2), BGSL n (3), BGSR n ,BGE n (1), BGE n (2), BGE n (3), BGX n (1), BGX n (2), BGX n (3) BGDL n…….ISPRING NSPRINGBGSL NSPRING (1), BGSL NSPRING (2), BGSL NSPRING (3), BGSR NSPRING , BGE NSPRING (1), BGE NSPRING (2), BGE NSPRING (3), BGX NSPRING (1), BGX NSPRING (2), BGX NSPRING (3)BGDL NSPRINGDBOTTM, CBOTTM, FBOTTMBKTIME, IFLAG6NOUTNEOUT(1), ..., NEOUT(NOUT)where NBVP is the total number of essential (displacement) boundary conditions in the leg. If it is not 0, for each essential boundary condition, users need to define: IBVP(1)--the nodal number where boundary condition applies; IBVP(2)--the direction of the boundary condition (1 in x, 2 in z direction for 2D analysis; 1 in x, 2 in y and 3 in z direction in 3D analysis); IBVP(3)--is 1 if the boundary condition is to define the position of the node and 2 if the boundary condition is to define the unit tangential vector of the node; BVP-- the value of the boundary condition. For example, an input of 8,1,1,20.0 means that the x coordinate of node 8 is 20.0 (in global coordinate system). As an another example, the input of 23,3,2 0.7 means that the z component of the unit tangential vector at node 23 is 0.7. Similarly, NBVS is the total number of natural (force) boundary conditions. If it is not 0, for each natural boundary condition, users need to define: IBVS(1)--the nodal number where boundary condition applies; IBVS(2)--the direction of the boundary condition (1 in x, 2 in z direction for 2D analysis; 1 in x, 2 in y and 3 in z direction in 3D analysis); IBVS(3)--is always 1; BVS--the value of the external force applied. For example, an input of 11,2,1,2.0e5 means that 2.0e5 (Newton or lb.) of force is applied in the y direction on node 11.For the nodes at the two ends of the leg, special boundary conditions are considered in addition to the ones described above. The ends can be connected to the boundary or rigid body (platform) by linear translational and rotational springs, and linear translational dampers. NEND1 is a flag for the first end (node 1). It is 1 if the end is connected to the boundary by springs and damper, and 0 if not. If NEND1=1, one needs to input the following information of spring and damper connection: GSL--the linear spring stiffness (unit: force/length) in each direction; GSR--the rotational spring stiffness (unit: moment/radian); GE--the unit vector defining the rotational spring reference axis. A moment will be created if there is a non-zero angle between the unit tangential vector at the end of the leg and the reference axis GE. GX is the position of the point on the boundary where the linear springs are connected; GDL is the damping coefficient (unit: force/velocity) of the damper. A damping will be created if there are relative translational velocities between the end of the leg and connection point on the boundary. NEND2 is a flag for the other end of the leg (node NELEM+1) which is connected to a rigid body (platform) or fixed boundary. It is 1 if the end is connected to a fixed boundary by springs and dampers (this condition is for uncoupled analysis only), 2 if the end is connected to a rigid body (platform), and 0 if there is no spring and damper connections at this end (this condition is for uncoupled analysis only). For option 1 or 2, the user should input the information about the spring and damper which are defined in a similar manner to those of the other end. GX will be the position of the connecting point on platform, and GE will be the unit vector fixed on the rigid body (platform) defining the rotational spring reference axis. When NEND1=0 or NEND2=0, springs and dampers are not used for the boundary condition. For instance, hinged and fixed joints can be easily modeled by using only essential boundary conditions without springs and dampers. However, for a coupled analysis, the spring and damper connections must be used between platform and legs. The user can vary the stiffness to simulate different types of connections. For instance, a large translational spring stiffness and zero rotational stiffness can simulate a hinged connection between platformand legs; Sufficiently large values for both translational and rotational stiffness can simulate a fixed connection. In these cases, it is recommended that the values of GSL and GSR should be about 1,000~10,000 times higher than the axial (stretch) stiffness (GAE) and bending stiffness (GEI) of the legs, respectively.If NEND2=3, the pneumatic tensioner, a kind of nonlinear spring, is applied at the top boundary instead of the linear spring connection. The tensioner is only in vertical (z) direction and the other directions keep the linear spring coupling. For the pneumatic tensioner, the tension at the top is calculated by the following equation()001/n T T z z =+ , wherez : stroke of the piston with upward direction positive,0z : effective length of gas in the associated accumulator,0T : Initial top tension at stroke z=0.n : gas constant.The connection is the additional one from version 3.6. The new format for the pneumatic tensioner is as follows.Input for the pneumatic hydraulic top tensioner (NEND(LEG,2)=3)GSL(1), GSL(2), GSL(3), GSRGE(1), GE(2), GE(3)GX(1), GX(2), GX(3)GDLDYMU, SZ, GASN, where GSL(1) and GSL(2) are the linear horizontal spring constant and GSL(3) is the initial top tension at zero stroke (). GSR is the rotational spring stiffness (unit:moment/radian). GE is the unit vector defining rotational spring reference axis and the pneumatic tensioner aligned axis as well. Usually, the tensioner is aligned in vertical direction, thus, GE(1)=0, GE(2)=0 and GE(3)=1. GX is the position of the point on the boundary where the pneumatic tensioner and linear horizontal springs are connected; GDL is the damping coefficient (unit: force/velocity) of the damper only for horizontal direction.A linear damping in horizontal direction will be created if there are relative translational velocities between the end of the leg and connection point on the boundary.0T In the vertical direction a kind of Coulumb friction damping mechanism is applied to model the friction of the piston. Thus the dynamic friction coefficient DYMU is input at the next line for the nonlinear damping calculation. SZ is the effective length of the gas . GASN is the gas constant (). The pneumatic tensioner is fully nonlinear model and is only applied to the time domain simulation.0z nAnother additional modules for the version 3.6 is the gap contact between TTR/Buoyancy Can and Spar guider. Without this model, the effect of the buoyancy can is only taken into account by the truncated TTR at the keel joint and the additional pitch/roll hydrostatic stiffness. However, the guide and the can may not contact when thebuoyancy can does not affect the hull motion. When they are contact, there can be a relative motion between them(SLIP), or there can be no relative motion(STICK). For the SLIP condition, the Coulumb friction between the can and the guide makes a roll as a damping force in heave direction. The real system is more complicated than the truncated TTR model. Thus we need the gap contact module which can exactly model the realistic contact and STICK/SLIP phenomena.At first, the program reads the NGAP which indicates the number of the gap contact or the pneumatic tensioner in the leg. The group of the input will be repeated NGAP times to specify the properties of each connection model.The first line of the group for each connection consists of NOD and NSPG which are, respectively, node number of the connection and the order of the spring in the direction of the contact. NSPG can be 1, 2 or 3 to choose linear, quadratic and cubic spring, respectively.GGSL is the spring constant of the gap contact. In the direction of the gap contact, the spring constant given is applied only when the guide and the can is in contact. The spring constant normal to the gap contact should be given as zero(0.0). Otherwise the program considers the connection as a simple linear spring coupling. GGAP defines the distance between the riser center line and the guider at which the riser and the guide are contacting with zero reaction force. If GGAP is set to zero, it is assumed that riser is always in contact with the guide. GGX is the gap contact point on the platform, namely the point where the guide is positioned. GGCS is the static friction coefficient which is used when the riser is stick to the guide. GGCD is the dynamic friction coefficient used when the riser is slipping contacting the guide. GPLOAD is the preload with which the guider initially grabs the buoyancy can or the riser.When there is small sub-buoys connected to a mooring line, those can be considered as a point force acting on the nodal point of the line element. The sub-buoys can be out of the water, surface piercing or fully submerged due to the mooring tension effect. NBUOY is the number of the buoy on the mooring line nodes. The NBUOY g roups of the input for the buoy properties are followed. IBUOY is the node number to which the sub-buoy belongs. BUOY(1) is the heave hydrostatic stiffness defined by (Water Plane Area)×(Water Density) × (Gravitational Acceleration). BUOY(2) is the total length of the buoy to determine the position of the buoy relative to the free surface. BUOY(3) is the initial draft of the buoy. Every time step, the program checks the vertical position of the buoy to determine if the buoy is out of the water, surface piercing or submerged. If it is out of the water, all the buoy related forces are removed. While it is surface piercing, the vertical spring whose stiffness is BUOY(1) is attached to the node point to model the hydrostatic stiffness due to buoy. If it is submerged, the constant buoyancy force due to the displacement(BUOY(1) ×BUOY(2)) of it is applied.NSPRING is the number of the spring at nodal points of the line element followed by the NSPRING groups of different spring properties. It is used when a nodal point of the line element is connected to a fixed point in space by a spring. The spring can be in x-, y- or z- direction and rotational one as well. If we give only the vertical spring fixed on the free surface with the other spring stiffnesses zero, the model can represent the surface piercing。