Strojniški vestnik
Journal of Mechanical Engineering
Strojniški vestnik - Journal of Mechanical Engineering (SV-JME)
Aim and Scope
The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis.
The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue.
The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s).
Editor in Chief
Vincenc Butala
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Technical Editor
Pika Škraba
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Founding Editor
Bojan Kraut
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office
University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www. sv-jme.eu
Print: Grafex, d.o.o., printed in 310 copies Founders and Publishers
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
University of Maribor, Faculty of Mechanical Engineering, Slovenia
Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association
President of Publishing Council
Branko Širok
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
International Editorial Board
Kamil Arslan, Karabuk University, Turkey
Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia
Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia
Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia
Anselmo Eduardo Diniz, State University of Campinas, Brazil
Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia
Imre Felde, Obuda University, Faculty of Informatics, Hungary
Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia
Imre Horvath, Delft University of Technology, The Netherlands
Aleš Hribernik, UM, Faculty of Mechanical Engineering, Slovenia
Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan
Julius Kaplunov, Brunel University, West London, UK
Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany
Jernej Klemenc, UL, Faculty of Mechanical Engineering, Slovenia
Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia
Peter Krajnik, Chalmers University of Technology, Sweden
Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia
Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia
Thomas Lubben, University of Bremen, Germany
Janez Možina, UL, Faculty of Mechanical Engineering, Slovenia
George K. Nikas, KADMOS Engineering, UK
José L. Ocana, Technical University of Madrid, Spain
Miroslav Plančak, University of Novi Sad, Serbia
Vladimir Popovič, University of Belgrade, Faculty of Mech. Eng., Serbia Franci Pušavec, UL, Faculty of Mechanical Engineering, Slovenia Bernd Sauer, University of Kaiserlautern, Germany Rudolph J. Scavuzzo, University of Akron, USA Arkady Voloshin, Lehigh University, Bethlehem, USA
Vice-President of Publishing Council
Jože Balič
University of Maribor, Faculty of Mechanical Engineering, Slovenia
Strojniški vestnik Journal of Mechanical "''i^'' Engineering
Cover:
Front cover shows a measurement and welding robotic application for a product called protector installed in cast-iron cooking plates. The robotic cell enables following of few manufacturing parameters of the protector and calculation of a welding point for each individual protector. Besides all measurement and pose transformation calibrations are performed automatically in regular intervals.
Image Courtesy:
Laboratory of robotics, Faculty of Electrical Engineering, University of Ljubljana, Slovenia
Photo: Jure Rejc
ISSN 0039-2480
General information
Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue).
Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor.
Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http://en.sv-jme.eu/.
You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content.
We would like to thank the reviewers who have taken part in the peerreview process.
© 2016 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website.
The journal is subsidized by Slovenian Research Agency.
Strojniški vestnik - Journal of Mechanical Engineering is available on http://www.sv-jme.eu, where you access also to papers' supplements, such as simulations, etc.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)12 Contents
Contents
Strojniški vestnik - Journal of Mechanical Engineering volume 62, (2016), number 12 Ljubljana, December 2016 ISSN 0039-2480
Published monthly
Papers
Jure Rejc, Marko Munih: Robust Visual Touch-Up Calibration Method in Robot Laser Spot Welding
Application 697
Weiping Wang, Bo Wang: An Energy-Saving Control Strategy with Load Sensing for Electro-Hydraulic
Servo Systems 709
Tjasa Duh Coz, Andrej Kitanovski, Alojz Poredos: Primary Energy Factor of a District Cooling System 717 Cheng Gu, Xinbo Chen: A Novel Universal Reducer Integrating a Planetary Gear Mechanism with an
RCRCR Spatial Mechanism 730
Aijun Yin, Juncheng Lu, Zongxian Dai, Jiang Li, Ouyang Qi: Isomap and Deep Belief Network-Based
Machine Health Combined Assessment Model 740
Nnamdi Onochie Chibueze, Chinwuba Victor Ossia, John Umunna Okoli: On the Fatigue of Steel
Catenary Risers 751
Yu Dai, Liping Pang, Lisong Chen, Xiang Zhu, Tao Zhang: A New Multi-Body Dynamic Model of a
Deep Ocean Mining Vehicle-Pipeline-Ship System and Simulation of Its Integrated Motion 757
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)12, 697-708 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2016.3708 Original Scientific Paper
Received for review: 2016-05-03 Received revised form: 2016-09-29 Accepted for publication: 2016-10-07
Robust Visual Touch-Up Calibration Method in Robot Laser Spot Welding Application
Jure Rejc* - Marko Munih
University of Ljubljana, Faculty of Electrical Engineering, Slovenia
The article is describing the use of visual touch-up calibration method for defining the mathematical transformations used in a shop-floor measurement and welding robot cell in the protector assembly process. The presented system is designed as a robust and cheap solution, using only the equipment needed for the production tasks in the robot cell. The main goal of the presented system is to use vision measurement system for measuring and calibration procedures and laser welding equipment to weld two protector assembly parts together where positioning tolerances are very narrow. These narrow tolerances forced us to implement auto-checking and auto-calibration procedures for all necessary mathematical aspects in the robot cell, based on the robust visual touch-up method. To demonstrate adequate solution in the measurement, calibration and also the production sequences, the graphs show production statistical results over a one year production period.
Keywords: robot welding, visual touch-up, calibration, kinematic error, transformations Highlights
• The presented system enables accurate distance measurements in the industry.
• The presented system enables robust visual touch-up robot calibration method.
• The proposed calibration procedures take into account all robot kinematic errors.
• All calibration procedures are automatic, enabling short production line dead times.
• The system uses a calibration procedure that is not well described in the literature.
0 INTRODUCTION
Robot spot welding is nowadays present in several industries all around the globe. These systems increase production efficiency [1] and increase the quality of the products. The physical burden on the human workforce is relieved as well as the stress on their health [2]. Many of these systems are present in automotive industry [3] where contact spot welding is mainly used [4]. Besides automotive industry the robotic spot welding is present also in other industries [5] and [6].
Industrial robots have high position repeatability, but have at least a grade worse absolute position accuracy [7] to [9]. The robots are mainly programmed on-line where all robot points are defined or recalculated in regard to the base coordinate system of the robot. However, improvements in technology enable off-line robot programming to be used more and more nowadays. This type of programming saves the robot points in regard to the virtual robot coordinate systems. When these points are transferred to the real system on the shop-floor, usually a point position difference is present and the literature [10] specifies this error as positional absolute error or kinematic position error. The same problem occurs when machine or robot vision [11] systems are used
[12], where points in camera coordinate system need to be transformed into the robot base coordinate system. When transforming coordinates from vision system to the robot system, usually the ideal robot kinematic model is used, but real kinematic parameters differ. For this reason an absolute calibration procedure is a must to accurately position the robot on proper position defined by the vision system.
To reduce or eliminate the absolute error, manual calibration of the robot system is usually used. But this conventional approach requires a large amount of calibration points, which results in a long calibration time and is therefore not suitable for shop-floor production. In the field of robotic automatic absolute error calibration procedures the reader can find several approaches using 1D and 2D vision calibration systems [13] to [16], calibration with laser trackers [17], image comparison [18], visual touch-up [19] and hybrid sensors using Kalman filters [20].
Most of the presented work for kinematic calibration of the robot system in the previous paragraph was tested in laboratory and used expensive dedicated measurement equipment or additional equipment needed to be installed that limits robot working space. These academic approaches have a large influence on calibration methods development, but are usually not implemented in shop-floor
*Corr. Author's Address: University of Ljubljana, Faculty of Electrical Engineering, Tržaška 25, 1000 Ljubljana, Slovenia, jure.rejc@robo.fe.uni-lj.si
697
production facilities [21]. All presented drawbacks forced us to develop the robot automated visual inspection (AVI) and measurement system for measuring and laser welding cell in a way that used previously installed equipment for measuring, welding and for all necessary calibration procedures. Among the presented approaches from the literature the simplified visual touch-up approach was implemented. The article presents the robot cell design, the automatic position error eliminating procedures and the results of installed approach. Presented is statistics for the whole year production analysis.
1 THE ROBOT CELL
1.1 The Protector
The robot cell workpiece called protector (Fig. 1) is a control and safety element incorporated into classic cast-iron cooking plates manufactured in different diameters and nominal power. The task of the protector is to turn off the power supply of the heating winding hobs in the event of overheating when the temperature reaches 400 °C ± 50 °C.
The basic components of the protector (Figs. 1 and 2) are: ceramic housing, 1.2 mm thick bimetal with the set screw, limiter, toggle element, electrical switch and electrical contacts for connecting wires. When the temperature of the cooking plate is rising, the bimetal bends in the protector and exerts force via the set screw on the limiter. When the pressure on the limiter is high enough, it triggers the toggle element, which represents half of the electrical switch.
Fig. 1. The protector and its assembly parts
The protector manufacturer was forced to change its design in 2013 for two reasons. The set screw was previously fixed to the bimetal with special glue paint. This solution was practical, but it sometimes happened that the set screw was not fixed enough and the protector switch-off temperature moved outside
the tolerances. Also, this special paint was expensive, which called for a cheaper and more reliable solution. The protector redesign declared that the set screw is bonded with the bimetal by laser welding these two assembly parts together.
Fig. 2. The field of view of the visual inspection system
1.2 The Robot Workspace and Attached Equipment
The robot cell (Fig. 3) is installed in the fourth of the five stages of the rotary table, where the assembly process of the protectors is finished. In this stage the task is to inspect two dimensions called A and B in the protector and to determine the intersection position of the bimetal and the set screw where laser welding of these two parts must be performed. The position of the set screw is set in the previous stage of the rotary table and is not important for the article.
Fig. 3. Robot cell as one part of the five stage rotary table
The selected robot is Epson G6 650 with 4 DOF and superb repeatability specifications: ±15 ^m for the first and second horizontal axis together, ±10 ^m for the vertical axis and ±0.005 ° for rotational axis. On the end of the robot two independent systems are installed, both necessary for all the production tasks in the fourth stage of the rotary table. The first system is a video camera inspection and measurement system. It consists of a video camera and appropriate optics. Installed camera is a type DMK 41AG02, monochrome with resolution of 1280x960 pixels, produced by The Imaging Source. We have chosen optics from the Keyence company, type CA-LM0510. It is specified as macro lens with C-mount connection. The field of view of this video system is approximately 6 mm x 4 mm, marked with a dotted square in Fig. 1. The second attached system is a laser welding system with appropriate optics. The type of laser welding system is TruPulse 44, manufactured by the company Trumpf with a wavelength of 1064 nm and average power of 40 W. The laser optics is BEO D35 with a focus distance of f = 100 mm and 1 mm laser beam spot. This parameter combined with the bimetal width of 1.2 mm defined the laser welding point (WP) tolerance to ±0.2 mm. The laser welding system is equipped with an additional video camera of the same type as the one in the measurement system and is used for calibration procedures of both systems. The camera attached on the welding system share the same optics, making the laser beam and the video camera visual path coaxial.
The production process of the fourth stage of the rotational table starts as follows: six protectors, set in a cluster, are rotated into the robot working space at once. After the cluster is positioned, additional mechanism positions the dedicated LED illumination for all protectors in a cluster. Then the robot positions the video measurement system over the first protector in a cluster, the dedicated image acquisition and image processing software captures the image, which is then processed during the motion to the next protector. When the last protector image in a cluster is processed, the robot moves in the opposite direction from protector to protector and positions the welding optics over the welding point according to the information from the measurement system. The laser welding is not performed if there is an error in image analysis or the measured dimensions A or B are not in defined tolerances.
In Fig. 3 two very important parts of the whole system are also seen. The measurement system first needs to be checked and calibrated to ensure accuracy. In our opinion the best object to perform
the measurement calibration procedures is a precisely known object that is also measured in the robot cell. That is why three calibration protectors are set in the robot working space in a special chamber protected from laser welding dust as much as possible. The height of the optics in regard to these protectors was set by the same robot vertical Z axis distance as by protectors fixed in a cluster rotated by the rotary table. This is possible because both production and calibration protector clusters are physically set to the same height in the production line. This simple approach minimizes the influence of camera intrinsic parameters error and also the optics distortion error. Both dimensions A and B on all three calibration protectors were previously measured with the certified profile projector measuring system, type Mitutoyo PV500. Each protector has different dimensions A and B. The captured image with all important parameters can be seen in Fig. 2, except that the welding point is not defined during the calibration of the measurement video system. The first two calibration protectors are used to gain the transformation information used to recalculate distances from pixels to millimeters and the last calibration protector is used to check the measurement accuracy. The checking is performed every 10,000 pieces and if the accuracy is inside predefined tolerances of ±0.1 mm for both measured distances A and B then the production line continues. Otherwise the robot moves the measurement system over the first two calibration protectors and a new transformation function is calculated. Then the accuracy checking is repeated on the third protector. If the measured values are still outside tolerances, the production line is stopped with error message and an operator must check the situation.
Fig. 4. Calibration coordinate system fixed in robot working space
The second important system in the robot working space in Fig. 3 is very simple, yet very efficient as we will show in the article. The 3D model of this assembly is shown in Fig. 4. The assembly consists
of an L stand that has a 5 mm hole on the bottom side. It is meant for a 5 mm LED. For better contrast we have chosen a green LED. On the top side of the stand two small holes in the 5 mm LED area are drilled. The bigger hole has a diameter of 0.8 mm and the smaller one a diameter of 0.5 mm. These two holes or dots represent one axis of fixed coordinate system (Fig. 5) in robot working space. The larger dot represents the origin (O) of coordinate system and the smaller a dot on the X axis of the coordinate system. Both dots are separated by 3 mm.
Fig. 5. Threshold image of two LED dots with drawn coordinate system
2 THE VISUAL TOUCH-UP METHOD
The visual touch-up method is a non-contact version of a standard touch calibration method and can be used in many robot calibration areas. The non-contact method can be connected with contact method if the term of virtual pin is introduced (Fig. 6). Virtual pin is a virtual connection from the robot to the target position. The literature is very poor in the field of visual touch-up method used for robot calibration purposes and only Watanabe et al. [19] published a contemporary research article in this field, which was used as the basis for our approach. Watanabe et al. used a single camera attached on the robot end-effector. The calibration target object is a perfect circle with its center point drawn in the robot working space. The size of the circle is predefined and is used to define geometric relations, where the center point is the target point. The authors state that the drawbacks of this approach are unidentified camera-intrinsic parameters and the distortion of the lens that can both affect calculations.
Visual touch-up method can be used for purposes of robot new tool calibration, robot absolute accuracy calibration and also for calibration of several robots carried vision systems as in our case. The non-contact method can use several sensors for calibration procedures: from laser distance sensors based on triangulation [22] and conoscopic holography [23], inductive or capacitive sensors and especially video
cameras as Watanabe et al. is presenting. In our case we have chosen the video camera approach, because both on robot attached systems are vision based and also the planar robot movement simplifies the calibration procedures approach (Fig. 6). The reference objects of our visual touch-up approach are small round green dots presented in Figs. 4 and 5.
Fig. 6. Visual touch-up system represented on 3D model
3 IMPLEMENTED VISUAL TOUCH-UP METHOD CALIBRATIONS
As described in the previous sections, the measurement vision system is checked for accuracy and calibrated by using the calibration protectors. With the implementation of this system, the welding point (PX, PY) is defined in camera image coordinate system [C, independent in regard to any other coordinate systems of the robot cell (Fig. 2). But to be able to transform the welding point (PX, PY) from [C in to the robot reference coordinate system [R and to position the welding optics to the proper position (WPX, WPY) several coordinate systems need to be defined automatically via the visual touch-up method.
These calibration procedures calibrate the system only in X and Y axis, where the Z axis is fixed. The focal distance of the welding laser and its attached camera is 11 cm from the welding optics to the observed or welded object, defined with manual calibration stick, provided by the laser manufacturer. The same analogy valid also for the measurement camera, where the focus distance is also near 11 cm, making the stand-off distance also fixed. In other case the captured images are blurred.
3.1 Robot Self-Calibration of the New Welding Laser Optics Tool
The robot can position the welding laser optics to the calculated coordinates WPX and WPY if the welding laser optics frame [L] regarding to the robot default end [e] E is defined as a transformation TL or as a predefined robot new tool (Fig. 7), also called tool center point (TCP). In production facilities this is frequently done by manually defining a new tool attached to the end of the robot with a special Epson wizard for manually defining new TCP's named as ToolN (e.g. Tooll). This procedure requires a certain amount of time of at least a few minutes to be finished by the operator and is therefore too inconvenient for a quick recalibration and totally improper for high volume production line.
Fig. 7. The positioning of the welding laser optics according to the WPX, WPY welding point
In our system the robot makes fast movements that can shift attached equipment, or the operator accidentally hits the laser optics during maintenance of the robot cell. For these reasons an automatic tool calibration procedure for defining the center of the welding laser as a new tool, called Tooll, was developed.
Fig. 8. Transformations for new TCP calculation
For the automatic tool calibration procedure a bigger LED dot from Fig. 5 is used. The calibration procedure requires positioning of the laser welding optics |L in the center of a reference point (bigger
dot), marked as a filled circle in Fig. 8, in two different robot configurations and can be described with Eq. (1) where N is 1 and 2. To make all further figures transparent, only the x axis of the coordinate frame is marked and the z axis points out of the plane. The y axis is set respectively to the right-hand coordinate system.
TR _ TR Tl _ TEN '
(1)
In this calibration procedure the robot first moves the laser optics to the predefined position (EXN, EYN, aN) in robot reference frame [R, saved in the previous calibration procedure, where the center of the laser optics and the center of the reference point should align. The decision whether the welding laser optics center is aligned with the center of the reference point is made by the dedicated software by implementing a circular Hough transform [24] to [26] on the captured laser optics video camera image. The task of the Hough transform method is to search for objects of different shapes (lines, circles, ellipse) in an image by a voting approach in parameter space. Within this space the objects are gained as local maximum in an accumulator space.
Unlike in Watanabe et al. [19] where the error between the reference point and the captured image point is calculated and used in further calculations, we implemented a simple step position controller to reduce the position error inside the predefined tolerance area of 0.05 mm if the movement is necessary. The movement of the robot is in steps of 0.015 mm in both planar axes. The tolerance area can be specified in millimeters because the width of the bigger dot is known and the result of the Hough transform is the radius of the circle in pixels. At this point the current robot TCP position [E] is saved as a new point (EXN, EYN, aN) for the next calibration attempt. With this information a new transformation TR (Eq. (2); N = 1) is set.
TeRn = Rot (z,aN ) • Trans (Em, Em, 0).
(2)
In order to calculate the transformation TL between the end of the robot (E) and welding laser optics (L) a second transformation is needed (Eq. (1); N = 2). It defines the new configuration (Eq. (2); N = 2) of the robot, pointing with the center of the welding laser optics in the same, bigger LED reference point. The procedure is the same as described before, only the robot initial pose is different EXN, EYN, aN where N equals 2.
From the matrix in Eq. (3) only position coordinates ELX and ELY are needed and can be
expressed as Eq. (4). Coordinates ELX and ELY define the new TCP (saved as Tooll) representing the welding laser optics center (L) relative to the end of the robot (E).
Ter = Rot (z,v) ■ Trans (EX, EY ,0).
(7)
Tle = Trans (ELX, ELY, 0),
ELX =(EX1 - EX 2 )-(cosa1 - cosa2) + + (EY1 - EY2 )-(sina1 - sina2),
ELy = (Ey1 -EY2)-(cosat -cosa2)-
-(EX l - EX 2 )'(sin«l - sin
4.2 Welding Point Calculation
(3)
(4)
Fig. 9. Transformations for welding point calculation
Fig. 9 shows the homogenous transformation relations that are important for calculation of the welding point in the reference robot frame [R. The welding point, marked as P, is determined in the measurement camera frame [C]. This transformation can be written as a homogenous transformation matrix (Eq. (5)) marked as with the same orientation as the measurement camera frame [cl.
Tt>
(5)
= Trans (PX, PY ,0).
The camera [C] is physically fixed in regard to the robot end [E] described by the transformation T^ (Eq. (6)). This transformation is defined with the distance (ECX, ECy) from robot end frame [E] to the measuring camera frame [c] and with a rotation angle 9 around the Z axis.
TCe = Rot (z,3)-Trans (ECX, ECY ,0).
(6)
The robot end frame |Ej pose in robot reference frame [R can be written as homogenous transformation T^ (Eq. (7)) with parameters (EX, Ey,
0 0, if ^<0'
(2)
where is the flow gain coefficient of the servo valve, kx is a positive constant, ps is the supply pressure of the pump, pt is the tank pressure, p is the oil density, p1 and p2 are the chamber pressures of the cylinder, Q1 is the supply flow rate to the forward chamber and Q2 is the return flow rate from the return chamber, ud is the input signal of PDV.
The pressure dynamics of the actuator can be written as:
[À = h1[Q1 - A1Xp - CtPl + CiP2] 1 p2 = h2[-Q2 + A2xp - CtP2 + CiPl ]
(3)
where xp is the piston position, Ct=Ci+Ce is the total leakage coefficient [17], Ce is the external leakage coefficient, Ci is the internal leakage coefficient, h1 = fîe / Vj , h2 = fie / V2 , Vj, V2 are the control volumes of the two chambers, respectively. fie is the effective bulk modulus of the system.
The force balance equation of the system is expressed as:
mxp = Pl Al - Pl 4 - Bc Xp - FL - Pf*
(4)
where m is the equivalent mass of the load, xp is the displacement of the cylinder rod, Bc is the coefficient of viscous damping, FL is the external load force, Ff is the Coulomb friction force.
Define the system state variables as Eq. (5):
[ X1, x2, X3, X4] _ [ Xp, pi, p2]
(5)
In order to make the system fall into the strict feedback form, the state variables are reconstructed as Eq. (6), where a = iA1 .
[xPx2,x3]T =[x,xv,pi -ap2f
(6)
Considering the uncertainties and disturbance in the model, the entire system can be expressed in a state space form as Eq. (7):
■ - a- _ B _ FL_Fl d
—2 — —3 —2 d^i
—3 — f1—2 /2—3 + /3 —4 + f4Uà ++ d2
(7)
y — — 1
/ = h a + h2 Aa
/2 = hCt + h2Cia f = hC + h2Cta f4 = h kq kx g + h2kq kx ^
(8)
wheref ~f4 are shown in Eq. (8) and dx, d2 represent the lumped modelling error including unmodelled dynamics and external disturbance. It can be seen that the system is highly nonlinear because it suffers from the lumped modelling error and parametric uncertainties.
The control task is summarized as follows: given the desired motion trajectory xd , the objective is to synthesize a control input of the PDV such that the output xj tracks xd as closely as possible. Meanwhile, The PRV is utilized to regulate the supply pressure ps in order to reduce the energy consumption.
2 ENERGY-SAVING STRATEGY
For a fixed displacement pump, there always exists excess pressure and overflow loss. The VSPC can reduce the excess pressure loss of the system. However, the overflow loss is not considered in these studies. The load-sensing pump can adjust the flow and pressure per the demand. Thus, the overflow loss can be eliminated by employing a load-sensing pump. However, the pressure margin of the load-sensing pump, usually fixed approximately j4 bar to 30 bar across the valve, is not optimized in studies involving load-sensing pumps. Thus, if the VSPC could be applied to the load-sensing pump, the system efficiency would be significantly improved.
2.1 Configuration of the Load-Sensing Structure
The load-sensing structure of the proposed system is shown in Fig. 2. The structure comprises the load-sensing pump, the PRV, and the throttle valve. If the pump pressure is larger than the sum of the preset pressure margin pm and feedback pressure pLs , the cylinder will be driven to reduce the angle of the swash plate 6 until the pressure is balanced on both
sides and vice versa. The pressure margin pm is preset by the spring inside the load-sensing pump. However, this margin can be indirectly regulated by adjusting the feedback pressure pLs , which is controlled by the PRV. The function of the throttle valve is only to establish the pressure difference between the pump and the PRV. The power consumed in the throttle valve is neglected because the orifice of the throttle valve can be adjusted to a small value and the pressure difference across the valve is equal to pm when the system is stable.
PC PRVThrottle Valve
Load sensing pump
Fig. 2. Schematic diagram of load sensing
2.2 Variable Supply Pressure Control
For a certain demand flow, the smallest possible pressure drop will occur if the valve is fully opened. However, in that case, the PDV cannot be adjusted in a certain range to compensate for the tracking error. It can be concluded that if the opening of the PDV remains at a relatively high level during the forward and backward motion while the pump only provides the demand pressure, the pressure loss across the valve will be minimized. The VSPC is introduced as follows.
For the forward motion (Xd > 0), the desired inlet pressure drop across the PDV Apj = ps -pj can be calculated as follows:
A1 Xd = kq kxaUmax4| "Aft,
Api =
PAi2 xd2
2£q2ka2uL
(9)
(10)
where Mmax is the maximum input voltage of PDV during the forward motion, a is the desired normalized input signal of PDV determined by the users/designers to indicate the importance of energy saving. The value of a can adjust the desired spool position of the PDV, and a higher value of a can improve the energy efficiency.
The pressure drop between pump and valve can be obtained as:
Ap2 = X
d 2
(11)
where X is the layer resistance coefficient, l is the pipe length, d is the pipe diameter.
Because of uncertainties in the plant, the design of the desired pump pressure must allow for a margin of safety. Thus, op is added to the value of desired pump pressure [3]. From Eqs. (9), (10) and (11), the desired pump pressure, the indirect feedback pressure and the input signal of PRV can be written as:
Psdi = Pi par- + ^P2 + aiPi' (12) Ik^a wmax d 2
PlS1 = = Psdl - Pm ^ Uvl = Psi1' Pm , (13)
where krv is the PRV gain, urvl is the control input of the PRV, psd1 is the desired pump pressure in the forward motion, pLs1 is the indirect feedback pressure governed by the PRV in the forward motion. Because the supply pressure always exceeds the cracking pressure, the dynamic of the PRV can be simplified as Eq. (13) [3].
As with the case of forward motion, the desired pump pressure, and control voltage of the PRV when xd < 0 are expressed as:
Psd2 = P2 + PA2f\ P2 + (14)
2kk a w„„ d 2
Pls2 = ^K = Psd2 - Pm ^ "rv2 = ^ ^ ■ (15)
Remark 1: It should be noted that the pump pressure will remain as the preset pressure margin pm when the desired pump pressure is less than the pm. Furthermore, although the pump pressure is not zero when xd = 0 , the input signal of PDV varies slightly around zero. Consequently, the pump only provides a little flow mainly for compensating for the internal leakage. Therefore, the energy consumption, in this case, is also very small.
e n , tanh(—)e2 + m ' e,
ftX2 /2*3 f~hX4 + fhUi + ^2 -*3d)
- e2 (-02x2-0l)-p,0t0i - p,0,0,. (28)
The actual control is designed as:
Ud = f ( f-X2 + f2X3 — f3X4 + X3d — k3e3 — J4
e A -S2 tanh(—)e3--- e2),
(29)
where k is a positive constant. The adaptation law is chosen as:
(30)
In case the parameter estimate di varies significantly, which may result in a very large input voltage, 0, is updated using the following projection type adaption law [16],
Projs (• )
0 if e i >6^ and • > 0 0 if ei < 6wini and • < 0 i = 1,2. (31) • otherwise
Substituting the Eqs. (29) and (30) into Eq. (28) yields:
V3 = -ke - k2e2 - k3e3 + \e2 -
- 8X tanh(—)e2A2e3 - S2 tanh(-^)e3 < 0. (32)
Remark2: Because the proposed control strategy is a set point strategy, it is necessary to verify the stability of the system. However, the stability of the internal dynamics of the system is difficult to ensure.
A method in [3], which utilized a new state variable and zero dynamics to validate the stability, can also be used in this study. The extended state variable x5 and the dynamic of the pump pressure, which are deduced by Eqs. (12) and (14), respectively, are as follows
X5 = hJp(PS - Pi) - h2)jpp(P 2 - Pt) X1 ^ 0
X5 = hjp(Pi - Pt) - Ps - P2) X1 < 0
, (33)
P = p2 + 72
pA2 xv xv +x-Pi d 2
k2k2a2u2 q x max
pA2 xp xp d2
k2k2a2u2 q x max
. (34)
Casel. Zero output: In this case, the output of the system should be kept at zero. This means that xd = xd = xd = xd = xl = xl = x = 'x[ = 0. Thus, the control input approaches zero in this case according to Eq. (29). Accordingly, the flows of the two chambers are also near zero. Neglecting the leakage, we can get P1 = P2 = 0 from Eq. (3). Thus, the chamber pressure will be regulated at constant values. Consequently, the pump pressure is also a constant value and ps = 0 according to Eqs. (12), (14), and (34). As a result, x5 will tend to be a constant value and x5 = 0.
Case2. Constant output: In this case, xd = x ^ 0 , xd = xd =xd = xi = xi =x1 = 0 . Based on Eqs. (29) and (34), the control input and ps will approach zero again, and the rest of the analysis will be the same as Case1, which means that in a set point strategy, the internal dynamics of the system will be stable.
4 EXPERIMENTAL RESULTS
The setup is presented in Fig. 3. The position of the cylinder is measured by a displacement sensor (LS 628C). The pressures in the two cylinder chambers and supply pressure are measured by pressure sensors (HDP702). A dSpace1104 data acquisition data board is used to acquire the feedback signals from sensors and to generate control signals. The variable pump is a Danfoss load-sensing pump. The parameters of the system are shown in Table 1.
Fig. 4 depicts the harmonic reference tracking test for the proposed method by selecting a as 0.86. It can be observed that the actual trajectory follows the desired position well. The pump pressure and chamber pressures are shown in Fig. 5, in which the maximum pressure margin is 7.5 bar. It is less than the
2
preset pressure margin 14 bar in the pump. Moreover, with the increase in the flow, the pressure margin is boosted, and the decrease of the pressure margin is associated with decreases in flow. A higher pressure margin always appears at the middle stroke of the cylinder in both forward and backward motion, which could verify the effectiveness of the energy-saving method.
overshoots due to the low-pressure difference and the chattering problem. As can be seen from Fig. 7, a significant tracking error of 3 mm occurs in the peaks of the reference signal. Moreover, the average tracking errors of the proposed method and the sliding mode controller are 0.75 mm and 0.9 mm, respectively. Hence, the proposed method has a better tracking
Table 1. Parameters of the system
Total leakage coefficient 3x10"11 Oil bulk modulus 1x109 Pa
Piston/rod diameter 40/30 mm equivalent mass of the load 15 kg
Oil density 890 kg/m3 Viscous damping 800 N/m
kqk-x 8.8x10"7 0.005
1 0.2
02 0.05 ki 100
k2 200 k3 135
Fig. 3. Schematic of the equipment
3 6 9
Time (s)
Fig. 4. Tracking performance
The parameter estimates are shown in Fig. 6, and the tracking errors of three different methods are shown in Fig. 7. It can be observed that the proposed method exhibits a better tracking performance than PID control with a maximum tracking error of 1.2 mm. If sliding mode control is used to achieve the same level of tracking performance, high feedback gain must be employed, which may result in
6 9
Time (s)
Fig. 5. Pump pressure and cylinder chamber pressures
e ¡il
3i 210 -1-2-3
Time (s>
Fig. 6. Estimation results
Proposed Control Strategy Siliding Mode control PID
3 6 9 12
Time (s)
Fig. 7. Tracking errors (a = 0.86)
- Proposed control strategy
- Sliding mode control
3 6 9 12
Time (s)
Fig. 8. Tracking errors (a = 0.7)
performance than the sliding mode control regarding the average error and maximum error. To further verify the tracking performance of the proposed method, another experiment was carried out by selecting the a as 0.7. As can be seen from Fig. 8, the proposed method demonstrates a better tracking performance compared with the sliding mode controller. It can also be found that although the tracking error of the proposed method is slightly smaller when a is 0.7, it is not obvious. Meanwhile, the energy-saving performance will degrade by choosing a smaller value of a. The suitable value of a is determined by trial-and-error and eventually chosen as 0.86 by taking both energy saving and position tracking performance into account.
To evaluate the energy-saving performance of the proposed method, it is compared with a fixed displacement system. The system pressure of the fixed displacement system is assumed to be 50 bar, and the flow provided by the pump is assumed to be 16 l/min, which is slightly higher than the maximum demand flow. The energy-saving performance of the proposed system was also evaluated by comparing the energy efficiency with that of the fixed displacement system with VSPC, which is similar to the existing VSPC system. Neglecting the energy loss across the throttle valve, the supply energy of the three methods are shown in Fig. 9. The energy is calculated by PsQL , in which QL is the flow rate of the pump. It can be observed from Fig. 9 that the supply energy of a fixed displacement system is 20 kJ, whereas the supply energies are 14.8 kJ and 7.5 kJ in the fixed displacement system with VSPC and the proposed system, respectively. Thus, it can be concluded that the proposed method is capable of saving 62.5 % energy in comparison with a fixed displacement system and approximately 49 % energy in comparison with a fixed displacement system with VSPC.
- Proposed system
- Fixed displacement system with VSPC Fixed displacement system
6 9
Time (s)
Fig. 9. Supply energy of the system
To further verify the proposed method, multistep reference tests were conducted by selecting a as 0.82. The results are depicted by Figs. 10 to 12. It can
be observed from Fig. 10 that the proposed method succeeded in realizing position tracking and the dynamic behaviour of the presented system is good. The pump pressure and chamber pressures are shown in Fig. 11. It can be observed that the maximum pressure margin is 10 bar and the pump pressure remains at 14 bar when the desired pump pressure is lower than the preset pressure margin. The supply energy during the whole motion cycle is shown in Fig. 12. The system pressure of the fixed displacement system is assumed to be 45 bar. It is shown in Fig. 12 that the supply energy of the fixed displacement system is 28 kJ. Compared with the fixed displacement system, the fixed displacement system with VSPC saved approximately 62 % energy, and the proposed system saved approximately 90 % energy.
0.16-, 0.140.120.100.080.060.040.02
-Reference
— Actual
V 1 1 \ I'l fl 1 1
3 6 9
Time (s)
Fig. 10. Tracking performance
3 6 9 12 15
Time (s)
Fig. 11. Pump pressure and cylinder chamber pressures
- Proposed system Fixed displacement system with VSPC
- Fixed displacement system
1)
6 9
Time (s)
Fig. 12. Supply energy of the system
5 CONCLUSIONS
A VSPC based load-sensing structure is presented in order to reduce the overflow loss and excess
pressure loss across the PDV. Thus, the flow and pressure of the EHSS can both meet the control demand. Additionally, an adaptive backstepping sliding mode control is introduced for position tracking.
2) Compared with the fixed displacement system, the proposed method is capable of saving 62.5 % and 90 % of supply energy in harmonic reference tests and multi-step reference tests, respectively. The results also demonstrate that the presented method has a better energy-saving performance in comparison with [3] and [4] because the overflow loss can be reduced in this study.
3) The proposed method exhibits a good tracking capacity with a maximum tracking error of 1.2 mm in the harmonic reference tests. Furthermore, the tracking performance and dynamic behaviour for the multi-step reference tests are also good. This method could also be combined with the
independent metering to save even more energy in mobile machinery and other applications because the load-sensing pump presented here is widely used in engineering machinery.
6 ACKNOWLEDGEMENTS
The project is funded in part by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD), the 2016 annual general university graduate research and innovation program of Jiangsu Province, China (KYLX16_0525). We also give our great thanks to Keivan Baghestan, Amirkabir University of Technology, Tehran, Iran, for his valuable suggestions.
7 REFERENCES
[1] Lu, X.C., Chen, Q.B., Zhang, Z.J. (2014). The electric vehicle routing optimizing algorithm and the charging stations' layout analysis in Beijing. International Journal of Simulation Modelling, no. 13, no. 1, p. 116-127, D0l:10.2507/ IJSIMM13(1)C04.
[2] Casoli, P., Pompini, N., Ricco, L. (2015). Simulation of an excavator hydraulic system using nonlinear mathematical models. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 10, p. 583-593, D0I:10.5545/sv-jme.2015.2570.
[3] Baghestan, K., Rezaei, S.M., Talebi, H.A., Zareinejad, M. (2015). An energy-saving nonlinear position control strategy for electro-hydraulic servo systems. ISA Transactions, vol. 59, p. 268-279, D0I:10.1016/j.isatra.2015.10.012.
[4] Tivay, A., Zareinejad, M., Rezaei, S. M., Baghestan, K. (2014). A switched energy saving position controller for variable-pressure electro-hydraulic servo systems. ISA Transactions, vol. 53, p. 1297-1306, D0I:10.1016/j.isatra.2014.04.010.
[5] Guo, K., Wei, J., Fang, J., Feng, R., Wang, X. (2015). Position tracking control of electro-hydraulic single-rod actuator based on an extended disturbance observer. Mechatronics, vol. 27, no. 4, p. 47-56, D0I:10.1016/j.mechatronics.2015.02.003.
[6] Mohanty, A., Yao, B. (2011). Indirect adaptive robust control of hydraulic manipulators with accurate parameter estimates. IEEE Transactions on Control Systems Technology, vol. 19, no. 3. p. 567-575, D0I:10.1109/TCST.2010.2048569.
[7] Mohanty, A., Yao, B. (2011). Integrated direct/indirect adaptive robust control of hydraulic manipulators with valve deadband. IEEE/ASME Transactions on Mechatronics, vol. 16, no. 4. p. 707-715, D0I:10.1109/TMECH.2010.2051037.
[8] Micheal, J., Rahmat, M.F., Abdul, N., Lai, W.K. (2013). Feed forward linear quadratic controller design for an industrial electro hydraulic actuator system with servo valve. International Journal on Smart Sensing and Intelligent Systems, vol. 6, no. 1, p. 154-170.
[9] Niksefat, N., Sepehri, N. (2001). Designing robust force control of hydraulic actuators despite system and environmental uncertainties. IEEE Control Systems, vol. 21, no. 2, p. 66-77, D0I:10.1109/37.918266.
[10] Gao, B., Shao, J., Yang, X. (2014). A compound control strategy combining velocity compensation with ADRC of electro-hydraulic position servo control system. ISA Transactions, vol. 53, no. 6, p. 1910-1918, D0I:10.1016/j.isatra.2014.06.011.
[11] Chiang, M.H., Yeh, Y.P., Yang, F.L., Chen, Y.N. (2005). Integrated control of clamping force and energy-saving in hydraulic injection moulding machines using decoupling fuzzy sliding-mode control. The International Journal of Advanced Manufacturing Technology, vol. 27, no. 1, p. 53-62, D0I:10.1007/s00170-004-2138-z.
[12] Cho, S.H., Noskievič, P. (2012). Position tracking control with load-sensing for energy-saving valve-controlled cylinder system. Journal of Mechanical Science and Technology, vol. 26, no. 2, p. 617-625, D0I:10.1007/s12206-011-1032-5.
[13] Lovrec, D., Kastrevc, M., Ulaga, S. (2009). Electro-hydraulic load sensing with a speed-controlled hydraulic supply system on forming-machines. The International Journal of Advanced Manufacturing Technology, vol. 41, no. 11, p. 1066-1075, D0I:10.1007/s00170-008-1553-y.
[14] Axin, M., Eriksson, B., Krus, P. (2014). Flow versus pressure control of pumps in mobile hydraulic systems. Journal of Systems and Control Engineering, vol. 228, no. 4, p. 245-256, D0I:10.1177/0959651813512820.
[15] Minav, T., Hanninen, H., Sinkkonen, A., Laurila, L., Pyrhonen, J. (2014). Electric or hydraulic energy recovery systems in a reach truck - a comparison. Strojniški vestnik - Journal of Mechanical Engineering, vol. 60, no. 4, p. 232-240. D0I:10.5545/sv-jme.2013.1581.
[16] Lu, L., Yao, B. (2014). Energy saving control of a hydraulic manipulator using five cartridge valves and one accumulator. IEEE Transactions on Industrial Electronics, vol. 61, no. 12, p. 7046-7054, D0I:10.1109/TIE.2014.2314054.
[17] Kovari, A. (2015). Effect of leakage in electrohydraulic servo systems based on complex nonlinear mathematical model and experimental results. Acta Polytechnica Hungarica, vol. 12, no. 3, p. 129-146, D0I:10.12700/APH.12.3.2015.3.8.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)12, 717-729 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2016.3777 Original Scientific Paper
Received for review: 2016-06-06 Received revised form: 2016-10-14 Accepted for publication: 2016-10-18
Primary Energy Factor of a District Cooling System
Tjasa Duh Coz* - Andrej Kitanovski - Alojz Poredos
University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
The primary energy efficiency for various energy-related processes can be calculated using the primary energy factor (PEF). In this paper, the PEFs of district cooling systems (PEFDC) for different types of cold production are derived. These concern cold production with an absorption chiller driven by different available sources and cold production with a compressor chiller driven by different types of engines and related energy sources. Based on the fundamental definition of the PEF, a mathematical model for calculating the PEFDC for different types of cold production was developed. The results in this study reveal that the PEFDC can be significantly improved in the case of combined cooling and power generation. The PEFDC in the case of combined cooling and power generation is lower than when cooling with electrically driven compressor chillers when the energy efficiency of the electricity generation in thermal power plant (qj is low or the PEF of the electricity (PEFe) is high. In cold production technologies where coal is used as the primary energy source more primary energy is consumed compared to other primary energy sources (i.e. natural gas, waste heat, etc.). Keywords: primary energy factor, combined generation, district cooling, chiller
Highlights
• Primary energy factors (PEFs) appear in many legislations and research reports.
• PEF can be used to determine the primary energy consumption of cold production.
• This work provides mathematical model for calculation of PEFs for different cooling technologies.
• Moreover this work provides results of analysis of PEFs of district cooling systems.
• The PEFs of district cooling systems can be used for the selection of the most primary energy efficient cold production technology.
0 INTRODUCTION
In the past decade, the EU has been taking a more active role in the field of improving energy efficiency, reducing energy consumption and exploiting renewable energy sources. In order to define the actual primary energy efficiency of various energy-related processes, the primary energy factor (PEF) can be used as a tool. The PEF enables a comparison between the input primary energy to the system and the energy delivered to the consumer. Its evaluation involves the energy required for the extracting, processing, storing and transporting to a power plant, energy conversion, transmission, distribution and the losses associated with these processes. The primary energy in this particular case includes the energy contained in the raw fuels as well as other forms of energy received as the input to the energy-supply system. It covers both renewable and non-renewable energy sources and the definition here is in accordance with EN 153164-5 [1], which states that '... waste heat, surplus heat and regenerative heat sources are included with the appropriate primary energy factors.'
A set of directives has been approved in order to reduce the energy consumption, increase the efficiency and exploit renewable energy sources. The PEF has been used as a significant metric in order to calculate the actual primary energy efficiency of different
processes in several legislations, i.e., in Directive 2012/27/EU on energy efficiency [2], Directive 2010/31/EU on the energy performance of buildings [3] and Directive 2009/28/EC on the promotion of the use of energy from renewable sources [4]. The essential requirements from those directives appear in the relevant European standards in order to implement the calculation of the energy efficiency and the PEF. Standard EN 15203/15315:2006 [5] presents the different values of the PEF provided by the Swiss Federal Institute of Technology. In the standards EN 15316-4-5:2006 [1] and EN 15603:2008 [6] the different values of the PEFs for various types of electricity generation are recommended. All the standards listed above use fixed values of the PEFs that are commonly not updated. However, PEFs are variable, because they depend on the specific mix of primary energy sources and the efficiency of the processes of generation, storage and transportation.
Beretta et al. [7] presented a method that provides a dynamic calculation of the PEF depending on the variation of time and geographical location. Laverge an Janssens [8] attempted to use empirical data to calculate the PEFs for a set of countries in a specific year. However, his study did not follow a complete and replicable methodology and it did not include the losses associated with generation, storage and transportation. Wilby et al. [9] claimed that fixed
*Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Askerceva c. 6, Ljubljana, Slovenia, tjasa.coz@fs.uni-lj.si
717
values of PEFs are not capable of representing the evolution of the energy system. He proposed an application to calculate dynamic PEFs based on empirical data. Among the 14 countries that were included in the study, 9 countries had real PEFs of electricity below the European average of 2.5 (Finland 1.32, Norway 1.39, Denmark 1.7, etc.) and the rest of the countries had values above the European average (Poland 2.91, France 3.4, etc.). Schicktanz et al. [10], Fumo and Chamra [11] and Mago and Chamra [12] used the PEF to calculate the primary energy consumption and energy costs of a combined heating, cooling and power system (trigeneration). Different PEFs from different types of district heating and cooling systems are presented by Wirgentius
[13]. The renewable part of the primary energy was excluded from the study. Dalla Rosa and Christensen
[14] presented the primary energy performance of a network that was designed for low-temperature operation. The calculations hypothesized PEFs of 2.5 for electricity and 0.8 for district heating. Prandin [15] used the PEFs to calculate the exergy consumption of a space heating system. Zyczynska [16] explored the methodology for determining the PEF of a specific urban heating system. The methodology was based on real measurements. A PEF of 0.5 (heat from cogeneration) was obtained in this study. Jungbauer et al. [17] presented measurements of the PEF of a district cooling system in Barcelona (0.6), Copenhagen (0.93), Helsinki (0.18), Gothenburg (0.44) and Vienna (0.62). Riepl et al. [18] used PEFs to calculate the primary energy savings using solar thermal cooling and heating plants with absorption chillers compared to electric vapour-compressor chiller systems.
The aim of this study is a development of a general mathematical model for the calculation of the PEFdc. This mathematical model can be used as a tool in a feasibility study of a new district cooling system in order to select the cooling technology with the lowest PEFDC. In this paper several different types of cold production are included in the calculation of the PEFdc. This study concerns cold production with an absorption chiller driven by the heat from different available sources as well as cold production with a compressor chiller driven by different types of engines and related energy sources. A study of PEFDC in such a wide range of district cooling technologies has not yet been found in the literature.
1 DEFINITION OF THE PEF OF A DISTRIC COOLING SYSTEM
The primary energy efficiency of a district cooling system can be defined by evaluating the PEFDC,
which is the ratio between the primary energy input QP and the cooling energy at the primary side of the substation QDC. The cooling energy QDC in the Eq. (1) presents the sum of the cooling energy of all consumers connected to the district cooling system. The basic definition of the PEF for a district cooling system can be shown using the following expression:
PEFdc = QpiQd
(1)
In the case of a trigeneration system only part of the fuel input is used for the cold production. The rest of it is used for the production of heat and electricity. Therefore, in the case of a district cooling system related to a trigeneration plant, the expression for PEF in Eq. (1) can be rearranged as [19]:
Z Qf ■ PEF F
pefdc = -
Z Qd
• PEFd +Z Qdh ,, • PEFa _i__
Z QDC Jk
(2)
In the following text, the PEFDC for different types of cold production and different types of fuel input are determined. In this case the cooling-energy consumption of an absorption chiller is defined as:
QDC ,abs = QH ■ COPabs '"HdN.
(3)
For the compressor chiller, the cooling-energy consumption is defined as:
Qd
' COPcom 'nDN.
(4)
Fig. 1 illustrates the different types of cold production that are analyzed in this paper. Cold can be produced either by a compressor or by absorption chillers. Different types of absorption chillers can be used: a single-effect absorption chiller (SEAC) driven by hot water, a double-effect absorption chiller (DEAC) driven by steam and a direct-fired absorption chiller (DFAC) driven by natural gas. Waste heat (WH) from an industrial process can also be used to drive an absorption chiller. Compressor chillers are divided into two groups: compressor chillers driven by an electric motor (CC) and compressor chillers driven by an internal combustion engine (ICE). In this study, the primary energy consumption of the cooling towers is not included in the calculation of the PEFDC.
The technical parameters of the chillers are presented in Table 2.
el
Fig. 1. The different types of cold production included in this study
1.1 Cold Production with an Absorption Chiller Driven by the Heat from a Boiler (CB)
The heat from the combustion of fuel in a boiler drives the absorption chiller. A single- or double-effect absorption chiller can be used to produce cold. The primary energy consumption is:
qp,cb = ' peff .
(5)
The PEF of a district cooling system for this particular case can be defined by using Eqs. (3) and (5) in Eq. (1):
( 1 A
PEFn,
1
COPabs 'VON ncB
PEFb
(6)
Coal and natural gas are considered as fossil fuels in this study. According to [1], the PEFcoa[ = 1.3, and for natural gas, the PEFNG = 1.1 (see also Table 1).
1.2 Cold Production with a Direct-Fired Absorption Chiller Driven by Natural Gas (DFACNG)
A direct-fired absorption chiller is driven by the heat from the combustion of natural gas. According to Makita [20], the COP of a direct-fired absorption chiller is defined as:
COP
DFAC
= QJQl
(7)
The primary energy consumption of a direct-fired absorption chiller is calculated as:
qp,ng qlhv ' pefng ■
(8)
Using Eqs. (1), (7) and (8), the PEF of a district cooling system can be defined as:
( 1 A
PEFn,
1
COPDFAC '^DN
pefm,
(9)
The PEFNG is given in Section 1.1 (see also Table
1).
1.3 Cold Production with an Absorption Chiller Driven by Waste Heat from an Industrial Process (WH)
In the case of a district cooling system consisting of an absorption chiller driven by the waste heat from an industrial process, the primary energy consumption is calculated as follows:
Q
'P ,WH QH ' PEFWH.
(10)
The heat delivered to the absorption chiller can be effluent, hot water or steam. Using Eqs. (1), (3) and (10), the PEF of the waste heat from the industrial process is:
( 1 A
PEFn,
1
COP„b
■Von
PEF,,r
(11)
The waste heat includes the heat from municipal incineration and industrial surplus heat. The PEF of the waste heat is defined according to Refs. [21] and [22] as PEFwh=0.05 (see also Table 1). Using waste heat
to drive an absorption chiller avoids the use of fossil fuels and make use of the energy flows that otherwise would be lost. Hence, the value of the PEF is almost 0. When calculating the PEFDC that uses waste heat to drive an absorption chiller, it is necessary to take into account that heat with different parameters (see Table 2) is used.
1.4 Cold Production with a Compressor Chiller Driven by the Electricity from the Grid (CCmix)
Different energy sources can be used to generate electricity. This is then further used to drive the motor of the compressor chiller.
The primary energy consumption of an electrically driven compressor chiller is:
Qp = We,. PEFeI,mic. (12)
Using Eqs. (4) and (12), Eq. (1) can be rewritten
as:
PEFn,
1
COP
■ PEF,
(13)
com l DN y
The PEF of the electricity from the grid is defined according to Ref. [14] as PEFe,mix=2.5. This value represents the European average. Each country has the right to set a different value for its PEFe, mix, providing its choice is adequately justified (see also Section 0).
1.5 Cold Production with a Compressor Chiller Driven by the Electricity Generated in a Thermal Power Plant (CCTPP)
produced in a thermal power plant) is calculated using Eqs. (1), (4) and (14):
( 1 A
PEF,
1
COPcom -VdN nd
PEF
(15)
The PEFf are given in Section 1.1 (see also Table
1).
1.6 Cold Production with a Compressor Chiller Driven by an Internal Combustion Engine (CC,ce)
This kind of system consists of a compressor chiller that is driven by the mechanical energy produced with an internal combustion engine (ICE). According to Yingjian [23], the COP of an engine-driven compressor chiller is defined as:
COP rCE
com,ICE
= QcWm
(16)
The primary energy consumption of an engine-driven compressor chiller is:
QP
■ PEFf .
(17)
Using Eq. (4), (16) and (17), Eq. (1) can be rewritten as:
( , A
PEFn,
1
COP„.
' '^DN
PEF
(18)
The ICE can be driven by diesel, petrol, kerosene or natural gas; the values for the PEF of the fuel are defined according to [1] and [24] (see also Table 1).
Since the electricity from thermal power plants (TPP) represents only a part of the electricity mix, the electricity produced in the thermal power plant is discussed separately.
In this subsection a compressor chiller driven by the electricity from a thermal power plant (TPP) is studied. Two different scenarios are discussed: a compressor chiller driven by the electricity generated in a generator, connected to a steam turbine (combustion of coal), and a compressor chiller driven by electricity, where the generator is connected to a combined cycle of gas and steam turbine (combustion of natural gas).
The primary energy consumption of the electrically driven compressor chiller is calculated as:
Q =-e-. PEF
TPP 1 -^-i F
Vel
(14)
The PEF of a district cooling system with an electrically driven compressor chiller (electricity
1.7 Cold Production in Combined Cooling and Power Generation (CCP)
A trigeneration system (combined cooling, heating and power generation) is a combination of a cogeneration plant and an absorption chiller (see Fig. 1). In such a system where the cooling energy is produced with absorption chillers driven by heat, the steam should be extracted from an extraction-condensing turbine using the parameters that are required for the chiller's normal operation (see the parameters in Table 2). Consequently, the generation of electrical energy is reduced, and more primary energy is consumed in order to produce the same amount of electrical energy as in the case of a condensing turbine regime. In this subchapter the combined cooling and power generation is considered. It is assumed that all the available heat from a thermal power plant is used to produce cold. Therefore, the district heat consumption Qdh in Eq. (2) is equal to 0.
To calculate the PEFDC based on combined cooling and power generation (CCP) the Eq. (2) can be rearranged to:
PEFn,
MQ, ■ PEF,
cop„
1
■Vdn (n
(19)
where AQP presents the difference between the energy consumption of fuel, when both, electricity and cold are produced (i.e. CCP) and when only electrical energy is produced in thermal power plant (TPP) (see the Eq. (20)). In both cases the amount of produced electrical energy remains the same (Wei = Wel,TTP = Wel,ccP).
( , , \
AQf = Qf ,CCP - Qf
1
J_
nd
(20)
The PEF of a district cooling system where all the available heat from plant is used in an absorption chiller to produce the cooling energy is calculated using Eqs. (3), (19) and (20):
PEFn,
where
= PEFF '(I - * ) _ C0Pabs nDN ' (c
1
x=nc
il n
(21)
(22)
is the ratio between the energy efficiency of the electricity generation in a combined cooling and power generation system nCCP,el (heat and power production; heat is extracted from an extraction-condensing turbine) and the energy efficiency of the electricity generation in a thermal power plant when operation is related only to the production of electricity nel (power production; condensing turbine). The total energy efficiency of a cogeneration system where all the extracted heat is used for cold production is nCCP,total. Only natural gas and coal are considered as the fuels in this study. The PEFf are given in Section 1.1.
2 PARAMETERS USED IN THIS ANALYSIS
Based on the equations from Section 1 a mathematical model was developed. The parameters used in this model are divided into two categories:
• Primary energy factors (PEFDC, PEFet PEFf (PEFcoal, PEFng, PEFwh)).
• Technical parameters (ncB , Vdn , nel , nccp.ei ,
nCCP, totah COPabs, COPDFAC, COpCom, COpComjCE).
The PEFs of the fuels considered in this study are shown in Table 1.
Table 1. Primary energy factors of fuels
Fuel PEF Ref.
Coal 1.3 [1]
Natural gas 1.1 [1]
OH 1.1 [1]
Heat from combined gas and steam turbine 0.5 [25]
Heat from steam turbine 0.8 [14] and [25]
Waste heat from industrial process 0.05 [21] and [22]
Electricity (EU average) 2.5 [14]
Diesel, petrol, kerosene 1.19 [24]
Table 2 shows the technical parameters of the chillers considered in this study. The calculations for the PEFdc were made for a single-effect absorption chiller driven by hot water, a double-effect absorption chiller driven by steam, and a direct-fired absorption chiller driven by natural gas or the waste heat from an industrial process. The electrically driven and engine-driven compressor chillers were also included in this study.
In the literature the COPs of different types of chillers were presented at given conditions: for gas absorption and engine-driven compressor chiller the chilled water temperature was 7 °C. The cooling water inlet temperature was 32 °C for DFAC and 30 °C for ICE. For all other types of chillers the chilled water temperature was 6 °C and cooling water inlet temperature 32 °C.
Table 2. Data for the selected chillers
Chiller Cooling power [kW] Chilled water temperature [°C] Cooled water inlet temperature [°C] Driving energy COP Ref.
Steam absorption (DEAC) 2000 6 32 8 bar steam COP = = Qc / Qh = 1.39 [26]
Hot water absorption (SEAC) 2000 6 32 98 °C hot water COP ■- = Qc / Qh = 0.79 [26]
Gas absorption (DFAC) 282 to 2462 7 32 Natural gas (LHV=34 MJ/Sm3) COP = Qc / Qlhv = 1.42 [20]
Compressor chiller 2000 6 32 Electricity COP = Qc / Wel - = 6.6 [26] and [27]
Engine-driven compressor chiller (ICE) 140 7 30 Fossil fuels COP = Qc / =2.7 [24] and [28]
The following energy efficiencies were used in this study [29] to [34]: the total energy efficiency of a combined cooling and power generation system (Vccp.totai), the boiler efficiency (nCs) and the energy efficiency of a district network (nDN). For the purpose of the study they are all assumed to be 0.9. Energy efficiency of the network considers heat gains which affect the primary energy consumption; therefore, more primary energy has to be used to provide the required cooling demand.
3 RESULTS AND DISCUSSIONS
Based on the parametric analysis, introduced in Section 1, this section provides the values for the PEF of a district cooling system (PEFDC) for different types of cold production as a function of:
• the COP of absorption or compressor chillers
(COPabs, COPcom),
• the PEF of the electricity (PEFel),
• the energy efficiency of the electricity generation when operation is related only to the production of electricity (nel).
In this study the ratios x in Eq. (22) are set as 0.95 when single effect absorption chiller is used for cold production and 0.92 or 0.90 when double effect absorption chiller is used. The values of the PEF for the electricity mix considered in this study are 2.5 [17] and 1.3 [9]. The PEFDC in Sections 3.1 and 3.2 were calculated for a constant energy efficiency of the electricity generation when operation is related only to the production of electricity: nel = 0.37 for a steam turbine and nel = 0.55 for a combined cycle of gas and steam turbine.
The COP is highly dependent on operating conditions. The operating conditions are considered in
this study indirectly in Sections 3.1 and 3.2, where the PEFDC as a function of COPs is presented.
3.1 PEFdc as a Function of the COP for the Absorption Chiller
The PEF of a district cooling system (PEFDC) as a function of the COP of an absorption chiller (COPabs) is presented in Fig. 2. In this case, different types of cooling with an absorption chiller were evaluated:
• cooling with an absorption chiller driven by heat from a boiler, where coal (Cscoal) and natural gas (CBng) are used as the fuel,
• cooling with a direct-fired absorption chiller where natural gas (DFACNG) is used as the fuel,
• cooling with an absorption chiller driven by waste heat from an industrial process,
• cooling with combined cooling and power ge^tron, where c°al (CCpcoal,x=o.95nel=0.37) and natural gas (CCPNG x=0 95 0 55 ) are used as the fuel. ' ' '" '
In this Section the ratio x remains the same (x = 0.95) for all the values of the COPci>s. In the following subsections the different values of the ratio x, as a consequence of the different parameters of the heat extracted from the extraction-condensing turbine, are taken into account.
As can be seen from Fig. 2, the PEFDC decreases with an increase of the COPabs (since less heat is required to drive the absorption chiller). The highest value of the PEFDC is achieved when the absorption chiller is driven by the heat from a boiler and when coal is used as the fuel (CBcoal). For this particular case the PEFDC=2.53 for the absorption chiller with the COP = 0.6 and the PEFDC = 1.09 for the absorption chiller with the COP = 1.4.
Fig. 2. Primary energy factor of a district cooling system as a function of the COP for an absorption chiller
The PEFdc can certainly be improved in the case of combined cooling and power generation (CCP - a CHP plant where all amount of heat is used to drive absorption chillers). For instance, when coal (i.e. case CCPcoal) is used as the primary energy source in such a plant, the PEFDC = 0.22 (for the COP=0.6) and the PEFdc = 0.09 (for the COP = 1.4). When natural gas (i.e. case CCPNG) is considered as the primary energy source, the PEFDC=0.27 (for the COP = 0.6) and to PEFdc = 0.12 (for the COP=1.4). The lowest value of PEFdc is achieved when the absorption chiller is driven by waste heat from the industrial process.
3.2 PEFdq as a Function of the COP for the Electrically Driven Compressor Chiller
The PEFdc as a function of the COP for an electrically driven compressor chiller (COPcom) is shown in Fig. 3. In this case different types of cooling with a compressor chiller were evaluated:
• cooling with a compressor chiller driven by electricity mix from the grid (CCmix),
• cooling with a compressor chiller driven by electricity generated in a thermal power plant,
where coal ( CCTPPcan=0.37 and cctppcoaln=0.45 )
and natural gas are used as the fuel
( CCTPP,NG,nel =0.55 and CCTPP,NGn =0.60 ).
From the results in Fig. 3 it is clear that the lowest values of the PEFDC are achieved when a compressor chiller driven by electricity mix from the grid (PEFei = 1.3) is used. For this particular case PEFdc = 0.36 (for the COPcom = 4) and PEFDC = 0.18 (for the COPcom = 8). In the case when the compressor chiller is driven by electricity from a thermal power plant (combined cycle of gas and steam turbine; natural gas is used as the fuel), PEFDC=0.56 (for the COPcom = 4) and PEFDC = 0.28 (for the COPcom = 8)
when nei=0-55. When nei=0.60, the PEFDC for this case decrease to PEFDC=0.51 (for the COPcom=4) and PEFdc = 0.25 (for the COPcom = 8). The highest values of PEFDC are achieved when the cold is produced by a compressor chiller driven by electricity generated in a thermal power plant with nel = 0.37 where coal is used as the primary energy source
( CCTPP,coaltfel =0.37 ).
3.3 PEFdc as a Function of the PEF for Electricity
Fig. 4 shows the PEFDC depending on the PEF of electricity (PEFel) for the typical COPs for the absorption and compressor chillers given in Table 2. Different types of cold production were considered in this analysis:
• the cold production by an electrically driven compressor chiller (CCmix); the PEFelmix ranges from 1.3 to 3.5,
• the cold production by an engine-driven compressor chiller (CCICE),
• the cold production by single- and double-effect absorption chillers; absorption chillers driven by the heat from a steam turbine; (CCPSEAC,coal, CCPDEAC,coal ); the net of a steam turbine ranges from 0.32 to 0.45 and consequently the PEFet of a steam turbine ranges from 2.9 to 4,
• the cold production by single- and double-effect absorption chillers driven by the heat from a combined cycle of gas and steam turbine (CCPseac,ng, CCPdeac,ng); the nei ranges from 0.35 to 0.60 and consequently the PEFet of a gas and steam turbine ranges from 1.8 to 3.1.
From the results in Fig. 4 it is clear that the PEFDC system decreases when increasing the PEFe of the thermal power plant in the case when absorption chillers are used to produce the cold. The PEFe of
Fig. 3. Primary energy factor of a district cooling system as function of the COP for an electrically driven compressor chiller
Fig. 4. Primary energy factor of a district cooling system as a function of the primary energy factors of electricity for different chillers
the thermal power plant is inversely proportional to the energy efficiency of the electricity generation (PEFel = PEFf/ rjel); therefore, by decreasing the nei , the PEFdc decreases. In the case of electrically driven compressor chiller the PEFDC increases when increasing the PEFeh
The PEFd does not affect the PEFDC when an engine-driven compressor chiller (CCICE) is used (because this particular compressor chiller is directly driven by mechanical energy, i.e., no electricity is used). In this particular case the PEFDC can be considered as constant (PEFDC,ICE=0.49). Therefore, when the PEFei > 2.9, the primary energy consumption is lower in the case of an engine-driven compressor chiller, compared to an electrically driven compressor chiller (pEFc^e < PEFCC >m&).
The PEFDC of the electrically driven compressor chiller (CCmix) increases from PEFDC=0.22 (at PEFel = 1.3) to PEFDC=0.59 (at PEFel = 3.5).
The lowest value for the PEFDC is achieved when a double-effect absorption chiller (x=0.92), driven by the heat from a combined cycle of gas and steam turbine, is used (CCPDEAC,NGx=092). For this example PEFDC = 0.24 (at nei = 0.6 and PEFel = 1.8) and PEFDC = 0.13 (at nei=0.35 and PEFel = 3.1). If ratio x for double effect absorption chiller is lower (x = 0.9), the PEFDC increases and is even higher compared to PEFDC when single effect absorption chiller is used.
The energy efficiency of the electricity generation using the steam turbine, when operation is related only to the production of electricity, ranges from 0.32 to 0.45. The lower the nel, the lower is the PEFDC. For this particular case the lowest PEFDC=0.18 is achieved at net=0.35 when double effect absorption chiller (x=0.92) is used.
3.4 PEFdc as a Function of the Energy Efficiency of Electricity Generation in a Thermal Power Plant
The aim of this subsection is to provide information about the PEFDC as a function of the energy efficiency of electricity generation in a thermal power plant, when operation is related only to the production of electricity (nel). In Fig. 5a, the condensing steam turbine (coal is used as the primary energy source) is illustrated when a compressor chiller is driven by the electricity from the thermal power plant, and an extraction-condensing steam turbine is illustrated when absorption chillers are used. In Fig. 5b the combined cycle of gas and steam turbine is taken into account (natural gas is used as the primary energy source).
As is clear from both examples in Fig. 5, the PEFDC is the lowest when a double-effect absorption chiller (CCPDEAC,x=0.92) is used to produce the cold. The PEFDC, with a compressor chiller driven by electricity from the grid (CCmix,), does not depend on the energy efficiency of the thermal power plant (nel). In this particular case, it remains constant (PEFDC = 0.42 at PEFel=2.5 and PEFDC = 0.22 at PEFel = 1.3). By increasing the value of the nel the PEFDC decreases when a compressor chiller driven by electricity from a thermal power plant (CCTPP) is used.
From the results in Fig. 5a, the PEFDC with a single-effect absorption chiller (CCPSEAC,coal,x=0.95) and double-effect absorption chillers (CCPDEAC,coal,x=0.92 and CCPDEAC,coal,x=0.90) are lower, compared to electrically driven compressor chillers (CCmix and CCTPP), for any value of r)ei at given ratios x.
Analysing the results from Fig. 5b, the PEFDC with a compressor chiller driven by electricity from
Fig. 5. Primary energy factors for a district cooling system as a function of the energy efficiency of the electricity generation in a thermal power plant: a) steam turbine (coal), b) combined cycle of gas and steam turbine (natural gas)
a thermal power plant is lower, compared with the one driven by the electricity mix (PEFel=2.5) when nei> 0.44. In the case of cold production with a single-effect absorption chiller (CCPSEAC,NG,x=0.95), the PEFdc is lower compared to electrically driven compressor chiller (PEFel=1.3) until nei> 0.57. Using the double-effect absorption chiller with x=0.92, the PEFdc are lower compared to any electrically driven compressor chiller for any value of neh When double effect absorption chiller with x=0.90 is used, the PEFdc is lower compared to the compressor chiller driven by electricity from the grid (PEFel = 1.3) when nel < 0.55.
4 CONCLUSIONS
In this paper an analysis of the primary energy factors of district cooling systems (PEFDC) for different types of cold production is presented. Several different types of chillers were included in the study: a single-effect absorption chiller, a double-effect absorption chiller, a direct-fired absorption chiller and a compressor chiller driven by electrical energy or by mechanical energy
from an internal combustion engine. The PEFs of the
different fuels were taken from earlier studies. The results of the study reveal that:
• The most primary energy of all the cases discussed in this study is consumed for cold production with an absorption chiller driven by the heat from a boiler (see Fig. 2).
• Cold production with an engine-driven compressor chiller (ICE) compared to an electrically driven compressor chiller (CCmix) is rational only if the PEFel > 2.9 (see Fig. 4).
• Primary energy consumption can be reduced by using heat from CHP plant to drive absorption chillers (heat and electricity are produced in a CHP plant; all of the heat is used for the cold production). The lower the energy efficiency of electricity generation in a thermal power plant (nel), the lower is the PEFDC (see Fig. 5).
• If the energy efficiency of electricity generation in a thermal power plant with a combined cycle of gas and steam turbine is higher than 0.44, the PEFdc is lower for a compressor chiller driven by electricity from a thermal power plant (CCTPP,NG)
, compared to a compressor chiller driven by electricity from the grid with PEFel = 2.5 (see Fig. 4). The PEFdc of a thermal power plant with a combined cycle of gas and steam turbine achieve its lowest value when nel=0.6 (PEFDC = 0.31).
• The heat extracted from the turbine affects the energy efficiency of electricity generation in a CHP plant. When a double-effect absorption chiller is used to produce cold, the energy efficiency of the electricity generation is lower compared to the energy efficiency of electricity generation when the heat to drive a single-effect absorption chiller is extracted. The ratio x, between the energy efficiency of electricity generation in a CHP plant (the heat is extracted from the extraction-condensing turbine) and the energy efficiency of electricity generation in a thermal power plant (condensing turbine) has a significant impact on the PEFDC (see Fig. 5).
• The comparison between the primary energy consumption of different cooling technologies and the cooling with electrically driven compressor chiller when PEFel=2.5 is presented in Appendix A.
The general mathematical model developed in this study can be used as a tool for selecting the most primary energy efficient type of cooling technology in a feasibility study of a new district cooling system implementation. In the future work this mathematical model has to be improved by considering the production of the electricity, heat and cold at the same time (trigeneration system). A very complex general mathematical model for the calculation of the PEFdc will have to take into account the ratio between the heat used for heating and the heat used for cold production. Moreover, it will consider the increase of the consumption of the primary energy as a consequence of heat extraction from the turbine and lastly, the primary energy consumption as a consequence of different parameters of the extracted steam has to be taken into the consideration.
5 NOMENCLATURE
COP coefficient of performance, [-]
PEF primary energy factor, [-]
Q energy, [kWh]
W work, [kWh]
n energy efficiency, [-]
Abbreviations
abs absorption chiller coal coal
com compressor chiller
el electricity
G generator
H heat
F fuel
me mechanical energy
mix electricity mix
P primary energy
trigen trigeneration
6 REFERENCES
[1] Standard EN 15316-4-5:2006. Heating systems in buildings - method for calculation of system energy requirements and system efficiencies - part 4-5 space heating generation systems, the performance and quality of district heating and large volume systems. European Standardization Organizations (CEN), Brussels.
[2] Directive 2012/27/EU of the European Parliament and of the Council of 25 october 2012 on energy efficiency, amending directives 2009/125/EC and 2010/30/EU and repealing directives 2004/8/EC and 2006/32/EC text with EEA relevance, oj 1315, 14.11.2012, p. 1-56.
[3] Directive 2010/31/EU of the European Parliament and of the Council of 19 may 2010 on the energy performance of buildings, oj l 153, 18.6.2010, p. 124-146.
[4] Directive 2009/28/EC of the European Parliament and of the Council on the promotion of the use of energy from renewable sources and amending and subsequently repealing directives 2001/77/EC and 2003/30/EC, oj l/eut l 140, 23 april 2009, p. 16-47.
[5] Standard EN 15203/15315:2006. Energy performance od buildings-overall energy use, CO2 emissions and definition of energy ratings. European Committee for Standardization (CEN), Brussels.
[6] Standard EN 15603:2008. Energy performance of buildings-overall energy use and definition of energy ratings. European Committee for Standardization (CEN), Brussels.
[7] Beretta, G.P., lora, P., Ghoniem, A.F. (2012). Novel approach for fair allocation of primary energy consumption among cogenerated energy-intensive products based on the actual local area production scenario. Energy, vol. 44, no. 1, p. 11071120, D0l:10.1016/j.energy.2012.04.047.
[8] Laverge, J., Janssens, A. (2012). Heat recovery ventilation operation traded off against natural and simple exhaust ventilation in Europe by primary energy factor, carbon dioxide emission, household consumer price and exergy. Energy and Buildings, vol. 50, p. 315-323, D0I:10.1016/j. enbuild.2012.04.005.
[9] Wilby, M.R., Rodríguez González, A.B., Vinagre Díaz, J.J. (2014). Empirical and dynamic primary energy factors. Energy, vol. 73, p. 771-779, D0I:10.1016/j.energy.2014.06.083.
[10] Schicktanz, M.D., Wapler, J., Henning, H.-M. (2011). Primary energy and economic analysis of combined heating, cooling and power systems. Energy, vol. 36, no. 1, p. 575-585, D0I:10.1016/j.energy.2010.10.002.
[11] Fumo, N., Chamra, L.M. (2010). Analysis of combined cooling, heating, and power systems based on source primary energy consumption. Applied Energy, vol. 87, no. 6, p. 2023-2030, DOI:10.1016/j.apenergy.2009.11.014.
[12] Mago, P.J., Chamra, L.M. (2009). Analysis and optimization of CCHP systems based on energy, economical, and environmental considerations. Energy and Buildings, vol. 41, no. 10, p. 1099-1106, DOI:10.1016/j.enbuild.2009.05.014.
[13] Wirgentius, N. (2006). Primary resource factor- for systematic evaluation of heating and cooling options. Euroheat & Power Industry and Utility Forum, Vilnius.
[14] Dalla Rosa, A.; Christensen, J.E. (2011). Low-energy district heating in energy-efficient building areas. Energy, vol. 36, no. 12, p. 6890-6899, DOI:10.1016/j.energy.2011.10.001.
[15] Prandin, M. (2010). Exergy analysis at the community level: Matching supply and demand of heat and electricity in residential buildings. MSc. Project, KTH-Royal Institute of Technology, Stockholm.
[16] Zyczynska, A. (2013). The primary energy factor for the urban heating system with the heat source working in association. Eksploatacja i Niezawodnosc - Maintenance and Reliability, vol. 15, no. 4, p. 458-462.
[17] Jungbauer, J., Serrano Garcia, D., Wallisch, A., Dalin, P., Terouanne, D., Wirgentius, N. (2011). Measurements of individual chiller systems compared to district cooling solutions. ECEEE Conference Proceedings, Toulon.
[18] Riepl, M., Loistl, F., Gurtner, R., Helm, M., Schweigler, C. (2012). Operational performance results of an innovative solar thermal cooling and heating plant. Energy Procedia, vol. 30, p. 974-985, DOI:10.1016/j.egypro.2012.11.110.
[19] Dalin, P., Ivancic, A., Tvarne, A., Penthor, A., Martin, B., Ricci, F. (2006). District cooling-cooling more with less. Euroheat Power, vol. 32, p. 26-27.
[20] Makita, K. (2008). Waste heat energy application for absorption chillers. 3rd International District Cooling Conference & Trade Show, Dubai,
[21] Werner, S. (2006). Guidelines for assessing the efficiency of district heating and district cooling system. EUROHEATCOOL Work package 3, vol. 3, p. 13.
[22] Pout, C. (2011). Proposed carbon emission factors and primary energy factors for SAP 2012. Technical papers supporting SAP 2012, Building Research Establishment, Watford.
[23] Li, Y. (2012). Operation proposal and efficiency analysis of direct-fire absorption chillers biogas produced in the brewer. World Renewable Energy Forum, Denver.
[24] Borjesson, P., Prade, T., Lantz, M., Bjornsson, L. (2015). Energy crop-based biogas as vehicle fuel—the impact of crop selection on energy efficiency and greenhouse gas performance. Energies, vol. 8, no. 6, p. 6033-6058, DOI:10.3390/ en8066033.
[25] Sweden Green Building Council (2014). Treatment of scandinavian district energy systems in LEED. Energy models gor LEED EA credit 1, Sweden Green Building Council, Stockholm.
[26] Poredos, A.; Kitanovski, A. (2011). District heating and cooling for efficient energy supply. International Conference on Electrical and Control Engineering Conference proceeding, p. 5238-5241, Yichang, DOI:10.1109/iceceng.2011.6058201.
[27] Kawasaki Thermal Engineering Co., L. (2014). Absorption Chillers Application and Efficiency, from http://www.khi.co.jp/ corp/kte, accessed on 2016-06-16.
[28] Ponomarev, P., Minav, T., Âman, R., Luostarinen, L. (2015). Integrated electro-hydraulic machine with self-cooling possibilities for non-road mobile machinery. Strojniški vestnik -Journal of Mechanical Engineering, vol. 61, no. 3, p. 207-213, DOI:10.5545/sv-jme.2014.2017.
[29] Lozano, D.E, Martinez-Cazares, G., Mercado-Solis, R. D., Colas, R., Totten, G.E. (2015). Estimation of transient temperature distribution during quenching, via a parabolic model. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 2, p. 107-114, DOI:10.5545/sv-jme.2014.1997.
[30] Rosen, M.A., Le, M.N., Dincer, I. (2005). Efficiency analysis of a cogeneration and district energy system. Applied Thermal Engineering, vol. 25, no. 1, p. 147-159, DOI:10.1016/j. applthermaleng.2004.05.008.
[31] Rezaie, B., Rosen, M.A. (2012). District heating and cooling: Review of technology and potential enhancements. Applied Energy, vol. 93, p. 2-10, DOI:10.1016/j.apenergy.2011.04.020.
[32] Patel, C.T., Patel, B., Patel, V.K. (2013). Efficiency with different gcv of coal and efficiency improvement opportunity in boiler. International Journal of Innovative Research in Science, Engineering and Technology, vol. 2, no. 5, p. 1518-1527.
[33] Durkin, T.H. (2006). Boiler system efficiency. ASHRAE Journal, vol. 48, no. 7, p. 51-57.
[34] Dobersek, D., Goricanec, D. (2009). Optimisation of tree path pipe network with nonlinear optimisation method. Applied Thermal Engineering, vol. 29, no. 8, p. 1584-1591, DOI:10.1016/j.applthermaleng.2008.07.017.
APPENDIX A
The ratio between the primary energy input of different cooling technologies and the primary energy input for the case of cooling with electrically driven compressor chillers (PEFei=2.5) is presented in Fig. A1. The primary energy input is calculated by rearranging the Eq. (1):
QP = PEFdc ' Qdc ,
and the ratio is calculated as:
Qp
PEFn,
Qp
PEFn,
PEFdc at given conditions for different cooling technologies are given in Table A1.
Fig. A1. The primary energy input of different cooling technologies compared to the primary energy input when cooling with electrically
driven compressor chillers (PEFei=2.5)
Table A1. PEFDC for different types of cooling technologies
Type of cooling technology
General Eq. for PEFDCi
CB
PEFn,
1
COPabs -VdN nCB
PEF
DFACng PEFd
DC,NG I rnp „I ^ NG
DFAC IDN y
■ PEF„,
WH
PEF,
DCWH~lcopabs n™
■ PEFW.
PEF,
DC
PEF = I DC ,CB-Pm ] .
The main procedures can be presented in Fig. 1.
(3)
Fig. 1. Flow chart of Isomap space transformation
As described in Section 2.1, in this paper, a high-dimensional «x38 original features array is obtained from the vibration signal. To obtain reasonable dimensionality, a maximum likelihood estimation (MLE) algorithm is introduced to evaluate the intrinsic quantity of the projected dimensionalities before the array is transformed to an n*m one (m< 38) by Isomap.
1.3 DBN-Based Assessment Model
As a probabilistic generative model that can effectively capture the typical information from raw data with various non-linear transformations and approximate complex non-linear functions, DBN is highly appropriate for machine health assessment. Next, the construction and training of DBN model will be introduced.
1.3.1 DBN Architecture
DBN is constructed by stacking a sequence of restricted Boltzmann machines (RBMs) layer by layer [8], as shown in Fig. 2, layer 1 (input layer V) and layer 2 (hidden layer H1) form RBMj, layer 2 (hidden layer Hj) and layer 3 (hidden layer H2) form
RBM2, and so on. Each RBM consists of a hidden and a visible layer, both of which are composed of binary stochastic units that only connect with the units from different layer but do not connect with the ones within the same layer.
The energy of a joint configuration of the visible and hidden units is given:
E (v, h,e) = -Y2wijvihj -jTbv -fatj, (6)
i=i j=1 i=i j=i
where matrix w denotes the weights between visible and hidden layer, vectors a and b are the biases of hidden units hj and visible units v , while 6 = {w, b, a} are the model parameters.
Fig. 2. Schematic architecture of DBN
To further clarify the working process of DBN, take single RBMj as an example. Vi, j,vt e{0,l}, hj e {0,1} . Each component of input array X corresponds to a node of the visible layer, by this way X is input into a RBM1. After a series of calculations, an output array Y of which each component corresponds to a node of the hidden layer is obtained.
1.3.2 DBN Model Training
The distinctive architecture makes it feasible to train DBN by training a series of RBMs with the contrastive divergence (CD) algorithm. The primary training procedure can be summarized as: each RBM layer is trained by using the activation probabilities of the sub-network RBM as the input training data, while the output serves as the input to the next RBM layer. After the unsupervised RBM pre-training, the first layer is fed with the raw input data and is Gaussian-binary restricted Boltzmann machines (GB-RBM) for
real-valued input; the others are binary or Bernoulli-Bernoulli RBM. Finally, the updating rule for the parameter is given:
W ^ W + 8w (vihj0 - vhj ),
-«. « j ' -< y ),
M '-< v. >'
a ^ a + i
b ^ b + £,,
(7)
(8)
where ew, ea and eb denote the learning rates of weights, hidden bias and visible bias, respectively. Details of the training process can be found in literature [24].
For machine health assessment tasks, after the generative pre-training, other typical discriminative, learning procedures which can effectively fine-tune the weights will be combined to improve the performance of DBN. A very effective way that has been successfully confirmed for implementing discriminative fine-tuning is to add an extra layer of variables that denote the desired labels after the last RBM. Then, the back-propagation algorithm similar to that in the standard back propagation neural network will be introduced to adjust all the network weights.
Multiple samples for training and testing are obtained by dividing the dataset, which are further divided into small "mini-batches" for efficient training [22]. To meet the application requirements, the first 2000 min data NL of the transformed dataset ML in projected space, which are all collected under normal conditions are used as training data to train and fine-tune the normal DBN model, that is:
Pi
Training data:
Nl =
P2
Pi
pN
Pi
p N
1.4 Assessment
The specific structure and nonlinear learning procedures of DBN make it very effective for obtaining the intrinsic characteristics from a large number of data. After obtaining the normal DBN model, the entire feature dataset ML that contains normal data as well as fault data are used as testing data to assess the machine whole life health state, that is:
Training data: M =
Pi P2
Pi
pM
Pi
p M
P
N
P
M
where NL e M .
The principal work for DBN-based machine health CAM can be summarized as follows: (1) Pre-process the original vibration signal data, such as eliminating the abnormal data. (2) Extract features of the datasets with Isomap and divide them into training and testing subsets. (3) Build the assessment model based on deep learning and DBN theories. (4) Train DBN with the training feature data NL to obtain a normal DBN assessment model. (5) Apply the normal DBN model and the entire dataset ML for machine health assessment. It should be noted here that steps (3) and (4) will be repeated until the number of epochs is reached. Since being trained with the feature data collected under normal conditions, when the entire dataset is input into the DBN assessment model, its output is the relative probability value of each input belonging to normal condition which is determined by DBN, and its changing trend can reflect the running condition of the tested machine.
Fig. 3 shows the procedures of the proposed CAM.
2 EXPERIMENTS AND ANALYSIS 2.1 Test Rig and Data
The experiment executed bearings run-to-failure tests un-der permanent load was conducted on a special test rig, as shown in Fig. 4, in order to make comparisons with the existing methods and highlight the advantages of the proposed method in this paper; the test data from literature [25] are adopted for the following further research. There are four double row bearings that are force lubricated in the test rig, and the rotation speed of the shaft is applied onto a radial load was kept constant by an alternating current (AC) motor coupled to it via rub belts.
The parameters and operation condition of the bearings are displayed in Table 2. To collect the vibration signal effectively, each bearing was installed with a PCB 353B33 high sensitivity quartz ICP accelerometer on the bearing housing. There are four channels: channels 1 to 4 correspond to bearings 1 to 4, respectively.
Collection of the data was facilitated with NI DAQ Card 6062E. When the test-to-failure
experiment was completed, a dataset which consists of 984 individual vibration signal data files was obtained. The files of ASCII format contain all the information of the four channels. Each file contains 20480 data points with the sampling rate of 20 kHz, and the recording interval of the files is 10 minutes. At the end of the test experiment, failure occurred in the outer race of bearing 1.
Fig. 3. Flowchart of the proposed CAM
Table 2. Bearing parameters and operation condition
Type Number Ball diameter [mm] Contact angle [deg] Rotation speed [RPM] Radial load [kN]
ZA-2115 4 10 0 2000 26.671
Fig. 4. Test rig 2.2 Features Extraction
In view of the outer race failure occurred in bearing 1 at the end of the test experiment, the data of channel 1 were chosen for subsequent experiments. Firstly, pre-processing the data of the bearing vibration signal datasets was carried out, during which three documents that were detected to have abnormal data points were removed. Then the techniques of time domain analysis, frequency domain analysis, and WPT were employed to extract the 38 original features, and a 9810x38 original features array was obtained. Next, the intrinsic quantity of the dimensionalities was calculated with the algorithm of MLE, and the answer was 6. Then the Isomap non-linear dimensionality reduction method was employed to further extract the more representative features; thus, the 9810x38 high-dimensional original features array was converted
into a 9810x6 low-dimensional one, which contains more representative features. Fig. 5 indicates that four features (features 1, 2, 3 and 5) appeared to become irregular at about 5300 and all these 6 features have a significant mutation around the point of 7000. The phenomena described above may suggest that the bearing tested in the experiment could have a fault in the vicinity of these abnormal points.
Although extracted from the same test dataset, the six feature trends in Fig. 5 still exhibit significant inconsistency; therefore, these features are not supposed to be applied individually for machine health assessment as they are not able to provide the consistent and accurate degradation pattern. A more feasible method is to construct an effective model that can fuse all the information of these representative features to obtain reliable and accurate machine health assessment results.
In the next experiments, the 9810x6 features array obtained above will be used to achieve two functions: the former 2000x6 subpart will be applied for training and tuning the assessment model, while the whole dataset will be used to test the model and get the health assessment results of the machine during the entire test period.
2.3 Design of DBN Structure
The numbers of input, hidden nodes, and hidden layers are the most critical parameters in the design of neural networks. As for the DBN-based assessment model, the number of input nodes corresponds to the number of the dimensionality of the array composed of extracted features data, of which the value is 6. Since the one-step-ahead assessment method is chosen, one output node is enough to meet the requirements.
Time [mm] Time [mm]
Fig. 5. The 6 features extracted by Isomap-based dimensionality reduction
The number of hidden nodes determine the ability of neural networks to capture nonlinear patterns from the input data. On the one hand, neural networks with too few hidden nodes may not be competent enough to model the data. On the other hand, too many hidden nodes may bring about over-fitting problems and result in inferior assessment performance. As a generative neural network model with multiple hidden layers, DBN is more powerful in acquiring the complicated relationship from the input data. Therefore, the changeable number of hidden layers and hidden nodes provides the designer of DBN a significant amount of freedom, which also brings many problems. To date, there has been no mature means in theory for determining the number of hidden layers and nodes, which remains an intractable task for the design of DBN structure.
In this paper, the numbers of hidden layers and nodes of the DBN were selected by a series
of experiments. The data for training DBN were divided into two parts: 80% of them were applied for training, while the rest for validation. Given that the numbers of hidden layers of DBNs [13] that had been successfully applied to solve various problems are all more than two, hence, the constructed DBN was initialized with three hidden layers. After several attempts, it was discovered that when the hidden nodes of the three hidden layers are set as 100-50-50, respectively, the evaluation result would be better than other design patterns. Analogously, after comparing with the same analysis method, the architectures of the DBNs with 4, 5 and 6 hidden layers that realized the most satisfactory evaluation results are 100-10050-20, 100-100-50-50-10 and 100-100-100-5020-20, respectively. The assessment outputs of these optimal architectures of the DBNs described above are displayed in Fig. 6, which shows that the DBN with 4 hidden layers could acquire better assessment
= 0.8
0.6
0 -0.2
2.5 -
Mutational point:700
1.5
Abnarmality beginning point:530
Assessment result of 100-50-50 DBN
200
400 600
Time [10min]
800
Assessment result of 100-100-50-50-10 DBN
Mutational point:700
Abnarmality beginning point530
400 600
Time [10min]
1000
1.4
0.9
0.9
0.6
0.4 -
400 600
Time [10min]
Mutational point:700
Abnarmality beginning point:530
Assessment result of 100-100-50-50-20-20 DBN
0
200
400 600
Time [10min]
800
1000
Fig. 6. The assessment results of the 4 different optimal DBNs
1.3
1.2
3 0.4
<
1.1
0.2
200
800
1000
3
0.8
2
0.7
1.5
3 0.5
<
0.5
0.3
0
0
200
800
1000
result than a 3-hidden-layer one, for example, the waveform of the former is more stable than that of the latter before the mutational point, and its mutation at time point 700 is more obvious. While the DBN with the architecture of 100-100-50-50-10 performs best in the experiment by virtue of its smoothest, clear and reasonably trended curve, with the number of hidden layers increasing to 6, the performance of DBN had become worse.
The phenomenon above can be explained as that DBN with less than five hidden layers may not have enough power to model the data; however, when there are more than 5 hidden layers the over-fitting problems may result in inferior assessment performance. It is, therefore, the conclusion that the 5-hidden-layer DBN with the architecture of 100-100-50-50-10 is the best model for the work of machine health assessment. (It is worth noting that since the value of each point of the waveforms is a relative one, the direction the
waveforms trended (up or down) does not affect the quality of the assessment results.)
3 COMPARISON EXPERIMENTS
3.1 Comparison Experiments of Dimensionality Reduction
Isomap, adopted in this paper, is a nonlinear global dimensionality reduction method; next, the comparison experiments with linear dimensionality reduction method PCA and nonlinear local dimensionality reduction technique Laplacian Eigenmaps that have been used for feature extraction were conducted.
In order to make fair and reasonable comparisons, the most suitable DBN architectures for PCA and Laplacian Eigenmaps are constructed, which are 100100-50-50-5 and 100-100-50-50-10, respectively. The assessment results of DBNs with these two different dimensionality reduction methods are
1.15
1.05
1
0.9
0.85
2.4 2.2
1.8
1.4 1.2 1
0.8
200
400 600
Time [10min]
800
1000
Assessment result of Laplacian Eigenmaps-based dimensionality reduction
Mutational point:647
Abnarmality beginning point:530
0
200
800
1000
400 600
Time [10min]
Fig. 7. The assessment results of PCA and Laplacian Eigenmaps-based dimensionality reduction method
1.6 1.4 1.2 1
0.8 0.6 0.4 0.2
Mutational point:700
f "Hi**
Abnarmality beginning point:530
Assessment result of BPNN
400 600
Time [10min]
Mutational point:700
Abnarmality beginning point:530
Assessment result of HMM
0 200 400 600 800 100C
Time [10min]
Fig. 8. The assessment results of BPNN and HMM assessment models
1.1
0.95
0.8
0
0
0
200
800
1000
2.6
0
-1000 -
2
0) -2000 -
g -3000
1.6
-4000
-5000
displayed in Fig. 7, from which it can be determined that the resulting waveform of PCA-based method has the same abnormality beginning point (at time point 530) and mutational point (at time point 700) with that of the Isomap-based method, but the former fluctuates much more severely before abnormality beginning point and becomes disorganized at the end. While, the resulting waveform of Laplacian Eigenmaps -based method performs better for the two problems mentioned above, but the overly bizarre mutational point (at time point 647) also demonstrated the defect of this method. The analysis and comparisons above surely proved the correctness of the Isomap-based dimensionality reduction method.
3.2 Comparison Experiments of Assessment Models
Next, two other kinds of intelligent algorithms HMM and BPNN that have already been successfully and widely applied in fault diagnosis, data classification, are also introduced to conduct the similar experiments under the corresponding same conditions. As an extension of Markov chains, HMM that contains a finite number of states, where each state generates an observation at a certain time point, is a state-of-the-art technique for model recognition due to its elegant mathematical structure and the availability of computer implementation [1]. While a BPNN model is composed of many idealized layers of nodes and specified by the node characteristics (weights), the learning rules (transfer functions, always the sigmoid function), network interconnection geometry (different layers), and dimensionality (number of layers and nodes), of which the learning feeds back into the model to change the weights of nodes between layers in order to decrease errors between predicted and measured values [26].
In order to make the results of the experiments more persuasive, the dimensionality reduction method is kept the same as Isomap only the assessment models were replaced by BPNN and HMM, respectively. The assessment results of these two models are shown in Fig. 8, which indicates that BPNN can accurately identify the anomaly and mutation at time points 530 and 700, respectively, but it performs disorderedly after the mutational point when the waveform trends to the opposite direction. While the waveform of HMM is very smooth at the beginning, it drops too fast after the abnormality beginning point 530, which does not accord with the actual situation, and the transformation at mutational point 700 is too vigorous so that the waveform graph disconnected there; additionally, its tail waveform is also very confusing.
4 ASSESSMENT RESULTS
After a series of experiments, in which the comparisons with two popular dimensionality reduction methods (PCA and Laplacian Eigenmaps-based) and two intelligent assessment algorithms (HMM and BPNN-based), the proposed CAM in this paper, which employs Isomap and the DBN with architecture of 100-100-50-50-10, was verified to be very effective and accurate for machine health assessment work.
Bearing health assessment result of the proposed CAM
1.6 -
1.4 -
Sharp deterioration beginning point: 940
Slight fault beginning point: 530
200
800
400 600
Time [10 min]
Fig. 9. Assessment result of the proposed CAM
As shown in Fig. 9, after a long time of normal running, the proposed CAM detected the early slight degradation of the bearing at time point 530, indicating where the slight faults of wear, pitting or overheat may begin occurring. Then the CAM monitored a mutation signal (at point 700) which presents that crackle, fatigue spalling or some other serious faults might occur there. Finally, it can be observed from the end of the waveform that after the time point 940, the condition of the bearing began to deteriorate so sharply that it could not continue working.
Compared with the real test situations, the assessment results of the proposed method in this paper are very much in line with the actual running conditions of the bearings.
5 CONCLUSIONS
Machine health assessment is playing an increasingly important role, which can provide diverse benefits including improved safety, improved reliability and reduced costs for operation and maintenance of complex manufacturing systems. This paper presents a novel combined assessment method (CAM) for
2
1.8 -
1.2
0.8 -
assessing the health condition of the target machine (rolling-element bearings). To deal with the non-stationary property of the vibration signals collected from machines, three widely used techniques involving time domain analysis, frequency domain analysis, and WPT are adopted to extract 38 original features of the datasets. Then, the Isomap nonlinear global algorithm is adopted for dimensionality reduction and extracting more representative features. Next, the acquired low-dimensional features array is input into the well-trained DBN model to evaluate the performance status of the bearing. Finally, the bearing-accelerated degradation data from Cincinnati University were investigated for further research, through the comparison experiments with two popular dimensionality reduction methods (PCA and Laplacian Eigenmaps-based), and two intelligent assessment algorithms (HMM and BPNN-based), the proposed CAM is proved more sensitive to the incipient fault of rolling bearings and more effective for bearing performance degradation evaluation. In future work, CAM can be applied in other fields for evaluation or classification.
6 ACKNOWLEDGMENT
This work is supported by the National Natural Science Foundation of China (No. 51374264) and State Key Laboratory of Coal Mine Disaster Dynamics and Control Open Project.
7 REFERENCES
[1] Yu, J. (2012). Health condition monitoring of machines based on hidden Markov Model and con-tribution analysis. IEEE Transactions on Instrumentation and Measurement, vol. 61, no. 8, p. 2200-2211, D0l:10.1109/TIM.2012.2184015.
[2] Loutridis, S. (2006). Instantaneous energy density as a feature for gear fault detection. Mechanical Systems and Signal Processing, vol. 20, no. 5, p. 1239-1253, D0I:10.1016/j. ymssp.2004.12.001.
[3] Ozturk, H., Sabuncu, M., Yesilyurt, I. (2008). Early detection of pitting damage in gears using mean frequency of scalogram. Journal of Vibration and Control, vol. 14, no. 4, p. 469-484, D0I:10.1177/1077546307080026.
[4] Loutridis, S.J. (2008). Self-similarity in vibration time series: application to gear fault diagnostics. Journal of Vibration and Acoustics, vol. 130, no. 3, p. 569-583, D0I:10.1115/1.2827449.
[5] Yu, D., Yang, Y., Cheng, J. (2007). Application of time-frequency entropy method based on Hilbert-Huang transform to gear fault diagnosis. Measurement, vol. 40, no. 9-10, p. 823-830, D0I:10.1016/j.measurement.2007.03.004.
[6] Feng, Z., Zuo, M.J., Chu, F. (2010). Application of regularization dimension to gear damage assessment. Mechanical
Systems and Signal Processing, vol. 24, no. 4, p. 1081-1098, D0l:10.1016/j.ymssp.2009.08.006.
[7] Do, V.T., Nguyen, L.C. (2016). Adaptive empirical mode decomposition for bearing fault detection. Strojniški vestnik -Journal of Mechanical Engineering, vol. 62, no. 5, p. 281-290, D0I:10.5545/sv-jme.2015.3079.
[8] Yan, J., Guo, C. (2011). A dynamic multi-scale Markov model based methodology for remaining life prediction. Mechanical Systems and Signal Processing, vol. 25, no. 4, p. 1364-1376, D0I:10.1016/j.ymssp.2010.10.018.
[9] Tran, V.T., Yang, B.-S., Oh, M.-S., Tan, A.C.C. (2009). Fault diagnosis of induction motor based on decision trees and adaptive neuro-fuzzy inference. Expert Systems with Applications, vol. 36, no. 2, p. 1840-1849, D0I:10.1016/j. eswa.2007.12.010.
[10] Hinton, G.E., Salakhutdinov, R.R. (2006). Reducing the dimensionality of data with neural networks. Science, vol. 313, no. 5786, p. 504-507, D0I:10.1126/science.1127647.
[11] Lecun, Y., Bottou, L., Bengio, Y., Haffner P. (1998). Gradient-based learning applied to document recognition. Proceedings of the IEEE, vol. 86, no. 11, p. 2278-2324, D0I:10.1109/5.726791.
[12] Li, Z., Ma, Z., Liu, Y., Teng, W., Jiang, R. (2015). Crack fault detection for a gearbox using discrete wavelet transform and an adaptive resonance theory neural network. Strojniški vestnik - Journal of Mechanical Engineering, vol. 61, no. 1, p. 63-73, D0I:10.5545/sv-jme.2014.1769.
[13] Gan, M., Wang, C., Zhu, C. (2016). Construction of hierarchical diagnosis network based on deep learning and its application in the fault pattern recognition of rolling element bearings. Mechanical Systems and Signal Processing, vol. 72-73, p. 92104, D0I:10.1016/j.ymssp.2015.11.014.
[14] He, Q.B. (2013). Vibration signal classification by wavelet packet energy flow manifold learning. Journal of Sound and Vibration, vol. 332, no. 7, p. 1881-1894, D0I:10.1016/j. jsv.2012.11.006.
[15] Klein, R., Ingman, D., Braun, S. (2001). Non-stationary signals: phase-energy approach theory and simulations. Mechanical Systems and Signal Processing, vol. 15, no. 6, p. 1061-1089, D0I:10.1006/mssp.2001.1398.
[16] Baydar, N., Ball, A. (2001). A comparative study of acoustic and vibration signals in detection of gear fail-ures using Wigner-Ville distribution. Mechanical Systems and Signal Processing, vol. 15, no. 6, p. 1091-1107, D0I:10.1006/mssp.2000.1338.
[17] Wan, X., Wang, D., Tse, P.W. Xu, G., Zhang, Q. (2016). A critical study of different dimensionality reduction methods for gear crack degradation assessment under different operating conditions. Measurement, vol. 78, p. 138-150, D0I:10.1016/j. measurement.2015.09.032.
[18] Wang, H. (2012). Block principal component analysis with L1-norm for image analysis. Pattern Recognition Letters, vol. 33, no. 5, p. 537-542, D0I:10.1016/j.patrec.2011.11.029.
[19] Borchers, B., Young, J.G. (2007). Implementation of a primal-dual method for SDP on a shared memory parallel architecture. Computational Optimization and Applications, vol. 37, no. 3, p. 355-369, D0I:10.1007/s10589-007-9030-3.
[20] Hemmati, F., Orfali, W., Gadala, M.S. (2016). Roller bearing acoustic signature extraction by wavelet packet transform,
applications in fault detection and size estimation. Applied Acoustics, vol. 104, p. 101-118, D0l:10.1016/j. apacoust.2015.11.003.
[21] Hauberg, S: (2015). Principal curves on Riemannian manifolds. IEEE Transactions on Pattern Analysis and Machine Intelligence, p. 1915-1921.
[22] Benkedjouh, T., Medjaher, K., Zerhouni, N., Rechak, S. (2013). Remaining useful life estimation based on nonlinear feature reduction and support vector regression. Engineering Applications of Artificial Intelligence, vol. 26, no. 7, p. 17511760, D0I:10.1016/j.engappai.2013.02.006.
[23] Dijkstra, E.W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, vol. 1, no. 1, p. 269271, D0I:10.1007/BF01386390.
[24] Tran, V.T., AlThobiani, F., Ball, A. (2014). An approach to fault diagnosis of reciprocating compressor valves using Teager-Kaiser energy operator and deep belief networks. Expert Systems with Applications, vol. 41, no. 9, p. 4113-4122, D0l:10.1016/j.eswa.2013.12.026.
[25] Lee, J., Qiu, H., Yu, G., Lin, J. (2007). Rexnord Technical Services: Bearing Data Set. Moffett Field, University of Cincinnati. NASA Ames Prognostics Data Repository, NASA Ames from: http://data-acoustics.com/measurements/ bearing-faults/bearing-4/, accessed on 2016-10-20.
[26] Meng, Z., Xu, Y., Zheng, Y., Zhu, Y., Jia, Y., Chen, S. (2014). Inversion of lunar regolith layer thickness with CELMS data using BPNN method. Planetary and Space Science, vol. 101, p. 1-11, D0I:10.1016/j.pss.2014.05.020.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)12, 751-756 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.3060 Original Scientific Paper
Received for review: 2015-09-25 Received revised form: 2016-05-03 Accepted for publication: 2016-05-25
On the Fatigue of Steel Catenary Risers
Nnamdi Onochie Chibueze - Chinwuba Victor Ossia* - John Umunna Okoli
University of Port Harcourt, Offshore Technology Institute, Nigeria
Steel catenary risers (SCR), though an appropriate and cost effective tool for deepwater development, are fraught with difficult fatigue challenges derived from loads due to the actions of wind, wave and current. In this study, a model SCR-AB with end-A attached to floating production, storage and offloading vessels (FPSO) and end-B anchored to the seabed was created using OrcaFlex© finite element software. Also attached to this model is a 104 kg 6D buoy of 6 m height and volume 20 m3, which was fixed to the SCR at 900 m of the riser length. The various sections of the homogenous 2200 m long SCR pipe include 50 m flex joint section, 350 m strake length, 1300 m riser pipe length, and 500 m flowline section. Metrological ocean data, geotechnical data and SCR data were input into the model which was subjected to the Ochi-Hubble and Joint North Sea Wave Project (JONSWAP) wave spectra. Analyses to determine the influence of fatigue damage on the model was made using system simulation under both static and dynamic modes. The results for JONSWAP wave spectra showed the fatigue life of the SCR to be 1.8 years at the touch-down-point (TDP) while that of Ochi-Hubble wave spectra showed a fatigue life of 13.6 years at SCR TDP, differing by a factor of 7.5. Furthermore, the S-N curves obtained from both wave spectra corroborated an inverse relationship between Stress values S and the number of stress cycles to failure N on a log-log scale. Keywords: steel catenary risers, touchdown point, Ochi-Hubble, JONSWAP
Highlights
• Fatigue life at SCR TDP and SCR hang-off obtained are very low, indicative of higher stresses derived from FPSO dynamics, ocean current, wave spectra, and other environmental loadings.
• There is a sharp rise in the SCR tension (mooring line effective tension) from zero at FPSO hang-off point to a maximum of 30*106 kN at 5 m arc length due to the reasons mentioned above.
• JONSWAP and Ochi-Hubble wave spectra used for SCR design are used to determine the SCR fatigue life, tension, cost, etc. by the oil and gas industry operators.
• JONSWAP wave spectrum is often used by engineers for SCR designs in the Gulf of Guinea leads to overdesign with the attendant cost (time, material, logistics, etc.) and is compared to Ochi-Hubble wave spectrum, which is closer to the real-time wave spectra.
0 INTRODUCTION
0.1 Background
Oil & gas exploration and production activities in deep waters of the Gulf of Mexico, Gulf of Guinea, Brazil, and the North Sea have increased dramatically in recent years, nearly doubling, compared to the activities of a decade ago [1]. Regardless of the floating platform concept adopted for many offshore field development activities, there is always need for a riser system, which has played a massive role as offshore infrastructure [2]. There is, therefore, a need to ensure that the integrity of riser systems is not compromised; that is, ensuring high fatigue life in performance. In West African deepwater, for example, steel catenary risers (SCR) are commonly used because of the prevalent conditions, though fatigue remains a challenge [3]. SCR is a cost-effective alternative for oil and gas export and for water injection lines on deepwater fields, where large-diameter flexible risers present technical and economic limitations [4]. Among the riser concepts (e.g. flexible risers, hybrid risers), SCRs have enjoyed
widespread acceptability for deep and ultra-deepwater applications in recent years. This is due to their relatively low cost, conceptual simplicity, significant structural capacity, and ease of fabrication and installation. However, their fatigue response has been a concern to the industry, especially in West Africa [3]. SCR fatigue damage can be caused by vortex-induced vibration (VIV) and vortex-induced motion (VIM) at the touch-down-point (TDP) [3]. Fatigue of SCR at TDP has been one of the most challenging problems in deep water due to interactions with the seabed [5]. The fatigue damage induced by VIV and VIM can be fatal to the riser [6].
0.2 Wave Environment
Harsh environments influence motions of floating platforms, which in turn affect the riser system dynamics. For instance, the application of SCRs with semi-submersibles or floating production, storage and offloading vessels (FPSOs) in harsh deepwater environments presents design challenges due to large wave-induced motions on the platform, and large vessel offsets caused by wind, currents, and slow-drift
*Corr. Author's Address: Offshore Technology Institute, University of Port Harcourt, Port Harcourt, Nigeria; ossiacv@otiuniport.org
751
wave motions. The resulting large heave motions of the vessel cause buckling and fatigue-related issues of the riser at TDP [7]. The wave spectra adopted in SCR design becomes crucial. The two-parameter Pierson-Moskowitz wave spectrum (Eqs. (1) and (2)) is recommended by the American Bureau of Shipping [8] for an open ocean, such as the Atlantic.
, , 5H>p
r \4 m
or,
HC ( 2n
4 nco
exp
^4
c
(1)
. (2)
The Joint North Sea Wave Project (JONSWAP) spectrum (Eqs. (3) and (4)) modified the Pierson-Moskowitz spectrum [9] above to account for regions with geographical boundaries the limit the fetch in wave generation.
, , 5H>4 Sn(œ) = ^ 5P exP
16œ3
where: a = exp
f v œ
œ
y" (1 - 0.287ln/), (3)
)2
2H 1-1
(12)
where w is weight per unit length of riser and H is a constant representing horizontal component of tension.
1 MATERIALS AND METHODS
To study the influence of fatigue in SCR adequately, two wave application modes (Ochi-Hubble wave spectra and JONSWAP wave spectra) were adopted using Orcina OrcaFlex software. The software quality assurance is based on the requirements of BS 5750, 5887, and 5515, and on Defence Standard 00-16. Its validation checks include checking the program against mathematical theory, the real world, and against other institutional computer programs.
Larsen [14] reported results from static and dynamic analyses of standard flexible riser systems performed by eleven different institutions using their in-house software. OrcaFlex developers took part in the major comparison of flexible riser programs initiated by the International Ship and offshore Structures Congress (ISSC) "Slender Marine Structures" committee. The study contains a complete description of the test cases and results that provide a benchmark test for users and programmers of computer code for flexible riser analysis. The results showed good agreement between OrcaFlex and other established programs with a 10 % to 15 % coefficient of variation.
A large-scale test has been carried out in Lake Pend Oreille under a joint industry project by scientific marine services (SMS) [15] for PMB Engineering; Australia. OrcaFlex results showed very close agreement with the test results, including details of SCR behaviour at touchdown. A major validation exercise was carried out to compare OrcaFlex with model test results obtained from a joint industry project [16]. Good agreement was obtained.
In this study, an SCR [17] and [18] comprising a line with the end-A of the riser attached to the FPSO and end-B anchored to the seabed as in Fig. 1 is modelled by carrying out the following tasks using the OrcaFlex post-processor software: a) Set the SCR line length to be 2,200 m and as a
homogenous pipe with flex joint, strake, riser
Fig. 1. Typical steel catenary riser between an FPSO and seabed
b) A 6D buoy is modelled and anchored onto the SCR-line.
c) The design parameters and environmental properties, which included the wave train of both Ochi-Hubble, a two parameter wave spectra that fully describes the typical sea state (valid for the swell and wind-sea conditions) of the Gulf of Guinea and that of JONSWAP were inputs in the different cases.
d) The line modelled as a riser is given real life data-values. Typical SCR data imputed to this study include: pipe of outer diameter 273.1 mm, inner diameter 205.7 mm, and thickness 67.4 mm, density 7850 kg/m3, elastic modulus 200 GPa, second moment of inertia 8.67 x10-5, mass per unit length 128.2 kg/m, bending stiffness 3.27x106 kN, axial stiffness 1.73x104 kN, torsional stiffness 10 kNm2/rad, design life 20 years, bulk modulus 160 GPa, yield strength 483 MPa, tensile strength 562 MPa, kinematic viscosity 1.35 X10-6 m2/s; and the SCR-contents, which is a crude oil of density 800 kg/m3 and flow rate of 0.1 kg/s.
e) Global analysis: This describes the overall static and dynamic structural behaviour of the riser by exposing the system to a stationary environmental loading condition [19]:
• Static analysis investigates the equilibrium position of the model without considering
w
environmental forces, weight, buoyancy and drag forces. After running the static analysis of this model, there was a static convergence for the model after a maximum of 100 iterations.
• The dynamic analysis is a time simulation of the SCR-model motion for a specified period, starting from the static analysis position. The real life environment data for current, wind and seabed conditions were used as inputs. Typical environment data used as input parameters in this study include: (i) metrological Ocean data for 100 years return period: significant wave height 4 m, most probable maximum wave height 7.8 m, crest of maximum individual wave 4.5 m, peak spectral period 15.5 s, mean zero crossing 7.17 s, mean wave direction 189°, (ii) wind data 10 min wind speed 20.4 m/s, wind direction 180 Set the SCR line length to be 2,200 m and as a homogenous pipe with flex joint, strake, riser pipe and flowline lengths modelled on the line (in that order) with the height of 1000 m set above the seabed with an azimuth angle of 180°, (iii) current data: current speed (below 1000 m depth) 0.15, surface directions 180°, (iv) sea water data: sea surface temperatures 25 °C, near-bed sea temperatures 4 °C, sea water density 1025 kg/m3, depth of water 1000 m, (v) seabed data: soil friction coefficient 0.45 (nonlinear soil model was adopted). In running this dynamic analysis, the simulation time was first set at 16 s and observed before being run for 3 h, the severest sea state in the Gulf of Guinea.
f) Fatigue analysis was carried out on the SCR model for both Ochi-Hubble and JONSWAP wave spectra for comparison. Before the fatigue analysis was performed, a set of Orcaflex simulation files that model the system was prepared under possible real-time load conditions. This load case was assigned an exposure level of 20 years (175,200 h) riser design life.
2 RESULTS AND DISCUSSION
2.1 Fatigue Life - Arc Length Analysis
From Fig. 2 based on the JONSWAP wave spectrum,
it can be observed that:
• hang-off fatigue life, HOFL is 10 years;
• fatigue life range, FL; 500 years < FL < 5,000 years;
• TDP fatigue life, TDPFL = 1.8 years at 1828.8 m;
• maximum fatigue life beyond TDP FLmax = 10,000 years.
From Fig. 3, based on Ochi-Hubble wave spectrum, it can be observed that:
• hang-off fatigue life; HOFL ~ 7 years;
• fatigue life range, FL; 100 years < FL < 5,000 years;
• TDP fatigue life; TDPFL = 13.6 years at 1828.8 m;
• maximum fatigue life beyond TDP FLmax = 10,000 years.
^MkISj
Hang-off jDP 1.8
0 304.8 609.6 914.4 1219.2 1524.0 1828.8 2133.6
Length from hang-off point [m]
Fig. 2. Fatigue life versus SCR length from hang-off point in JONSWAP wave spectra
From Fig. 2, the fatigue life of the SCR when subjected to stress due to environmental loadings is 1.8 years, or approximately 2 years, which is very low. Low fatigue life at the SCR TDP zone is indicative of severe fatigue damage therein. From Fig. 3, the SCR fatigue life is showing a high value of 13.6 years, or approximately 14 years, though less than the proposed 20-year design life. However, Figs. 2 and 3 show that the least fatigue life is observed in the SCR hangoff and TDP section. Hence, from Figs. 2 and 3, the fatigue life of SCR in Ochi-Hubble wave spectra is 7.5 times that of SCR fatigue life in JONSWAP wave spectra. This implies that the common practice of SCR designs based on JONSWAP wave spectrum applied in the Gulf of Guinea, where the real-time wave pattern follows the Ochi-Hubble spectra is an overdesign and the cost implication is of order 7.5 high. Figs. 2 and 3 showed the same order of fatigue life beyond TDP arc length. The low hang-off fatigue life is corroborated by the sharp rise in SCR tension from (0 m, -30E6 kN) through (2.5 m, 0 kN) to (4.5 m, 29E6 kN) along the SCR arc length in Fig. 4.
ft /^Naí r
13.6
* Hang-off TDP
0 304.8 609.6 914.4 1219.2 1524.0 1828.8 2133.6
Length from hang-off point [m] Fig. 3. Fatigue life versus SCR length from hang-off point in Ochi-Hubble wave spectra
to
o
f ■>
o
CD 20E6
to
O 10E6
o
00
II
0
c o -10E6
(/)
c
(1)
ve -20E6
o
;
r
.......I.......— Minimum — Maximum — Mean — Allowable Tension Compression Limit |
| (24.00, 11.243E6) |
20 40 60 80
Arc length [m]
Fig. 4. SCR tension versus arc length curve
2.2 Stress (S) - Cycles (N) to Failures
The S-N curves in Figs. 5 and 6 showed that higher stress values (S) caused a lower number of cycles to failure (N), that is, lower endurance limits.
472.0, 92.9E3)
1000 10E3
Cycles to failure (N) [-]
The fall rate of stress with respect to cycles is higher with the JONSWAP spectra than the Ochi-Hubble spectra. Comparing the slopes of the S-N curves, it can be observed that the fall-rate for JONSWAP spectra is almost of order > 100 on a loglog scale relative to the Ochi-Hubble spectra. This implies that when the stresses due to environmental loading vary, it will have an effect on the riser life before failure.
1E6 1E7 1E8
Cycles to failure (N) [-]
Fig. 6. Stress S versus cycles to failure N curve for SCR in Ochi-Hubble wave spectra
3 CONCLUSION
Using an SCR model in Ochi-Hubble wave spectra gives a higher fatigue life than JONSWAP spectra, which are often used by design companies. With a higher SCR fatigue life in Ochi-Hubble wave spectra, there is no overdesign since it thoroughly describes the actual sea state in the Gulf of Guinea.
SCR modelled in JONSWAP wave spectra gave a misrepresentation of the SCR fatigue life leading to overdesign, which implies high cost both on the side of multinational companies and contractors.
Stress variations due to environmental loading contributed immensely to SCR fatigue damage.
The high tension observed at the SCR hang-off, close to the FPSO, with both wave-spectra can be attributed to the combined action of waves, wind and current, as well as, FPSO heavy motion.
4 NOMENCLATURE
Fig. 5. Stress S versus cycles to failure N curve for SCR in JONSWAP wave spectra
A C Ca Cd
Fx F
drag area, [m2]
viscosity damping factor, [Ns/m] added mass coefficient of body drag coefficient for the body body force in x,y,z directions, [N/m3] fluid force, [N]
Hs significant wave height, [m]
M mass, [kg]
P pressure, [Pa]
Sn wave spectral density, [m2/Hz]
T average zero up-crossing period, [s]
V volume, [m3]
^r fluid velocity relative to body, [m/s]
a acceleration, [m/s2]
ar fluid -body relative acceleration, [m/s2]
aw fluid acceleration rel. to earth, [m/s2]
g acceleration due to gravity, [m/s2]
j integer counter
k stiffness, [N/m]
u, v, w x, y, z axis velocity components, [m/s]
t time, [s]
a circular frequency of the wave, [rad/s]
ap peak frequency for peak spectrum; [rad/s]
Y peaked-ness parameter (1 to 7).
a parameter
P density, [kg/m3]
M dynamic viscosity, [Pas]
A Mass of fluid displaced by body, [kg]
1 Ochi-Hubble's parameters 5 & 6
5 REFERENCES
[1] Song, R., Stanton, P. (2007). Deepwater tieback SCR, unique design challenges and solutions. Offshore Technology Conference, paper OTC-18524-MS Houston DOI:10.4043/18524-MS.
[2] Lim, F. (2006). Installation of risers in deepwaters. 4th Petromln Deepwater and Subsea Technology Conference and Exhibition, Kuala Lumpur.
[3] Song, R., Stanton, P. (2009). Advances in deepwater steel catenary riser technology state of art: Part II-Analysis. ASME 28th International Conference on Ocean, Offshore and Arctic Engineering, Honolulu, p. 285-296, D0I:10.1115/0MAE2009-79405, DOI:10.1115/OMAE2009-79405.
[4] Bai, Y., Bai, Q. (2010). Handbook of Subsea Structural Engineering, 5th Edition, Gulf Professional Publishing, Burlington, p. 847-849.
[5] Li, F.Z., Low, Y.M. (2012). Fatigue reliability analysis of a steel catenary riser at the touchdown point incorporating soil model
Uncertainties. Applied Ocean Research, vol. 38, p. 100-110, DOI:10.1016/j.apor.2012.07.005.
[6] Bai, Y. (2001). Pipelines and Risers, vol. 3. Elsevier Ocean Engineering Books Series, Bhattacharya, R, McCormick, M.E. (eds.), Elsevier Sciences Limited, Kidlington, p. 84-110.
[7] Xia, J. (2009). Weight - Optimised Steel Catenary Risers and their Application in Harsh Deepwater Environment, PhD Thesis, University of Strathclyde, Strathclyde.
[8] American Bureau of Shipping (2010). Guide for Spectral-Based Fatigue Analysis for Floating, Production, Storage and Offloading (FPSO) Installations, update 2014, p. 8-12, Houston.
[9] Vyzikas, T. (2014). Application of Numerical Models and Codes. A Best Practice Report prepared as part of the MERIFIC Project - Marine Energy in Far Peripheral and Island Communities (MERIFIC), University of Plymouth, Plymouth, p. 56-60.
[10] Ochi, M.K., Hubble, E.N. (1976). Six-Parameters Wave Spectra. Coastal Engineering Proceedings, p. 301-328.
[11] Orcina Ltd (2013). OrcaFlex version 9.7. Software Technical Specification. Ulverston, Cumbria.
[12] Orcina Ltd (2013). Orcina Orcaflex, retrieved from: http://www. orcina.com/SoftwareProducts/OrcaFlex/index.php, accessed on 2013-12-22.
[13] Orcina Ltd (2013). OrcaFlex Manual: Version 9.7a. Orcina, Ulverston, Cumbria.
[14] Larsen, C.M. (1992). Flexible riser analysis - comparison of results from computer programs. Marine Structures, vol. 5, no. 2-3, p. 103-119, DOI:10.1016/0951-8339(92)90024-J.
[15] Scientific Marine Services, Inc (1998). Highly Compliant Riser Large Scale Model Test Joint Industry Project Data Reduction and Report SMS Project 97-504 PMB Engineering, retrieved from: http://www.bsee.gov/Technology-and-Research/Technology-Assessment-Programs/Reports/300-399/305AB/, accessed on 2013-12-21.
[16] Hartnup, G.C., Airey, R.G., Fraser, J.M., (1987). Model basin testing of flexible marine risers. Proceedings of the 6th International Offshore Mechanics and Arctic Engineering Symposium & Exhibition, Houston.
[17] Bai, Y., Bai, Q. (2005). Subsea Pipelines and Risers, (1st ed.). Elsevier, Amsterdam.
[18] Bai, Y., Bai, Q. (2012). Subsea Engineering Handbook. Gulf Professional Publishers (Elsevier), Waltham.
[19] DNV-OS-F201 (2010). Dynamic Risers, offshore standard. Det Norske Veritas, Oslo.
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)12, 757-763 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2015.3211 Original Scientific Paper
A New Multi-Body Dynamic Model of a Deep Ocean Mining Vehicle-Pipeline-Ship System and Simulation of Its Integrated Motion
Yu Dai123* - Liping Pang1 - Lisong Chen1 - Xiang Zhu1 - Tao Zhang1
1 Central South University, College of Mechanical and Electrical Engineering, China
2 Shanghai Jiao Tong University, State Key Laboratory of Ocean Engineering, China
3 Central South University, Shenzhen Research Institute, China
In this paper, a three-dimensional multi-body dynamic model of a seafloor mining vehicle is developed in the RecurDyn simulation program by integrating the developed terramechanics model of the seafloor sediment. The discrete element method (DEM) is employed to model a mining pipeline as discrete rigid elements linked by flexible connectors. An innovative user-defined subroutine for parametrically and rapidly automatically building a DEM model for such long pipeline is developed in the C# language in RecurDyn. A new 3D multi-body dynamic model of a total deep ocean mining vehicle-pipeline-ship system is established with the integration of developed models for subsystems. The integrated motion simulations are achieved and demonstrate that the synchronized stable motion of the total system can be well maintained during operation processes.
Keywords: deep ocean mining system, terramechanics model, multi-body dynamic model, discrete element model, user-defined subroutine, integrated motion simulation
Highlights
• A multi-body dynamic model of a seafloor mining vehicle has been developed by integrating an innovative user-defined subroutine for characterizing the terramechanics model of the seafloor sediment.
• An innovative user-defined subroutine for parametrically and automatically building a multi-body discrete element model of an extremely long mining pipeline has been developed.
• A new multi-body dynamic model of a deep ocean mining vehicle-pipeline-ship system has been established.
0 INTRODUCTION
Given the gradual reduction of terrestrial resources, the rich mineral resources of the deep seafloor have excited great interest and are likely to be commercially exploited in the near future. An integrated system with a seafloor-tracked mining vehicle, flexible hoses, a buffer, lifting pipes and a surface mining vessel is considered to be the most promising continuous mining system [1], as shown in Fig. 1.
Fig. 1. Schematic diagram of a 1000 m deep ocean mining system
To analyse the mobility of tracked mining vehicles on the seafloor, Dai et al. [2] and [3] built a seafloor-tracked mining vehicle as a single-body vehicle model with mesh elements of track-sediment interactions, as well as a multi-body dynamic model, and performed various operational simulations. Additionally, research has been conducted for modelling tracked vehicles on land surface soils [4] and [5]. Rubinstein and Coppock [5] used the program LMS-DADS to develop a multi-body model of tracked vehicles and applied user-defined forces to represent the trackterrain interactions. Gao et al. [6] and Wong et al. [7] and [8] developed two simulation programs, RTVPM and NTVPM, for parametrically modelling, analysing and evaluating tracked vehicles with long-pitch tracks and flexible tracks.
To evaluate the dynamics of the pipeline as well as the total ocean mining system, Hong et al. [9] developed a method for 3D nonlinear dynamic analysis of marine pipes. Rustad et al. [10] used the finite element method to build a mathematical model of an ocean pipeline system. Brink and Chung [11] carried out a computer simulation for the dynamic positioning control of a large ocean mining ship-pipe
system. Hong and Kim [12] performed analysis on the coupled dynamics of a tracked vehicle and flexible riser.
Thus far, however, the complex dynamic characteristics of the total deep ocean mining system have not yet been ascertained; therefore, its design scheme and continuous operation performance have yet to be evaluated, as it is costly and extremely difficult to perform in-situ ocean tests. The dynamic simulation is an effective solution. Accurate and efficient dynamic models for subsystems and then integrated into the total system have been developed in this paper.
1 TERRAMECHANICS MODEL OF THE SEAFLOOR SEDIMENT
As seafloor in-situ terramechanics tests are extremely difficult to perform, laboratory-simulated experiments were conducted. Fig. 2 shows the total experimental system.
A series of pressure-sinkage and shear experiments was performed, and corresponding experimental curves were obtained [3], as shown in Fig. 3.
The pressure-sinkage experimental curves can be characterized by Bekker's exponential equation, as follows [13]:
Simulated shear plate
Fig. 2. Photo and illustration of the terramechanics experimental system
b)
Fig. 3. Experimental curves and fitting curves: a) pressure-sinkage; b) shear stress- displacement
a)
p = ly+k
(1)
where p is the normal pressure; b is the width of the track; kc and k^, are the soil cohesive and friction modulus, respectively; z is the sinkage; and n is the soil deformation exponent.
The values kc, k^ and n in the Eq. (1) can be derived from the results of a minimum of two tests with two sizes of penetration plates, as illustrated in Fig. 3a. The tests produced two curves:
(2)
P1=( kc/ b1 + ^ P2 =( kcl b2+k$) z1 On the logarithmic scale, Eq. (2) can be rewritten
as:
[log pi = log (kjb + k ) + n log z
[log p2 = log [kj b2 + k ) + n log z
(3)
They represent two parallel straight lines with nearly same slopes on the log-log scale. The exponent of deformation n can be determined from the slopes of the straight lines. At the sinkage z = 1, the value of the normal pressure for the two sizes of plates are:
[(Pi L = hi bi + = K
l(P2 )z=i = kc/b2 + K2
(4)
In Eq. (4), b1 and b2 are known values as the widths of the penetration plates; K1 and K are measured values; and the only unknowns are kc and k^. Thus, two different sizes of plates were used in the pressure-sinkage experiments, numbered as penetration plates A and B, as illustrated in Fig. 3a.
On the log-log scale, the curves in Fig. 3a were expressed in Fig. 4:
The slopes of line A and line B in Fig. 4 were ca. 0.54 and 0.58. As mentioned above, the exponent of deformation n can be determined from the average value of the slopes of these two lines as n = 0.56.
It was also determined from Fig. 4 that when the sinkage z = 1 mm, the coordinate values of intersection points of line A and line B with the pressure axis were about 2.11 and 1.65, respectively. Converting the unit of the sinkage from millimetre to metre, while considering the exponents of deformation were 0.54 and 0.58, the following equations can be derived:
\kjbi + 87.96 (kPa/mn+2) |kc/b2 + k,= 90.67 (kPa/mn+2)'
(5)
10
Sinkage [m]
Fig. 4. Pressure stress-sinkage experiments curves on the logarithmic coordinate
Substituting the known values b1 = 0.15, b2 = 0.12 into Eq. (5), then kc =1.62 kN/mn+1, k^ = 77.1 kN/mn+2 can be obtained.
The observed shear property is different from most of the land surface soils [14]. However, Wong's equation can be fitted well [15]:
1
Kr (1 - O
-1
J-PK
(1 - -j(6)
where rmax is the maximum shear stress, Kr is the ratio of the residual shear stress to Tmax, and Km is the shear displacement when Tmax occurs. Observed from the shear experimental curves in Fig. 3b, the above parameters can be obtained directly.
2 MULTI-BODY DYNAMIC MODEL
OF A SEAFLOOR MINING VEHICLE
A multi-body dynamic simulation model of a seafloor-tracked mining vehicle has been established with the RecurDyn simulation program [3], as shown in Fig. 5.
Fig. 5. 3D multi-body dynamic simulation model of a seafloor-
tracked mining vehicle
A user-defined subroutine for characterizing the specific terramechanics models of the seafloor
T=Tmax Kr +
sediment was developed using the C language in the Visual Studio.Net environment and then integrated into the RecurDyn/Track simulation environment.
3 MULTI-BODY DEM MODEL OF A MINING PIPELINE
The pilot ocean mining pipeline is composed of a 900 m long rigid lifting pipe and a 400 m long flexible riser. Its nonlinear dynamic characteristics, such as large displacement and small deformation, should be taken into account. There were three main modelling methods for dynamic analysis of a pipeline: finite element method (FEM), lumped parameter method (LPM) and discrete element method (DEM). For the FEM, its computation is relatively large and computational time is extremely long. The LPM has a clear shortcoming in simulating the anti-bending and anti-torque performances of a pipeline. However, the DEM can describe well the large displacement and large rotation of a pipeline, and it can significantly reduce the computation as its stiffness matrix is a highly sparse matrix.
The long mining pipeline was discretized into a certain number of rigid elements linked by flexible connectors. The flexible connector was expressed by six spring elements. The external hydrodynamic loads were applied to the corresponding discrete rigid elements.
The stiffness coefficients for the flexible connector were the axial spring stiffness kx, lateral spring stiffness ky and kz, torsional spring stiffness krx, and rotational spring stiffness kry and krz, which can be expressed as:
\kx =EA/L; ky =12EIz/L3; kz =l2EIy/L krx=GID/L; kry =EIy/L; krz = EIz/L
(7)
y
where E and G are Young modulus and shear modulus; A and L are cross-section area and length; Iy and Iz are second moments of area; Ip is the polar moment of inertia.
The external hydrodynamic force can be calculated as:
1 I Tid1
Fh = - 1 Cd Pwd\vc - vp\(vc - vp )Cm Pw (vc - 2,9. Rabo primarne energije lahko zmanjšamo z uporabo toplote iz kogeneracijskih sistemov (proizvodnja električne energije in toplote), ko vso toploto izrabimo za pogon absorpcijskih hladilnikov. Predčasen odjem toplote iz odjemno-kondenzacijske turbine zmanjšuje električni izkoristek termo elektrarne. Zaradi višjega tlaka in temperature odjemne pare iz turbine v primeru uporabe dvostopenjskega absorpcijskega hladilnika je električni izkoristek nižji, kot v primeru enostopenjskega absorpcijskega hladilnika. Razmerje x nam pove, kakšno je razmerje med električnim izkoristkom v primeru odjema pare iz odjemno-kondenzacijske turbine (dva produkta: električna energija in toplota) ter električnim izkoristkom kondenzacijske turbine (para v turbini popolnoma ekspandira, imamo samo en produkt-električno energijo).
Splošni matematični model za izračun faktorja primarne energije sistema daljinskega hlajenja je uporaben za določitev energijsko najučinkovitejše tehnologije proizvodnje hladu v smislu minimalne rabe primarne energije.
V prihodnosti bi bilo potrebno ta model razširiti za primere sočasne proizvodnje in rabe električne energije, toplote za ogrevanje in hladu za hlajenje. V tem primeru je potrebno upoštevati razmerje med toploto, ki jo pozimi uporabimo za ogrevanje ter toploto, ki jo poleti uporabimo za pogon absorpcijskih hladilnikov. Potrebno je upoštevati tudi povečanje rabe primarne energije zaradi predčasnega odjema pare iz odjemno kondenzacijske turbine, pri čemer se raba primarne energije povečuje z večanjem tlaka in temperature odjemne pare.
Ključne besede: faktor primarne energije, proizvodnja hladu, daljinsko hlajenje, hladilnik
Integracija planetnega gonila in prostorskega mehanizma RCRCR v novem univerzalnem reduktorju
Cheng Gu12 - Xinbo Chen12* 1 Univerza Tongji, Center za čisto energijo v avtomobilski tehniki, Kitajska 2 Univerza Tongji, Šola za avtomobilske študije, Kitajska
Reduktorji, ki opravljajo naloge zmanjševanja vrtilne hitrosti in pretvorbe vrtilnega momenta, so pomemben podsestav sodobnih strojev. Težava pri klasičnih reduktorjih je v tem, da se ne prilagajajo spremembam kota med vhodno in izhodno gredjo. Zato se običajno vgradi par kardanskih zgibov, ki pa spremenijo razmerja med vrtenji gredi in povečajo prenosni sistem v aksialni smeri.
V članku je predstavljen predlog novega mehanizma reduktorja, ki se lahko prilagaja spremembam kota med gredema. Sestavljen je iz planetnega gonila z majhno razliko v številu zob in prostorskega mehanizma RCRCR (R - tečajni par, C - cilindrični par). Mehanizmu so dodani še trije pari povezovalnih gredi za podporo in bolj homogeno porazdelitev obremenitev. Prestavno razmerje značilnega planetnega reduktorja in število prostostnih stopenj mehanizma RCRCR sta izračunana po teoriji mehanizmov.
Opravljena je kinematična analiza prostorskega mehanizma RCRCR, pri čemer je razmerje med koordinatnimi sistemi izraženo z neposredno kosinusno matriko in načelom enakovredne geometrije. V programski opremi za samodejno dinamično analizo mehanskih sistemov (ADAMS) je postavljen model mehanizma za simulacijo. Kotne hitrosti vhodne gredi, izhodne gredi in planetnega gonila v reduktorju so izračunane za kota a = 135° in a = 160° v simulacijah z enakimi parametri, vključno s simuliranim časom in vhodno hitrostjo.
Končno je zasnovan tudi prototip po teoriji mehanizmov. Ugotovljeno je, da je mehanizem RCRCR prostorski mehanizem z eno prostostno stopnjo, ki zagotavlja zahtevano translacijsko gibanje planetnega gonila v univerzalnem reduktorju. Prestavno razmerje reduktorja ni odvisno od kota ali od razdalje med gredema in gre torej za univerzalen reduktor, ki ohranja stabilno prestavno razmerje in se lahko prilagaja različnim kotom. Kinematični pari mehanizma so v površinskem stiku, kar je ugodno za njihove protiobrabne lastnosti.
Pripravljeni so podrobni načrti za proizvodnjo. Za gladek in zanesljiv prenos je zelo pomembna natančnost proizvodnje in montaže. Opravljena je začetna trdnostna analiza glavnih komponent in podane so smernice za proizvodnjo in mazanje.
Delo še ni dokončano in kliče po nekaterih izboljšavah. Preizkusiti bi bilo treba mehansko učinkovitost predstavljenega univerzalnega reduktorja, za večjo uporabnost pa bi bila potrebna tudi lahka konstrukcija in trdnostna analiza za večjo kompaktnost in integracijo konstrukcije. Edinstvena možnost prilagajanja spremembam kota med gredema odpira priložnosti za uporabo v reduktorjih, ki obratujejo v težavnih pogojih, možnost kompenzacije relativnih odmikov pri transmisijskih gredeh pa je zaželena tako pri montaži kakor tudi pri obratovanju.
Ključne besede: integracija, samoprilagajanje, planetno gonilo, razlika nekaj zob, prostorski mehanizem RCRCR, univerzalni reduktor
SI 124
*Naslov avtorja za dopisovanje: Univerza Tongji, Center za čisto energijo v avtomobilski tehniki, No.4800 Cao'an highway, Shanghai, Kitajska, xinbo_chen@126.com
Kombinirani model za ocenjevanje stanja strojev na osnovi tehnike Isomap in globoke verjetnostne mreže
Aijun Yin1 - Juncheng Lu1* - Zongxian Dai2 - Jiang Li1 - Ouyang Qi3
1 Chongqing univerza, Kolidž za strojništvo, Državni laboratorij za mehanske prenose, Kitajska 2 Chongqing akademija za meroslovje in kakovost, Kitajska 3 Chongqing univerza, Kolidž za avtomatizacijo, Kitajska
Nadzor in vrednotenje trendov degradacije nekaterih ključnih strojnih delov kot so ležaji omogoča odpravo degradirane zmogljivosti ali napak še pred okvaro stroja. Ker pa količine zbranih podatkov o strojih postajajo vse obilnejše in ker so vse strožje tudi zahteve glede hitrosti in točnosti vrednotenja stanja strojev, tradicionalne metode ne jamčijo več za učinkovito delo.
Isomap kot tehnika za nelinearno globalno transformacijo dimenzionalnosti poišče rešitev preslikave z vrsto pretvorb, ki omogočijo predstavitev geodezične razdalje med vhodnimi točkami v izvirnem prostoru z evklidsko razdaljo v prostoru projekcije. Globoka verjetnostna mreža (DBN) kot probabilistični generativni model, ki uspešno zajema značilne informacije v surovih podatkih z raznimi nelinearnimi transformacijami in kompleksnimi aproksimativnimi nelinearnimi funkcijami, je primerno orodje za vrednotenje stanja strojev. Po primerjavi in analizi pomanjkljivosti obstoječih metod je v članku predstavljen novi kombinirani model vrednotenja (CAM), ki združuje WPT, Isomap in DBN za vrednotenje stanja obravnavanega stroja (oz. njegovih kotalnih ležajev).
Podatki o pospešeni degradaciji ležajev so bili zbrani s preizkušanjem ležajev pod stalno obremenitvijo do odpovedi na posebnem preizkuševališču, ki so ga izdelali na Univerzi v Cincinnatiju. Zaradi nestacionarnosti signalov vibracij, zbranih na strojih, in za odpravo naključnih dejavnikov so bile uporabljene klasične (analiza v časovni in frekvenčni domeni) in sodobnejše metode za določanje značilnosti (paketna valčna transformacija (WPT)). Po izboru je bilo izluščenih 11 originalnih značilnosti v časovni domeni (kot sta temenski faktor in vršna vrednost), 13 originalnih značilnosti kot sta standardna deviacija frekvence in koeficient sploščenosti, ter 14 originalnih značilnosti WPT. Nato je bila uporabljena tehnika Isomap za pretvorbo visokodimenzijskega originalnega polja značilnosti 9810*38 v nizkodimenzijsko polje 9810*6, ki vsebuje več reprezentančnih značilnosti. Sledilo je učenje postavljenih modelov DBN z RBM po algoritmu kontrastivne divergence (CD). Po vrsti eksperimentov se je pokazalo, da je najboljši delovni model za vrednotenje stanja stroja DBN s petimi skritimi nivoji in arhitekturo 100-100-50-50-50. Pridobljeno šestdimenzijsko polje značilnosti je bilo vneseno v naučeno DBN za vrednotenje delovanja ležajev. Za celovito primerjavo in osvetlitev prednosti predlagane metode sta bili za izvedbo istih eksperimentov uvedeni tudi dve drugi priljubljeni metodi za zmanjšanje dimenzionalnosti (PCA in Laplaceove lastne mape) in dva pametna algoritma (HMM in nevronska mreža z vzvratnim postopkom učenja (BPNN)). Primerjalni eksperimenti so dokazali, da je predlagani model CAM občutljivejši na začetne napake v kotalnih ležajih in učinkovitejši pri vrednotenju degradacije ležajev. Zaznal je prve znake degradacije v točki 530, razpoke in utrujenostno luščenje ter nekatere druge resne napake v točki 700, ter popolno odpoved po točki 940. To je v skladu z dejanskimi pogoji obratovanja ležajev.
CAM bi bilo mogoče v prihodnjih raziskavah uporabiti tudi pri vrednotenju ali klasifikaciji na drugih področjih, med drugim tudi za določevanje občutljivejših značilnosti vibracijskih signalov in za optimizacijo modela vrednotenja.
Ključne besede: Isomap, zmanjšanje dimenzionalnosti, globoko verjetnostno omrežje (DBN), stanje stroja, kombinirani model vrednotenja (CAM)
O utrujanju jeklenih visečih dvižnih cevi
Nnamdi Onochie Chibueze, Chinwuba Victor Ossia*, John Umunna Okoli
University v Port Harcourtu, Tehnološki inštitut, Nigeria
Obstoječa črpališča se počasi praznijo in iskanje nafte ter plina se zato seli v globoke vode in v zahtevnejša okolja. Vse večje potrebe po energiji tudi spodbujajo razvoj naprednejših tehnologij. Jeklene viseče dvižne cevi (SCR) so se uveljavile kot rešitev za globoka in ultragloboka črpališča predvsem zaradi svoje ekonomičnosti, obstojnosti proti tlaku, tehnično preproste montaže in možnosti priključitve na različne gostiteljske ploščadi. Zato obstaja potreba po raziskavah utrujenostne trajnostne dobe SCR v Gvinejskem zalivu, kjer uvajajo te cevi na osnovi zbranih podatkov iz Severnega morja in iz Mehiškega zaliva.
Jeklene viseče povezovalne cevi so občutljive na dinamiko gostiteljske ploščadi in so podvržene utrujanju. Uporaba SCR pri polpotopljivih ploščadih v globokih vodah prinaša posebne konstrukcijske izzive zaradi velikih amplitud gibanja ploščadi pod vplivom valov in velikih odmikov ploščadi, ki jih povzročajo vibracije zaradi vrtincev (VIV), gibanja zaradi vrtincev (VIM), veter in počasni tokovi. Močno zibanje ploščadi tako povzroča uklon in utrujanje v območju stika dvižne cevi z morskim dnom (TDP). Pri načrtovanju SCR je izjemno pomembna izbira ustreznega spektra valovanja.
Homogena, 2200 m dolga cev SCR, ki se začne 1000 m nad morskim dnom, je bila modelirana v programskem paketu za analize po metodi končnih elementov Orcaflex s 50 m dolgim gibkim spojem, 350 m dolgo oplato, 1300 m dolgo dvižno cevjo in 500 m dolgim pretočnim vodom. Konec A je pritrjen na FPSO, konec B pa je zasidran v morsko dno. Na oddaljenosti 900 m je na cev SCR pritrjena boja 6D z maso 104 kg, dolžino 6 m in prostornino 20 m3. Pri analizi utrujanja so bili uporabljeni oceanografski podatki, podatki o vetru, tokovih, morski vodi in morskem dnu. Za primerjavo sta bila uporabljena spektra valovanja Ochi-Hubble in JONSWAP. Pred analizo utrujanja je bil pripravljen set simulacijskih datotek OrcaFlex za različne realne obremenitvene pogoje. Temu obremenitvenemu primeru je bila dodeljena 20-letna izpostavitev, ki ustreza dobi trajanja dvižne cevi. Simulacija po 16-sekundnem preizkusu je trajala 3 ure.
Podana je primerjava rezultatov utrujanja SCR za spektra valovanja JONSWAP in Ochi-Hubble. Analiza utrujenostne trajnostne dobe po dolžini cevi na osnovi spektra valovanja JONWAP je pokazala, da je utrujenostna trajnostna doba v območju obešenja (HOFL) 10 let, območje utrujenostne trajnostne dobe (FL) od 500 do 5000 let, utrujenostna trajnostna doba v območju TDP 1,8 leta na globini 1828,8 m in najdaljša utrujenostna trajnostna doba za TDP (TDP FLmax) 10 000 let. Pri spektru valovanja Ochi-Hubble pa je bilo izračunano, da je utrujenostna trajnostna doba v območju obešenja (HOFL) 7 let, območje utrujenostne trajnostne dobe (FL) od 100 do 5000 let, utrujenostna trajnostna doba v območju TDP 13,6 leta na globini 1828,8 m in najdaljša utrujenostna trajnostna doba za TDP (TDP FLmax) 10 000 let.
Utrujenostna trajnostna doba v območju TDP ter v območju obešenja pri spektrih JONSWAP in Ochi-Hubble je zelo kratka in izhaja iz večjih napetostih, ki jih povzročajo dinamika FPSO, oceanski tokovi, spektri valovanja in druge obremenitve iz okolja.
Rezultata za utrujenostno dobo SCR v točki dotika z morskim dnom (TDP) po spektrih Ochi-Hubble in JONSWAP sta 13,6 in 1,8 leta, torej se razlikujeta za faktor 7,5. Iz krivulj S-N za oba spektra valovanja je razvidna inverzna povezava med napetostjo S in številom ciklov do odpovedi N na log-log skali. Daljša utrujenostna trajnostna doba SCR pri spektru valovanja Ochi-Hubble v primerjavi s spektrom JONSWAP, ki ga pogosto uporabljajo projektanti, pomeni, da spekter Ochi-Hubble dobro opisuje stanje v Gvinejskem zalivu in ne privede do predimenzioniranja. Uporaba spektra JONSWAP pa je nasprotno povezana s predimenzioniranjem in posledično z večjimi stroški.
Kratka utrujenostna trajnostna doba v območju TDP in v območju obešenja SCR izhaja iz večjih napetosti, ki jih povzročajo dinamika FPSO, oceanski tokovi, spektri valovanja in druge obremenitve iz okolja. Efektivna napetost v SCR zato hitro zraste od nič v območju obešenja do maksimalne vrednosti 30*106 kN na dolžini 5 m. Prispevek te študije so tudi obogateni oceanografski podatki za Gvinejski zaliv.
Ključne besede: jeklene viseče dvižne cevi, točka stika z morskim dnom, Ochi-Hubble, JONSWAP, Gvinejski zaliv, utrujanje
SI 126
*Naslov avtorja za dopisovanje: University v Port Harcourtu, Tehnološki inštitut, Port Harcourt, Nigerija; ossiacv@otiuniport.org
Nov model dinamike sistema več teles globokomorsko rudarsko vozilo-cevovod-ladja in integrirana simulacija njegovih gibanj
Yu Dal1-2-3-* - Liping Pang1 - Lisong Chen1 - Xiang Zhu1 - Tao Zhang1 1 Univerza srednjega juga, Kolidž za strojništvo in elektrotehniko, Kitajska
2 Univerza srednjega juga, Raziskovalni institut Shenzhen, Kitajska 3 Univerza v Šanghaju, Državni laboratorij za ladijsko strojništvo, Kitajska
Do danes še ni bil narejen opis kompleksne dinamike celotnih sistemov za rudarjenje na globokih morjih. Izvajanje preizkusov na oceanu bi bila draga in izjemno težavna naloga, zato s tem pristopom ni mogoče vrednotiti konstrukcijskih rešitev in njihove sposobnosti za trajno obratovanje. Učinkovit vpogled v te lastnosti pa je mogoče pridobiti s simulacijami dinamike sistema.
Članek opisuje razvoj točnih in učinkovitih modelov dinamike podsistemov, ki so združeni v model celotnega sistema. V jeziku C je bila pripravljena inovativna uporabniška podrutina za popis posebnega teramehanskega modela sedimentov na morskem dnu, ki je bila dodana v paket za simulacije RecurDyn/Track in nato validirana z laboratorijskimi preizkusi.
V paketu RecurDyn/Track je bil z integracijo pripravljenega teramehanskega modela sedimentov na morskem dnu razvit tridimenzionalni model dinamike več teles za vozilo za rudarjenje na morskem dnu. Metoda diskretnih elementov (DEM) je v primerjavi s tradicionalno metodo končnih elementov (MKE) in metodo združenih parametrov (LPM) učinkovitejše in primernejše sredstvo za hitro simulacijo dinamike pri zelo dolgih cevovodih. V modelu DEM cevovoda za rudarjenje na globokem morju so bili tako uporabljeni diskretni togi elementi, povezani z upogljivimi deli.
Brezmasni upogljivi povezovalni deli so popisani s šestimi vzmetnimi elementi, zunanje hidrodinamične obremenitve pa delujejo v masnih središčih teh elementov. Model vključuje veliko število diskretnih togih elementov, diskretne hidrodinamične sile in upogljive povezovalne dele, zato bi bilo zelo zamudno in povezano z napakami, če bi ga gradili na klasičen način, torej za vsak diskretni element posebej. Namesto tega je bila s pomočjo sekundarne razvojne platforme ProcessNet razvita inovativna uporabniška podrutina v jeziku C# za hitro in samodejno parametrično izgradnjo modela DEM v programskem paketu RecurDyn. Zgrajen je bil kinematični model ladje za rudarjenje s šestimi prostostnimi stopnjami, pri katerem je bil ključen popis zibanja.
Z združitvijo modelov podsistemov je bil postavljen nov tridimenzionalni model dinamike več teles za celotni sistem globokomorsko vozilo za rudarjenje-cevovod-ladja. Z integrirano simulacijo gibanj so bile pridobljene nekatere ključne lastnosti dinamike, kot so trajektorije gibanj, interakcije med podsistemi itd. Izkazalo se je, da je stanje gibanja vseh podsistemov stabilno znotraj dovoljenih meja ter da je mogoče zagotoviti sinhronizirano stabilno gibanje celotnega sistema med obratovanjem.
Predstavljena raziskava tako podaja dragoceno in učinkovito metodo modeliranja za namene optimizacije konstrukcije, vrednotenja zmogljivosti in integriranega nadzora delovanja celotnega sistema vozilo za rudarjenje na globokem morju-cevovod-ladja.
V prihodnjih raziskavah bo postavljen laboratorijski sistem za eksperimentalno preučevanje integriranega sistema za rudarjenje na morju. Pri tem bodo uporabljeni načelo podobnosti ter načrtovalska merila za preizkušanje integriranih gibanj in operativnega nadzora celotnega sistema. Na ta način bo omogočena celovita validacija rezultatov simulacij dinamike iz tega članka ter priprava predlogov za optimizacijo zasnove in praktične metode upravljanja sistemov za rudarjenje na globokih morjih.
Ključne besede: sistem za rudarjenje na globokih morjih, teramehanski model, model dinamike več teles, model z diskretnimi elementi, uporabniška podrutina, integrirana simulacija gibanja
Strojniški vestnik - Journal of Mechanical Engineering 62(2016)12, SI 128-129 Osebne objave
DOKTORSKE DISERTACIJE
Na Fakulteti za strojništvo Univerze v Ljubljani so obranili svojo doktorsko disertacijo: • dne 25. novembra 2016 Andrej ŠARC z naslovom: »Odstranjevanje bakterije Legionella pneumophila s pomočjo hidrodinamske kavitacije« (mentor: izr. prof. dr. Matevž Dular, somentor: prof. dr. David Stopar);
Visoka koncentracija bakterije Legionella pneumophila v vodovodnih inštalacijah predstavlja resno nevarnost zdravju človeka. Do sedaj ni poznanih učinkovitih mehanizmov, ki bi konstantno zagotavljali nizke koncentracije bakterij L. pneumophila v vodi. Trenutni ukrepi za preprečevanje rasti bakterij in preprečevanje možnosti izbruha so dragi in okoljsko nesprejemljivi ter v večini primerov ne zagotovijo varnosti pred morebitnimi ponovnimi izbruhi.
V okviru doktorske disertacije smo za odstranjevanje bakterije L. pneumophila preučili možnost uporabe kavitacije - fizikalnega pojava, ki označuje prehod iz kapljevine v plinasto fazo in nazaj v homogeno kapljevino. Razlog za nastanek kavitacije je krajevno zmanjšanje tlaka, kjer temperatura medija ostane približno konstantna. Ob ponovnem prehodu v kapljevito stanje v območju povečanega tlaka prihaja do velikih tlačnih in temperaturnih nihanj, ki so po eni strani nezaželeni (poškodbe naprav) po drugi strani pa so lahko uporabni pri različnih procesih.
Z meritvami rasti mikroorganizmov smo določili vplive različnih režimov kavitacijskega toka na preživljivost bakterije L. pneumophila. Pokazali smo, da različni tipi kavitacije omogočajo različne pogoje,
ki vodijo do uničenja bakterij, ter da je učinkovitost delovanja odvisna od tipa kavitacije in obratovalnih pogojev modelne naprave za uničevanje bakterij.
Rezultati dela so privedli do boljšega poznavanja vpliva kavitacije na bakterije v vodi in do razvoja naprave za generiranje kavitacije in njenih potencialnih aplikacij v vodovodnih sistemih; • dne 30. novembra 2016 Matija BRUMAT z naslovom: »Prostorsko dušenje prožnih struktur« (mentor: izr. prof. dr. Janko Slavič);
Doktorsko delo obravnava prostorsko dušenje prožnih struktur. Opravljen je bil pregled 22 obstoječih metod za identifikacijo prostorskega dušenja. Za poglobljeno eksperimentalno vrednotenje metod identifikacije v realnih pogojih se je razvilo namensko preskuševališče, ki omogoča analizo in nadzor nad vplivnimi parametri (modalni in prostorski odrez, razlike v identifikaciji med histereznim in viskoznim modelom dušenja in na razlike ob prisotnosti dodatne dušilke na strukturi). Za eksperimentalno analizo so bile izbrane štiri najobetavnejše metode: Lee-Kim, Chen-Ju-Tsuei, Fritzen IV in metoda lokalne enačbe gibanja, kjer se je ugotovilo občuten vpliv modalnega in prostorskega odreza na identificirano prostorsko matriko dušenja, ter da je pri znani strukturi bolje uporabiti metodo lokalne enačbe gibanja. V zadnjem delu naloge je bil vpeljan pristop za načrtovanje prostorskega dušenja prožne strukture na podlagi želenega frekvenčnega odziva, ki je validiran z eksperimentom na nosilcu.
v
V spomin profesorju dr. Iztoku Zunu
Vest o smrti profesorja dr. Iztoka Žuna nas je močno presenetila in globoko pretresla. Svojo bolezen je večkrat premagal, zato nismo dvomili, da bo tudi tokrat tako.
Profesor dr. Iztok Žun je diplomiral na Fakulteti za strojništvo v Ljubljani leta 1972. V želji po iskanju novega znanja je nadaljeval študij. Magister znanosti je postal leta 1974, nato pa je bil na študijskem obisku na Thayer School of Engineering na Dartmouth College v ZDA. Doktorsko disertacijo tehniških znanosti je obranil leta 1976. Njegova učiteljska pot se je pričela leta 1977, ko je na Fakulteti za strojništvo postal docent za termodinamiko in mehaniko fluidov. Leta 1982 je bil izvoljen v izrednega profesorja, leta 1987 pa v rednega profesorja. V svoji dolgoletni karieri je uspešno prenašal svoje znanje na mlade, saj se je vedno trudil študentom vzbuditi občutek radovednosti, da bi lahko premikali meje svojega razumevanja.
Primarni raziskovalni interes profesorja Žuna je bil dvofazni tok, pri čemer seje posebej posvečal mehurčastemu toku. Vedel je, da se izvirno gonilo tehničnega napredka prenosa toplote in snovi nahaja v razvoju razumevanja fizike večfaznih tokov. Videl je pomen teh raziskav za širok nabor industrijskih panog. Spoznanja profesorja Žuna so odmevna v mednarodnem merilu, saj je nekatere probleme s področja kompleksnih tokov rešil prvi v svetu. Prvi je pokazal, kako polje prečne vzgonske sile deluje na mehurčke v turbulentnem toku kapljevine in povzroča nehomogeno porazdelitev deleža plinske faze v vertikalnem kanalu. Odkritje je bilo pomembno za pojasnitev nehomogene porazdelitve parne oz. plinaste faze v jedrskih in kemijskih reaktorjih. Bil je prvi, ki je predstavil rezultate numerične simulacije evolucije kopičenja deleža plinske faze iz območja stene kanala v njegov osrednji del, kar je predstavljalo enega izmed glavnih problemov v analizah jedrskih reaktorjev. Nadaljne študije modeliranja na večih skalah hkrati v sodelovanju z dr. Tomiyamo iz Kobe Univerze na Japonskem so razkrile vzroke bočnih migracij mehurčkov, kar je nudilo trdno osnovo relevantnim računalniškim kodam, ki se danes uporabljajo. Znanstveno-raziskovalni opus profesorja Žuna obsega več kot 50 člankov v revijah in več kot 130 člankov v recenziranih zbornikih. Ob tem je imel več kot 80 vabljenih predavanj na konferencah, inštitutih, univerzah in v industriji v EU, na Japonskem in v ZDA. Ima več kot 1000 citatov ter 5 patentov.
Vzporedno z raziskavami dinamike dvofaznega toka je profesor Žun svoje znanje prenašal v proizvodno okolje. Leta 1989 je uspešno realiziral implementacijo nadzora, upravljanja in modeliranja slovenskega plinovodnega sistema. Učinkovito upravljanje sistema je bilo izvedeno z uporabo koncepta dinamične eksergije, ki jo je profesor Žun uvedel kot novost. Med pomembnejše aplikacije v industriji sodita še: projekt usmerjen v zmanjšanje emisij dizelskega motorja ter razvoj procesne naprave za oblaganje delcev v farmacevtskih tehnologijah.
Profesor Žun je leta 1983 ustanovil in od takrat vodil Laboratorij za dinamiko fluidov in termodinamiko. Leta 1997 je ustanovil katedro z istim imenom ter bil ves čas njen predstojnik. Bil je eden od glavnih pobudnikov in zagovornikov koncepta katedre kot osnovne izobraževalne, raziskovalne in razvojne enote, pri čemer je predstojnik odgovoren za njen razvoj. Od mnogih vodstvenih in drugih pomembnih funkcij doma in v svetu naj omenim le nekatere. Leta 1998 je bil ustanovitelj ter do 2015 predsednik Evropsko-japonskega združenja za dvofazne tokove. Leta 1999 je ustanovil in od takrat vodil nacionalni raziskovalni program Tranzientni dvofazni tokovi. Od leta 2008 je bil redni član Inženirske akademije Slovenije. Bil je eden od njenih ustanoviteljev, nekaj let tudi aktiven član njenega Izvršnega odbora. Od leta 2014 je bil podpredsednik Virtualnega mednarodnega razvojnega inštituta za dvofazni tok in prenos toplote.
Od nagrad, ki jih je prejel za svoj doprinos k znanosti in tehnološkemu razvoju, je potrebno izpostaviti Zoisovo priznanje leta 1999 in Nagrado Japonskega združenja za večfazne tokove leta 2003. To združenje mu je letos podelilo plaketo o častnem članstvu za izjemen prispevek k znanosti in tehnologiji na področju večfaznih tokov.
Pogrešali bomo njegovo neizmerno voljo do raziskovanja ter navdušenje nad iskanjem novih spoznanj. Njegova zapuščina je velika, z njo je ustvaril pogoje za nadaljevanje kreativnega raziskovalnega, razvojnega in izobraževalnega dela za svoje bližnje sodelavce in za Fakulteto za strojništvo v celoti.
Ohranili ga bomo v trajnem spominu.
doc. dr. Matjaž Perpar
Information for Authors
All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 10 pages. Longer contributions will only be accepted if authors provide juastification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu .
SUBMISSION:
Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process.
The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided.
SUBMISSION CONTENT:
The typical submission material consists of:
- A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references),
- Supplementary files:
• a manuscript in a WORD file format
• a cover letter (please see instructions for composing the cover letter)
• a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files)
• possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory
comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments.
COVER LETTER INSTRUCTIONS:
Please add a cover letter stating the following information about the submitted paper:
1. Paper title, list of authors and their affiliations.
2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03).
3. A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere.
4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose.
5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors.
FORMAT OF THE MANUSCRIPT:
The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format:
- A Title that adequately describes the content of the manuscript.
- A list of Authors and their affiliations.
- An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions.
- 4 to 6 significant key words should follow the abstract to aid indexing.
- 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it.
- An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated.
- A Methods section detailing the theoretical or experimental methods used.
- An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results.
- A Results section that should clearly and concisely present the data, using figures and tables where appropriate.
- A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.)
- A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract.
- Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research.
- Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature.
- References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript.
- Appendix(-icies) if any.
SPECIAL NOTES
Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf .
Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References.
Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing.
Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic.
REFERENCES:
A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary.
References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow.
Journal Papers:
Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code.
[1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnih - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, DOI:10.5545/sv-jme.2011.013.
Journal titles should not be abbreviated. Note that journal title is set in italics. Books:
Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication.
[2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken.
Note that the title of the book is italicized.
Chapters in Books:
Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages.
[3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordic, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576.
Proceedings Papers:
Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages.
[4] Štefanic, N., Martinčevic-Mikic, S., Tošanovic, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427.
Standards:
Standard-Code (year). Title. Organisation. Place.
[5] ISO/DIS 16000-6.2:2002. Indoor Air — Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva.
WWW pages:
Surname, Initials or Company name. Title, from http://address, date of access.
[6] Rockwell Automation. Arena, from http://www.arenastmuiation.com, accessed on 200909-07.
EXTENDED ABSTRACT:
When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters). The instruction for composing the extended abstract are published on-line: http://www.sv-ime.eu/mformation-for-authors/ .
COPYRIGHT:
Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/.
PUBLICATION FEE:
Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 240.00 EUR (for articles with maximum of 6 pages), 300.00 EUR (for articles with maximum of 10 pages), plus 30.00 EUR for each additional page. The additional cost for a color page is 90.00 EUR. These fees do not include tax.
Strojniški vestnik -Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu
http://www.sv-jme.eu
Contents
Papers
Jure Rejc, Marko Munih
Robust visual touch-up calibration method in robot laser spot welding application
Weiping Wang, Bo Wang:
An Energy-Saving Control Strategy with Load Sensing for Electro-Hydraulic Servo Systems
Tjasa Duh Coz, Andrej Kitanovski, Alojz Poredos: Primary Energy Factor of a District Cooling System
Cheng Gu, Xinbo Chen:
A Novel Universal Reducer Integrating a Planetary Gear Mechanism with an RCRCR Spatial Mechanism
Aijun Yin, Juncheng Lu, Zongxian Dai, Jiang Li, Qi Ouyang: Isomap and Deep Belief Network-Based Machine Health Combined Assessment Model
Nnamdi Onochie Chibueze, Chinwuba Victor Ossia, John Umunna OkoLi: On the Fatigue of Steel Catenary Risers
Yu Dai, Liping Pang, Lisong Chen, Xiang Zhu, Tao Zhang: A New Multi-Body Dynamic Model of a Deep Ocean Mining Vehicle-Pipeline-Ship System and Simulation of Its Integrated Motion
730
740
751
9770039248001