Abstract
This paper aimed to study the longitudinal vibration characteristics of the 5000 m mining pipe in the ocean under different working wind conditions, offset angle, damping, and ore bin weight. Based on the finite element method, the mining pipe is simplified into beam element and discretized, and the physical and mathematical models of the mining pipe system are established. The Wilson-θ direct integral method is adopted for numerical calculation. The results show that the longitudinal vibration of the mining pipe is irregular, which presents the phenomenon of oscillation. The vibration amplitude decreases first and then increases from top to bottom, the minimum vibration amplitude appears at 1000 m, and the maximum vibration amplitude appears at the top of the mining pipe. Under the same working wind condition, the overall longitudinal vibration amplitude of the mining pipe can be increased by increasing the ore bin weight and the offset angle, but neither of them can change the frequency of the longitudinal vibration. The closer the excitation frequency generated by different working wind conditions is to the natural frequency, the larger the mining pipe longitudinal vibration amplitude is. The closer the vibration frequency generated by the same excitation frequency is to the natural frequency, the stronger the vibration intensity is, and when damping is added, the vibration intensity decreases faster.
1. Introduction
With the rapid development of the global economy and industry, people’s demand for mineral resources is increasing. Due to the shortage of land resources and the deepening of people’s understanding of the ocean, people turn their attention to the ocean which is rich in mineral resources. The development of marine mineral resources has become a new strategic goal of all countries in the world [1]. Seabed mineral resources are rich and the environment is complex, and the exploitation of resources must be aided by the deep sea mining systems [2]. The deep sea mining system can collect mineral resources such as seabed polymetallic nodules, lift them, and transport them back to the land [3]. The deep sea mining system mainly includes mining vessel, lifting pump group, mining pipe, ore bin, flexible pipe, and mining machine [4]. As an important part of the deep sea mining system, the mining pipe system is affected by waves and currents in the deep sea, and a series of vibration deformation occurs [5]. The longitudinal vibration of the mining pipe is the most dangerous type of vibration, and it must be reduced as much as possible [6].
Korean scholars used the dynamic program OMPA to study the drag characteristics of the mining ship models under the conditions of 5000 m and 6000 m water depth based on the iterative method of concentrated mass [7]. Chung concentrated on the mechanical properties of the three-dimensional nonlinear, axial bending, and torsional coupling of the mining pipe [8]. Achouyab and Bahrar simplified the pipe into a plane beam and used Newmark algorithm to calculate the pipe transverse displacement [9]. Japanese research scholar Aso used the analytical method to study the dynamic characteristics of equal-diameter mining pipe under the variable speed towing conditions [10]. Cheng et al. regarded the mining pipe as a beam model and used three-dimensional beam element modeling to conduct dynamic analysis on the 1000 m mining pipe [11]. Chinese research scholar Tang established a vibration test system in the laboratory to test the vibration characteristics of the pipeline [12]. Dr. Liof Central South University established a spatial discrete model of the 1000 m mining pipe by using the multirigid body discrete element method and analyzed the overall dynamic simulation and performance of mining system based on multirigid body dynamic theory [13].
In summary, previous studies mainly focused on the dynamic response of 1000 m equal diameter mining pipe or the dynamic response of 5000 m equal diameter mining pipe under towing condition, and the longitudinal vibration characteristics of the 5000 m mining pipe were less studied. In this paper, the 5000 m ladder-shaped mining pipe is simplified into the nonlinear beam element. The finite element method is used for the dynamic theoretical analysis, and the algorithm is used for numerical calculation in MATLAB to obtain the longitudinal vibration response of the mining pipe under different damping, different ore bin quality, and offset angles. The main research contents can provide useful reference for the design of deep sea mining system, the research of heave compensation system, and the actual deep sea mining operations.
2. Physical Modeling
The physical structure of the deep sea mining system is shown in Figure 1 [14]. The system mainly consists of the mineral collection subsystem, mining pipe subsystem, and mining ship. The mineral collection subsystem includes mining machine and flexible pipe. The mining machine can automatically move on the seabed to collect the polymetallic nodules stored in the bottom of the sea and in the mining machine to clean desilting and crush polymetallic nodules. Then, the polymetallic nodules are transported to the middle ore bin through the flexible pipe. Under the action of the lifting pump set, the polymetallic nodules in the middle ore bin are transported to the surface mining ship through the mining pipe. Finally, the resources are transported back to the land.
The schematic diagram of the mining pipe subsystem physical structure is shown in Figure 2. The mining pipe subsystem mainly includes the mining pipe, lifting pump group, middle ore bin, shock absorber, etc. The gravity of all underwater equipment is transferred to the mining ship through the mining pipe; in order to reduce the load borne by the mining pipe, the mining pipe is in the shape of a ladder, which is divided into four step segments. The top end of the mining pipe is connected with the mining ship. The first mining pipe is installed with lifting pump group and vibration absorber device. The bottom end of the mining pipe is connected with the middle ore bin, and the ore bin is connected with the flexible pipe. The top of the pipe string is subjected to the harmonic heave motion of the mining ship.
3. Theoretical Analysis Method
The basic idea of the finite element method is to discretize a continuum into a finite set of unit combinations that are connected to each other in a certain way. Since the elements can be combined in different ways and the elements themselves have different shapes, the finite element method can solve the fields with complex geometric shapes [15]. The approximate function within a unit is usually expressed by the unknown field function or its derivative value at each node of the unit and its interpolation function, so a continuous infinite degree of freedom problem becomes a discrete finite degree of freedom problem.
3.1. Beam Element Analysis
Due to the wide application of beam structure [16], in the longitudinal motion, the mining pipe is simplified into a space beam unit and a small beam unit is taken in a three-dimensional space. As shown in Figure 3, let the nodes at both ends be i and j, and ,, represent the linear displacement component in three directions. Then, the unit node displacement vector is .
Hermite difference is used to obtain the interpolation function [17]:where , , , , , and are the interpolation functions of the beam element and l is the length of the beam unit.
Unit node displacement function:
The geometric equation corresponding to the beam unit:where is geometric matrix, expressed as follows:
The unit stress can be expressed simply as follows:
The unit stiffness matrix calculation formula is expressed as follows:
The unit mass matrix expression formula is expressed as follows:where is the per unit length weight of the mining pipe.
The unit stiffness and mass matrix can be obtained by calculating equations (12) and (13) as follows:
In equations (15) and (16):where is the mining pipe cross-sectional area.
3.2. Overall Matrix Analysis
Assemble the unit matrix obtained in the previous section into a global matrix; the 5000 m mining pipe is divided into 250 units, as shown in Figure 4; each unit is 20 m in length, with a total of 251 nodes and 753 degrees of freedom, and the overall stiffness and mass matrix are the square matrix of .
The nth element stiffness matrix is , n = 1, 2, …, 250. Its elements are , where r and s range from 1 to 6. The element in the nth element stiffness matrix is added to the element in the (n + 1)th element stiffness matrix; this is used as the anchor point to add the corresponding elements, and the overall stiffness matrix is obtained as follows:
In the same way, the element in the nth element mass matrix is added to the element in the (n + 1)th element mass matrix; this is used as the anchor point to add the corresponding elements, and the overall mass matrix is obtained as follows:
3.3. External Load Analysis
The mining pipe work in ocean environment and the forces are very complicated. The external load in the whole work mainly includes the marine fluid power provided by the combined action of wave and current, the gravity of the mining pipe subsystem itself and the buoyancy, etc. The middle ore bin is considered as the concentrated mass force. Here, only the longitudinal vibration of the mining pipe is considered under the combined action of wave and current after steady state balance the mining pipe. As shown in Figure 5, in the actual sea conditions, the mining pipe will have a slight deviation under the action of the current, and the vertical vibration will occur under the action of waves, which will convert the distributed marine hydrodynamic load into equal force added on the unit node [18]. The calculation formula of tangential force of the unit length mining pipe at depth x [19] (i = 1, x = 0; i = 2, x = 20; …; i = 251, x = 5000):
The normal force calculation formula is as follows:
In equations (53)–(60), is the tangential resistance coefficient; generally, the tangential resistance coefficient is 30 to 120 times smaller than the normal resistance coefficient [20]; is the normal resistance coefficient; is the horizontal velocity of fluid under combined action of wave and current; is the horizontal velocity of fluid under wave action; is the horizontal velocity of fluid under the action of ocean currents at different depths; x is the depth of water, the ocean surface is 0, and the downward value is negative; is the height of the wave; is the wave number; is the wave length; is the wave motion frequency; is the period of wave motion.
The external force matrix is expressed as follows:
The load force of the node is expressed as follows:
3.4. Equation Solving
Its dynamic equation is expressed as follows [21]:where C is the damping matrix. In general, the damping of the structure is difficult to be determined accurately; in practice, the approximate value of the damping is often given by the measured data, so the overall damping matrix of the structure can be directly given; here is the form of the product of the viscous damping coefficient c and the matrix [22].
The top of the mining pipe is subject to the harmonic heave motion of the mining ship:
The direct integration method is used to obtain the recursive formula of displacement and velocity of nodes as follows [23]:
At time t to , the iterative calculation is needed. The new state quantity is iterated according to equations (65)–(67). The displacement and velocity at the next moment can be obtained, and finally the displacement and velocity in the whole time process can be obtained.
4. Results and Discussion
The 5000 m mining pipe is divided into four step sections with lengths of 1000 m, 1000 m, 1500 m, and 1500 m, respectively. The specific parameters are shown in Tables 1 and 2. The vibration displacement under different weights, offset angles, and the vibration frequency characteristics under different conditions are analyzed, and the corresponding results are obtained.
4.1. Natural Frequency Calculation
The upper end of the mining pipe is connected with the mining ship and moves with the ship. The lower end is connected with the middle ore bin, which is connected with the mining machine through the flexible pipe, and the damping device is placed at the end of the first step pipe to provide variable damping c. In addition, the pipe is also subjected to the external viscous damping of seawater, which is smaller than the controllable damping, and the value is set as [24]. The external dimension of the middle ore bin is relatively smaller than that of the 5000 m mining pipe, so only the influence of its weight m on the mining pipe is considered. The stiffness of the flexible pipe is small and the buoyant material is attached, and the force on the mining pipe is negligible [25]. The lower end of the mining pipe is treated as the free end. The first three rows and columns of corresponding elements in matrix K and M were removed, and the stiffness matrix and mass matrix were obtained. The natural frequency calculation formula is as follows [26]:
After calculation, the first 20 order natural frequencies of the mining pipe are shown in Table 3.
4.2. Longitudinal Vibration Analysis of the Mining Pipe under Different Ore Bin Weights
Moore’s empirical formula is used to calculate the heave motion amplitude of the mining ship under the working conditions of level 6, and its value is 1.208 m [27], calculations in MATLAB.(i)The mass of the ore bin is 0 t, and the damping is . When the mining pipe is vertical to sea level, its longitudinal displacement vibration response is shown in Figure 6. In order to ensure that the longitudinal vibration displacement change effect is obvious, the change rule diagram of damping is selected here, as shown in Figure 6. The longitudinal displacement vibration response of the mining pipe is irregular motion, which presents an oscillation phenomenon. The maximum longitudinal vibration displacement of the mining pipe at 1000 m–5000 m increases gradually from top to bottom, and the maximum vibration amplitude at 1000 m is 0.1650 m and at 5000 m is 0.3103 m. The response time of receiving the vibration signal increases gradually from top to bottom.(ii)When the mass of the ore bin is 0 t, 30 t, 40 t, and 50 t and the damping is , , and , the maximum vibration amplitude at different positions of the mining pipe is shown in Figure 7.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
It can be seen from Figures 7(a)–7(c) that the overall longitudinal vibration amplitude of the mining pipe decreases first and then increases from top to bottom. The maximum value of displacement vibration occurs at the top of the mining pipe, which is stimulated by the mining ship. With the increase of the pipe depth, the vibration amplitude gradually decreases, the vibration amplitude decreases to the minimum at 1000 m and then gradually increases with the increase of the pipe depth. This is because at 0-1000 m, the contact area between the mining pipe nodes is large, the force is strong, the energy consumed by the acting force is large, so the amplitude tends to decrease. At 1000 m–5000 m of the mining pipe, the diameter of the pipe is reduced, the action of inertia force is enhanced, and the forces between nodes are weakened, so the amplitude of the pipe shows the upward trend. As the weight of the ore bin increases, the inertia force effect increases, and the overall vibration amplitude of the mining pipe will also increase. In addition, increasing damping can significantly reduce the overall vibration amplitude of the mining pipe.
4.3. Longitudinal Vibration Analysis of the Mining Pipe under Different Offset Angles
Under the working wind condition of level 6, when the offset angles are 0°, 0.15°, and 0.3°, the weight of the ore bin is 0 t, 30 t, and 40 t, and the damping is . The maximum vibration amplitude at different positions of the mining pipe is shown in Figure 8.
(a)
(b)
(c)
According to Figures 8(a)–8(c), when the offset angle of the mining pipe is increased, its overall movement trend remains unchanged, but increasing the offset angle can increase the overall vibration amplitude of the mining pipe. This is because when the mining pipe has the slight angle offset, the longitudinal exciting force transferred by the mining ship to the mining pipe can be considered to be all added to the axial direction of the mining pipe, while the components in other directions are very small and negligible. However, with the appearance and increase of the offset angle, the component force of marine fluid power along the pipe’s axial direction increases, and the overall resultant force increases, so the overall vibration amplitude of the pipe will rise. Similarly, when there is a slight offset angle, the weight of the ore bin increases, and the overall vibration amplitude of the mining pipe increases.
4.4. Vibration Frequency Analysis under the Working Wind Condition of Level 6
(i)Under the working wind condition of level 6, the offset angle of the mining pipe is 0°, the weight of the ore bin is 0 t, 30 t, and 40 t respectively, and the damping is . The longitudinal vibration frequency of the mining pipe is shown in Figure 9. As shown in Figure 9, the longitudinal vibration intensity of the mining pipe is enhanced with the increase the ore bin weight, but the increase in the ore bin weight cannot change the vibration frequency of the pipe. The maximum vibration intensity occurs at the first-order vibration frequency (28.125 Hz). When the weight of the ore bin is 0 t, 30 t, and 40 t, the corresponding vibration intensity is 10.121 dB, 10.624 dB, and 11.109 dB.(ii)Under the working wind condition of level 6, the weight of the ore bin is 30 t, the offset angle of the mining pipe is 0°, 0.15°, and 0.3°, and the damping is . The longitudinal vibration frequency of the mining pipe is shown in Figure 10. As shown in Figure 10, the vibration intensity of the pipe increases with the increase of the offset angle. Similarly, the vibration frequency of the pipe cannot be changed with the change of offset angle. The maximum vibration intensity occurs at the first-order vibration frequency (28.125 Hz). When the offset angle of the pipe is 0°, 0.15°, and 0.3°, the corresponding vibration intensity is 10.121 dB, 10.624 dB, and 11.109 dB.(iii)The first three order vibration frequencies are selected as the research object, and the values are 28.125 Hz, 31.152 Hz, and 32.031 Hz, distance from the natural frequency. The first-order vibration frequency is the closest to the natural frequency, the third-order vibration frequency is the furthest from the natural frequency, and the second-order vibration frequency is in the middle position. Under the working wind condition of level 6, the weight of the ore bin is 30 t, the offset angle of the mining pipe is 0°, and the damping is , , and . The vibration intensity attenuation diagram of the mining pipe is shown in Figure 11.
(a)
(b)
(c)
As shown in Figures 11(a)–11(c), damping can reduce the intensity of the vibration, but cannot change the frequency of the vibration response. When the damping increases from to , the vibration intensity decreases greatly, and damping increases from to , and the vibration intensity decreases slowly. When the damping is , the vibration intensity of the first three order vibration frequencies are 10.624 dB, 9.896 dB, and 9.604 dB. When the damping of is added, the vibration intensity was reduced by 5.421 dB, 4.963 dB, and 4.678 dB, respectively. It shows that the closer the vibration response frequency is to the natural frequency, the greater the vibration intensity is. When the damping is added, the greater the attenuation degree of the vibration intensity is and the more obvious the damping effect is.
4.5. Vibration Frequency Analysis under Different Working Wind Conditions
(i)The weight of the ore bin is 30 t, the offset angle of the mining pipe is 0°, and the damping is , under the working wind condition of levels 4 and 6. The longitudinal vibration intensity response diagram of the pipe is shown in Figure 12. It can be seen from Figure 12 that the vibration intensity is different under different working conditions, and the frequency of vibration is also different. In general, the vibration intensity of the level 6 working wind condition is greater than that of the level 4 working wind condition. The vibration frequency of each order under the working wind condition of level 6 is less than the corresponding vibration frequency of each order under the working wind condition of level 4. It shows that the closer the frequency of wave excitation is to the natural frequency of the mining pipe, the larger the vibration is.(ii)The first three order vibration frequencies are also selected as the research object, and the values are 28.223 Hz, 31.25 Hz, and 32.129 Hz, distance from the natural frequency. The first-order vibration frequency is the closest to the natural frequency, the third-order vibration frequency is the furthest from the natural frequency, and the second-order vibration frequency is in the middle position. Under the working wind condition of level 4, its heave motion amplitude is 1.018 m, the weight of the ore bin is 30 t, the offset angle of the mining pipe is 0°, and the damping is , , and , respectively; the vibration intensity attenuation diagram of the mining pipe is shown in Figure 13.
(a)
(b)
(c)
As shown in Figures 13(a)–13(c), when the damping is , the vibration intensity of the first three order vibration frequencies are 8.514 dB, 7.942 dB, and 7.752 dB. Similarly, it can be obtained that the closer the vibration frequency is to the natural frequency, the greater the vibration intensity will be, and the vibration intensity will decay faster when damping is added. And the degree of vibration attenuation from to is greater than the damping from to . It can be seen from Figures 11(a) and 13(a) that when 400 damping is added, the vibration attenuation value is 5.421 dB under the working wind condition of level 6 and 4.142 dB under the working wind condition of level 4. It shows that the closer the vibration frequency is to the natural frequency of the mining pipe, the greater the vibration is generated, and the greater the vibration attenuation rate is after adding damping.
5. Conclusions
In this paper, the longitudinal vibration of the 5000 m stepped mining pipe is studied. Based on the finite element method, the pipe is simplified into beam element and discretized, and the theoretical model of the pipe longitudinal vibration is established. The direct integration method is used for numerical calculation. The research content has important guiding function for the design of deep sea mining system, the research of heave compensation system, and the actual deep sea mining operations. The main conclusions are as follows:(1)The longitudinal displacement vibration response of the mining pipe is irregular movement, which presents the phenomenon of oscillation. The vibration amplitude decreases first and then increases from top to bottom, and the minimum vibration amplitude appears at 1000 m, whose value is 0.1650 m. The overall maximum vibration amplitude appears at the top of the mining pipe.(2)Under the same working wind condition, the mining pipe offset angle and damping are unchanged. Due to the influence of inertia force, the larger the weight of the ore bin is, the larger the longitudinal vibration amplitude is.(3)Under the same working wind condition, the ore bin weight and damping are unchanged. When the mining pipe is slightly offset in seawater, the larger the angle of slight offset is, the larger the longitudinal vibration amplitude is.(4)The offset angle, the ore bin weight, and the applied damping cannot change the frequency of longitudinal vibration. The closer the vibration frequency generated by the same excitation frequency is to the natural frequency, the greater the vibration intensity is. When damping is added, the vibration intensity decreases faster.(5)The closer the excitation frequency generated by different working wind conditions is to the natural frequency, the larger the mining pipe longitudinal vibration amplitude is, the greater the influence of external damping is, and the more obvious the damping effect is.
Appendix
A. The Properties of Materials and Geometry
function [k, m] = FrameElement(prop) E = prop(1); u = prop(2); rho = prop(3); D = prop(4); d = prop(5); A = prop(6); G = prop(7);
B. Element Stiffness and Mass Matrix
(1)function y = spaceFrameElementstiffness(E, A, x1, y1, z1, x2, y2, z2) L = sqrt((x2 − x1) ∗ (x2 − x1) + (y2 − y1) ∗ (y2 − y1) +(z2 − z1) ∗ (z2 − z1)); kprime = (E ∗ A/L) ∗ [17/5 1 −1 (L − 11)/10 −1 1; 1 17/5 −6/5 0 (L − 11)/10 −L/10; −1 −6/5 17/5 0 L/10 (L − 11)/10; (L − 11)/10 0 0 (14 ∗ L^2 − 12)/15 0 0; −1 (L − 11)/10 L/10 0 (14 ∗ L^2 − 12)/15 −L/30; 1 −L/10 (L − 11)/10 0 −L/30 (14 ∗ L^2 − 12)/15]; y = kprime;(2)function y = spaceFrameElementmass (rho, A, x1, y1, z1, x2, y2, z2) L = sqrt((x2 − x1) ∗ (x2 − x1) + (y2 − y1) ∗ (y2 − y1)+ (z2 − z1) ∗ (z2 − z1)); = rho ∗ A ∗ L mass = ∗ L ∗ [113/105 7/20 3/20 11 ∗ L/210 − 4/35 3/20 7/20; 7/20 113/105 9/70 L/20 11 ∗ L/210 − 4/35 13 ∗ L/420; 3/20 9/70 113/105 −1/30 −13/420 11 ∗ L/210 − 4/35; 11 ∗ L/210 − 4/35 L/20 −1/30 L^2/105 + 12/35 L/30 −1/20; 3/20 11 ∗ L/210 − 4/35 −13/420 L/30 L^2/105 + 12/35 –L/140; 7/20 13 ∗ L/420 11 ∗ L/210 − 4/35 −1/20 –L/140 L^2/105 + 12/35]; y = mass;
C. Assembly Matrix
(1)function y = SpaceFrameAssemble(K, k, i, j) DOF(1) = 3 ∗ i − 2; DOF(2) = 3 ∗ i − 1; DOF(3) = 3 ∗ i; DOF(4) = 3 ∗ j − 2; DOF(5) = 3 ∗ j − 1; DOF(6) = 3 ∗ j; for n1 = 1 : 6 for n2 = 1 : 6 KK(DOF(n1), DOF(n2)) = KK(DOF(n1), DOF(n2)) + k(n1, n2); end end y = KK(2)clc clear k1 = spaceFrameElementstiffness(E, A, 0, 0, 0, L, 0, 0); k2 = spaceFrameElementstiffness1(E, A1, 0, 0, 0, L, 0, 0); k3 = spaceFrameElementstiffness2(E, A2, 0, 0, 0, L, 0, 0); k4 = spaceFrameElementstiffness3(E, A3, 0, 0, 0, L, 0, 0); K = zeros(n1, n2); K = SpaceFrameAssemble(K, k1, 1, 2); K = SpaceFrameAssemble(K, k2, 51, 52); K = SpaceFrameAssemble(K, k3, 101, 102); K = SpaceFrameAssemble(K, k4, 176, 177); K = SpaceFrameAssemble(K, k4, 250, 251);
D. Calculation Program
function [a a1 a2] = multi_wiltheta(M, K, C, d, e, f, P, theta, t) n = length(P(:, 1)); m = length(P(1, :)); u = zeros(n, m); u1 = zeros(n, m); u2 = zeros(n, m); uu = zeros(n, m); PP = zeros(n, m); u(:, 1) = d; u1(:, 1) = e; u2(:, 1) = f; l = 6/(theta ∗ t)^2; q = 3/(theta ∗ t); r = 6/(theta ∗ t); s = theta ∗ t/2; KK=K + l ∗ M + q ∗ C; for i = 1 : m − 1 pp(:, i) = P(:, i) + theta ∗ (P(:, i + 1) − P(:, i)) + M ∗ (l ∗u(:, i) + r ∗ u1(:, i) + 2 ∗ u2(:, i)) + C ∗ (q ∗ u(:, i) + 2 ∗u1(:, i) + s ∗ u2(:, i)); uu(:, i) = KK\pp(:, i); a2(:, i + 1) = u2(:, i) + l/theta ∗ (uu(:, i) − u(:, i)) − l ∗ t ∗u1(:, i) − 3/theta ∗ u2(:, i); a1(:, i + 1) = u1(:, i) + t/2 ∗ (u2(:, i + 1)+u2(:, i)); a(:, i + 1) = u(:, i) + t ∗ u1(:, i) + t^2/6 ∗ (u2(:, i + 1)+2 ∗u2(:, i)); end end
E. FFT Program
a = load(‘dianshu.m’); fs = ft;number = N; y = fft(a, number); n = 0 : length(y) − 1; f = fs ∗ n/length(y); subplot(211) plot(f, abs(y)); title(‘a’) xlable(‘Hz’) Re = real(y); Im = imag(y); grid on
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (allotment grant no. 51774193) and the Shandong Provincial Natural Science Foundation, China (allotment grant no. ZR2017MEE025).