Dynamic Analysis of Overhead Transmission Lines under Turbulent Wind Loading ()

Alok Dua^{1}, Mathias Clobes^{2}, Thomas Höbbel^{2}, Vasant Matsagar^{1}

^{1}Department of Civil Engineering, Indian Institute of Technology (IIT), Delhi, India.

^{2}Institut für Stahlbau, Technische Universität Carolo-Wilhelmina, Braunschweig, Germany.

**DOI: **10.4236/ojce.2015.54036
PDF
HTML XML
5,157
Downloads
7,460
Views
Citations

Transmission tower-line systems are designed using static loads specified in various codes. This paper compares the dynamic response of a test transmission line with the response due to static loads given by Eurocode. Finite element design software SAP2000 was used to model the towers and lines. Non-linear dynamic analysis including the large displacement effects was carried out. Macroscopic aspects of wind coherence along element length and integration time step were investigated. An approach is presented to compare the probabilistic dynamic response due to 7 different stochastically simulated wind fields with the response according to EN-50341. The developed model will be used to study the response recorded on a test line due to the actual wind speed time history recorded. It was found that static load from EN overestimated the strength of conductor cables. The response of coupled system considering towers and cables was found to be different from response of only cables with fixed supports.

Share and Cite:

Dua, A. , Clobes, M. , Höbbel, T. and Matsagar, V. (2015) Dynamic Analysis of Overhead Transmission Lines under Turbulent Wind Loading. *Open Journal of Civil Engineering*, **5**, 359-371. doi: 10.4236/ojce.2015.54036.

1. Introduction

Collapse of transmission tower-line systems is not a well understood phenomenon. These systems are subjected to various loads like wind, snow, icing and earthquake. Comparatively, wind loads are more complex for these structures due to high geometric non-linearity of cables and randomness of wind turbulence. This thesis is aimed at understanding the dynamic behavior of transmission tower-line systems under fluctuating wind loads. Present code recommendations are based on static loading. In this paper, [1] was considered for comparison of results. EN-50341 gives wind pressure including 2 second gusts with peak wind velocities, for conductors, insulators and towers. In the present design practice, towers and conductors are considered separately ignoring the coupling effect and static loads are applied individually.

Coupled transmission tower-line systems are highly complex in their behavior due to the interaction between non-linear conductors and stiff towers which results in closely spaced frequencies. [2] studied these systems considering the geometric non-linearity and aerodynamic damping; however the work was not in 3 dimensions. [3] presented linear response in a coupled system using a 3D finite element model. [4] recently showed that the design codes overestimate the strength of transmission towers. They brought out that a 3D finite element analysis is more accurate compared to linear analysis. Until now, not many researchers have studied the coupling effect on response of cables using non-linear dynamic analysis. Previous studies were linear and without the effects of large displacements. Secondly, the response of such systems is usually assumed as Gaussian for convenience in calculating the extreme response.

A 3D non-linear analysis including the large displacements was carried out to study the dynamic response of conductors. Effects of parameters like coherence along element length and integration time step were considered. The response of cables and insulators was found to be non-Gaussian. Methods lately published for calculating extreme values of a non-Gaussian process were used. These extreme values were then compared to the response due to wind pressure recommended in EN-50341. The chosen transmission line is in Rostock, Germany. It has 2 end towers and 2 suspension towers. Each of the 3 spans is about 400 m. A 3D finite element model of a real transmission line was modeled in finite element software SAP2000. Two models were studied: only conductors and conductors coupled with towers. The results were compared to show the importance of conductor-tower interaction.

2. Details of the Model

2.1. Chosen Test Line

The test line is a 380 kV line with 2 circuits and 3 phases. The conductors are 2 × 3 quad bundle made of aluminum conductor with steel core (ACSR). There are 2 end towers referred as WA15 and WA18 and 2 supporting towers referred as T17 and T16 (Figure 1). As part of cooperative research project between TU Braunschweig and BAM (National Institute for Materials Research and Testing, Berlin) the line is being monitored with anemometers, accelerometers, tension load cells and tilt sensors. The anemometers, accelerometers and load cells are installed on the second traverse of the towers and the tilt sensors are installed on T17. In this paper these measurements have not been discussed.

2.2. Modelling of Towers

SAP2000v15 OAPI (Open Application Programming Interface) was used to create the geometry of towers. A typical transmission tower can have close to 1500 members and 500 nodes. To recreate the towers with varied slenderness ratios can be time consuming. The developed VB code solved this problem. Parameters like slope of main columns, base and top width, height of the broader part of tower and total height were taken as input.

Figure 1. Layout of the test line.

In past, [5] - [10] have brought out the importance of connections in transmission towers. Only geometric and material non-linear analysis can fully depict the behavior. In this study, the towers were modeled as a beam-truss model. Rigid connections (with two or more bolts) were modeled using beam elements. Flexible connections (single bolt) were modeled using truss elements with moment released about appropriate direction. This is an approximation and ideally flexible connections should have some stiffness. The eccentricity in the connections and the load application point has been ignored. European norms recommend these structures to be in elastic range during service life. Hence, material non-linearity has not been considered in the model.

The created tower geometry was checked for disjointed nodes. The towers have been created in SAP2000 using open application programming interface (OAPI) by taking required height as input. Each member length is numerically calculated and geometry is created in SAP2000 v15. This sometimes results in common nodes falling out of a member while it was created with other members. These occurrences are very less and can be pin pointed by dead load analysis. After rectifying the disjointed nodes the sectional and material properties were defined for each member. Each member was discretized into 3 parts to ensure adequate accuracy. A dead load analysis was performed on these tower models. Characteristics of the 3 tower types are given in Table 1. The stiffness matrix from dead load analysis was used as initial condition for the modal analysis of towers. This was done to ensure that the member forces due to dead load are accounted for in the modal analysis. The first mode shape and frequencies are shown in Figure 2 for the 3 towers.

2.3. Modelling of Conductors

Each tower is at a different ground level and this has been accounted for. Distances and sag in each span is given in Table 2.

Figure 2. First mode shape and frequencies of the towers.

Table 1. Characteristics of towers.

Table 2. Distances and sag in span.

The conductors were modeled as tension-only linear elastic material. The non-linearity of these flexible cables was taken into account for dynamic analysis. Quad bundle of conductors was assumed as one conductor with 4 times cross sectional area. In this study, each conductor was divided into 20 parts for application of wind load time histories. The sectional and material properties for conductors are given in Table 3.

2.4. Damping in Conductors

The sectional and material properties for towers, insulator strings and conductors are as per the technical drawings. Realistic modeling of aerodynamic damping is very complex for such systems. Aerodynamic damping affects the conductors in a varied way. Aerodynamic damping is an aeroelastic phenomenon and opposes the action of cables depending on the direction of motion. This cannot be modeled in SAP2000v15 however a satisfactory model has been recently presented by [11] using commercial software ADINA. For this study the aerodynamic damping was incorporated as viscous damping. An equivalent viscous material damping ratio of 2% is sufficient to model the aerodynamic damping effects on the conductors [3] ; [8] . Apart from aerodynamic damping in conductors, 0.5% damping ratio was used to account for structural damping due to steel material, connections and the foundation.

2.5. Insulator Strings

In most of the numerical research works, either the insulator strings have not been modeled or have been assumed to be a beam element while towers have been neglected. [12] have discussed about the importance of insulator string in overhead transmission line under wind load. As per their study, the non-rigid insulator strings have to be modeled as per the real material, sectional and boundary properties. It is not accurate to assume the insulator strings as rigid beam elements. In this study, the insulator strings were modeled as cable elements to account for the local slackening effects in the flexible insulator strings. However, a better model as discussed by [12] may be used to get more accurate results for failure criteria of insulators.

2.6. Simplification of Towers

The processing time depends on the computing platform however with best commercially available platform also the software and degrees of freedom can be a restriction. The analysis time for complete SAP2000 model was about 60 - 75 hours due to large number of degrees of freedom. To reduce the analysis time, the towers were reduced to equivalent beams. As presented by [13] , latticed towers can be reduced to beams with equivalent stiffness and material. The towers were first divided into segments with same sectional properties. The equivalent axial, flexural, torsional, shear stiffness and mass were determined for each segment. Displacement of the complete towers and equivalent towers were compared and a good match was observed. The modal stiffness of the equivalent towers was also close to the values from the corresponding towers (Table 4). The comparison of displacements in first mode shape of lattice tower T17 and its equivalent tower is shown in Figure 3.

2.7. Reduced Model

Two models were made to compare the effect of interaction between towers and conductors. The coupled model with reduced tower is shown in Figure 4. Second model had only the cables from two spans of the test line. An insulator was modeled as the center support. The two end supports were considered to be fixed and the support for the insulator was modeled as a pinned support.

Table 3. Sectional and material properties for conductors.

Figure 3. Comparison of displacements in first mode of T17 and its equivalent tower.

Figure 4. Reduced model.

Table 4. Comparison of modal stiffness’s.

3. Wind Forces

Separate time histories of 300 seconds were generated for 19 points on each conductor. Distance between each point on the conductor is about 20 m. The mean wind speed at respective heights for generating the time history is from DIN-1055-4:2005-03 [14] . Wind turbulence is modeled using the weighted amplitude wave superposition (WAWS) model based on Shinozuka and Jan [15] . Details of simulation can be found in the work done by Clobes [16] . Von *Kármán* power spectral density function was used to characterize the power distribution of the turbulence in longitudinal direction, Kaimal for lateral direction and Busch and Panofsky spectrum for vertical turbulence. The cross-correlation of two neighboring points decreases with increasing distance between them. This point has to be kept in mind as the loads generated can be up to 10 - 20 times higher if each element length is very large. At low frequencies, the eddies have a large integral length scale and take long time to cross the structure. In this case the distribution of load equally over the length of element is justified. However, for eddies corresponding to higher frequencies, that are smaller than the element length, it is incorrect to consider the wind load fully coherent along the length of the element. The high frequency eddies actually compensate each other from one point on element to other. If the load along the element is considered same it could result in overestimation of the forces [17] . While conducting this work, initially a model was made with element length of 20 m and it was found that the tension in the conductor was almost 65% - 70% higher than the tension due to design loads. Equation (1) [17] gives the ideal element length, below which the forces can be considered to be coherent.

(1)

where C is the decay factor taken as 11 (Table 4.2, pg 51, [16] ), n is the frequency of eddy and is the mean wind speed at that height. It was calculated from equation 1, that wind loads generated at a distance of 20 m are not coherent and need to be rectified. This effect of coherence was corrected using the joint acceptance function given in Clobes [16] . One of the simulated longitudinal wind speed time history at a point on the cable is shown below (Figure 5).

Mean wind speed was then added to the turbulence generated for each point. The drag force coefficient (C_{d}) for conductors was taken to be 1 (EN-50341). A force time history was generated for each point using , where ρ is air density, U is the longitudinal wind speed, d is the diameter and L is the length of exposure. The solution could not converge for the first time step as the load at the first time step was suddenly applied. To avoid this, the time history was ramped for first 10 seconds.

4. Dynamic Analysis

There are 16 conductors in two spans and each divided into 20 parts, hence 304 time histories were simulated in one wind field. Seven such fields were simulated for probabilistic analysis. Each time history was applied to the model which was a time consuming task. To reduce the efforts of generating the loads and applying the time history, OAPI was used. Hilber-Hughes-Taylor time integration method (1977) was used for direct integration initially with α = −0.33. Integration time step was taken as 0.1 s at first to get results. To study the effect of integration time step a comparative study was done. Figure 6 shows the effect of integration time step on tension in wind ward cable in the longer span. The extreme values from a time step 0.005 and 0.01 were very close. Al-

Figure 5. Time history of wind speed and power spectral density plot (along-wind direction).

Figure 6. Effect of integration time step.

though 0.01 was the right time step, due to constraints of SAP2000’s memory usage, a time step of 0.02 s was used. With an integration time step of 0.01 s the whole time history could not be solved due to inadequate system memory. The system being used is a dual core processor with 64 bit OS and 24 GB RAM. We are of the opinion that the multi-threaded solver in SAP2000 could not recognize the 64 bit system so as to enable it to use the whole system memory (more than 4 GB). This error may be attributed to the .NET communication between SAP2000 and the 64 bit OS. However, a detailed investigation of this issue is still under process.

A separate static model was created on which static loads were applied. The design loads for overhead electrical lines exceeding 45 kV are taken from EN 50341. The variation of wind pressure along height is given as per Equations (2) and (3).

(2)

(3)

where h is the height above ground level in meters and is the reference wind pressure including gusts with peak wind velocities (2 sec gusts). The test line that was modeled is present in wind zone 2 for Germany and hence has been taken as 390 N/m^{2}. The wind loads on conductors, insulators and towers were calculated as per Equations (4), (5) and (6) respectively.

(4)

(5)

(6)

where q is the wind pressure as per Equations (2) and (3), G is the dynamic response factor given by 0.45 + 60/L for spans greater than 200 m and as 0.75 for spans lesser than 200 m, d is the diameter of the conductor, L is the length of conductor exposed, is the angle between wind direction and the cross arms, is the exposed area of insulator, is the area on the tower exposed to wind, is the drag force coefficient taken as 2.8 for lattice towers (EN-50341). A static non-linear analysis was done for this model with above design loads.

5. Results and Discussions

5.1. Stochastic Analysis

In time history analysis due to wind loads, extreme value of the response is generally found out. This is due to the fact that the applied wind loads are generated from a random process and the extreme value of the response will vary with each time history. In probabilistic analysis this is called as mean extreme value. It is the mean value of probability density function for extreme values given by Equation (7) [18] .

(7)

where ν is the mean frequency of occurrence of zero crossings with positive slopes only and is given by Equation (8).

and

(8)

here is the power spectral density function of the random process.

It has been shown by [19] that the mean extreme value can be given with an approximate relation (Equation (9)) derived from Equation (7).

(9)

here is Euler’s constant (0.5772).

Equation (9) is valid only for a Gaussian process. However, it was observed that the response of the structure is Gaussian (Figure 7) for none of the wind fields.

The skewness and kurtosis for the 4 response parameters for the 7 wind fields can be seen in Table 5 which shows that response is not Gaussian. [20] and [21] have presented approaches to calculate the peak factors for non-Gaussian stationary processes. These 2 methods were chosen to compare the extreme values for the response.

5.1.1. Hermite Moment-Based Method

Kareem and Kwon [20] expressed the Hermite moment-based non-Gaussian peak factor as given in Equation (10).

(10)

where is Euler’s constant (0.5772);; is the frequency of occurrence of zero crossings with positive slopes only (Equation (8)); parameters and control the shape of the distribution and is the scaling factor given by Equation (12) [22] .

Figure 7. Probability distribution of cable tension.

Table 5. Peak factors from considered approaches to evaluate peak factors.

(12)

5.1.2. Skewness Dependent Peak Factor

[21] studied the peak factor of mild non-Gaussian process. They recommended a simplified empirical formula for non-Gaussian peak factor dependent only on the skewness but effect of mild softening has been empirically calibrated (Equation (12)).

(13)

where the variables have the same meaning as in Equation (10).

5.2. Comparison of Results

Figure 8 shows the time histories of four selected parameters from the coupled model having both cables and the towers. The response is of 1 of the 7 wind fields.

Peak factors for these parameters using the above methods are shown in Table 5. Experimental works by [21] , [23] and [24] proves that the peak factors from [19] method not accurate as the process becomes non-Gaussian. The peak factors from [19] are always close to 3. Peak factors with [20] method are close to experimental peak factors from wind tunnel tests by [21] . However, if the kurtosis is too high or too low the peak factors from [20] are also higher and lower respectively than the experimental values. The same trend can be seen in the results in Table 5. Based on the above observation, for further comparisons the peak factor from [21] was chosen.

Extreme values based on the peak factors from [21] method are shown in Table 6. Extreme value of a particular parameter is a probabilistic value obtained by assuming the response as a non-Gaussian process. It can be seen that the extreme value changes for each wind field. It is a common practice to consider the highest extreme

Figure 8. Response time history of 4 parameters.

value from the set of extreme values obtained. The last column of the table shows these values that have been considered for comparison with the response due to static wind loads given by EN-50341.

In Table 7 these extreme values are compared with the results from static loads.

There is a difference of 12-30% between the two responses for various parameters. To investigate the effect of coupling between the towers and the cables, the responses of the two models were compared. The results are shown in Table 8.

There is a difference in the results from two models. Hence, it is more accurate to consider the coupled model for analysis of such structures. At the same time the coupled model increases the analysis time considerably even with the reduced towers. In addition the accuracy of peak factors is less as the time step used was 0.02 s.

6. Conclusions

In the interest of dynamic response of power transmission lines, the following main conclusions were made from this study:

Aerodynamic damping is an influential parameter in the response of transmission lines. For this study an equivalent viscous damping suggested on the basis of experimental results gave satisfactory results. However recently [11] have presented a new method to determine wind response of transmission lines using fluid-structure interaction. This method gives a more accurate representation of wind loads acting on moving conductors.

The effect of wind coherence along the element length in the simulated stochastic wind field is of great importance. It was found that if the element length is larger than the eddy size, the effective wind load on that element can be up to 45% larger. Two options are suggested based on this study: firstly, considering an element length which is smaller than the eddies with higher frequencies; or secondly, reducing the force on an element based on the joint acceptance function.

Table 6. Extreme values response of 4 parameters for 7 wind load time histories.

Table 7. Comparison of considered extreme values and response from static wind load.

Table 8. Comparison of response from two models.

Response of a transmission line to gust wind is non-Gaussian in nature. An appropriate method to calculate the peak factor for a non-Gaussian random process gave results that were higher than the davenport’s peak factors. A numerical verification of these peak factors can be envisaged as a future research task.

Static wind loads specified in EN-50341 overestimate the cable strength. Extreme values of 4 response parameters were found to be greater than the static design response.

The response for a coupled model is different from the response when only the cables are considered. The results from the model with only the cables gave wrong estimations for the displacements of lines and insulators. The swing angle of insulators can be larger in coupled model thereby resulting in flashovers. To account for the difference in displacement values, it is recommended to use an equivalent stiffness at supports, instead of the towers, in the model with only cables. This can reduce the analysis time to 3 - 4 hours and also satisfactorily account for the difference in the response as compared to the coupled model.

As a future scope for the project, it would be interesting to input the wind speed records from the test line and compare the numerical response with the recoded response. The insulator swing angle can be conveniently calculated from the numerical model and the same parameter is being recorded at the test line.

Conflicts of Interest

The authors declare no conflicts of interest.

[1] | European Committee for Standardisation (2010) DIN-EN-50341-3-4-VDE-0210-3, Overhead Electrical Lines Exceeding AC 45 kV. Part I—General Requirements—Common Specifications. European Committee for Standardisation, Germany. |

[2] |
Yasui, H., Marukawa, H., Momomura, Y. and Ohkuma, T. (1999) Analytical Study on Wind-Induced Vibration of Power Transmission Towers. Journal of Wind Engineering and Industrial Aerodynamics, 83, 431-441.
http://dx.doi.org/10.1016/S0167-6105(99)00091-4 |

[3] |
Battista, R.C., Rodrigues, R.S. and Pfeil, M.S. (2003) Dynamic Behavior and Stability of Transmission Line Towers under Wind Forces. Journal of Wind Engineering and Industrial Aerodynamics, 91, 1051-1067.
http://dx.doi.org/10.1016/S0167-6105(03)00052-7 |

[4] |
Rao, N.P., Légeron, F. and Prud’homme, S. (2012) Variation of Damping and Stiffness of Lattice Towers with Load Level. Journal of Constructional Steel Research, 71, 111-118. http://dx.doi.org/10.1016/j.jcsr.2011.10.018 |

[5] |
Robert, V. and Lemelin, D.R. (2002) Flexural Consideration in Steel Transmission Tower Design. In: Electrical Transmission in a New Age, ASCE, Omaha. http://dx.doi.org/10.1061/40642(253)11 |

[6] |
Albermani, F.G.A. and Kitipornchai, S. (2003) Numerical Simulation of Structural Behavior of Transmission Towers. Journal of Thin-Walled Structures, 41, 167-177. http://dx.doi.org/10.1016/S0263-8231(02)00085-X |

[7] |
da Silva, J.G.S., da Vellasco, P.C.G., de Andrade, S.A.L. and de Oliveira, M.I.R. (2005) Structural Assessment of Current Steel Design Models for Transmission and Telecommunication Towers. Journal of Constructional Steel Research, 61, 1108-1134. http://dx.doi.org/10.1016/j.jcsr.2005.02.009 |

[8] |
McClure, G. and Lapointe, M. (2003) Modeling the Structural Dynamic Response of Overhead Transmission Lines. Journal of Computers and Structures, 81, 825-834. http://dx.doi.org/10.1016/S0045-7949(02)00472-8 |

[9] |
McClure, G. and Lee, P.S. (2007) Elastoplastic Large Deformation Analysis of a Lattice Steel Tower Structure and Comparison with Full Scale Tests. Journal of Constructional Steel Research, 63, 709-717.
http://dx.doi.org/10.1016/j.jcsr.2006.06.041 |

[10] |
McClure, G., Jiang, W.Q., Wang, W.L. and Geng, J.D. (2011) Accurate Modeling of Joint Effects in Lattice Transmission Towers. Journal of Engineering Structures, 33, 1817-1827. http://dx.doi.org/10.1016/j.engstruct.2011.02.022 |

[11] |
Keyhan, H., McClure, G. and Habashi, W.G. (2013) Dynamic Analysis of an Overhead Transmission Line Subject to Gusty Wind Loading Predicted by Wind-Conductor Interaction. Computers and Structures, 122, 135-144.
http://dx.doi.org/10.1016/j.compstruc.2012.12.022 |

[12] | Yan, B., Lin, X.S., Luo, W., Chen, Z. and Liu, Z.Q. (2010) Numerical Study on Dynamic Swing of Suspension Insulator String in Overhead Transmission Line under Wind Load. IEEE Transactions on Power Delivery, 25, 248-259. |

[13] | Limongelli, M.P., Martinelli, L. and Perotti, F. (2003) A Reduced Model for the Dynamic Analysis of Power Transmission Lines with Truss Supporting Towers. Proceedings of the 5th International Symposium on Cable Dynamics, Santa Margherita, 15-18 September 2003, 125-132. |

[14] | DIN-1055-4:2005-03 (2005) Einwirkungen auf Tragwerke—Teil 4: Windlasten. Beuth Verlag GmbH, Berlin. |

[15] | Shinozuka, M. and Jan, C.M. (1972) Digital Simulation of Random Processes and Its Application. Journal of Sound and Vibration, 25, 111-128. |

[16] | Clobes, M. (2008) Identifikation und Simulation instationärer übertragung der Windturbulenz im Zeitbereich, in Fakultät für Architektur, Bauingenieurwesen und Umweltwissenschaften. Technischen Universität Carolo-Wilhelmina zu Braunschweig, Braunschweig. |

[17] | Denoël, V. (2005) Accounting for Coherence in Wind Forces in Finite Element Models. Département de Mécanique des matériaux et Structures, Université de Liège, Belgique. |

[18] | Cartwright, D.E. and Longuet-Higgins, M.S. (1956) The Statistical Distribution of the Maxima of a Random Function. Proceedings of the Royal Society of London, Series A, 237, 212-232. |

[19] |
Davenport, A.G. (1964) Note on the Distribution of the Largest Value of a Random Function with Application to Gust Loading. ICE Proceedings, 28, 187-196. http://dx.doi.org/10.1680/iicep.1964.10112 |

[20] | Kareem, A. and Kwon, D. (2009) Peak Factor for Non-Gaussian Processes Revisited. Proceedings of the 7th Asia-Pacific Conference on Wind Engineering, Taipei, 8-12 November 2009, 719-722. |

[21] |
Huang, M.F., Lou, W.J., Chan, C.M. and Bao, S. (2012) Peak Factors of Non-Gaussian Wind Forces on a Complex-Shaped Tall Building. The Structural Design of Tall and Special Buildings, 22, 1105-1118.
http://dx.doi.org/10.1002/tal.763 |

[22] | Winterstein, S.R. (1988) Nonlinear Vibration Models for Extremes and Fatigue. Journal of Engineering Mechanics, 114, 1772-1790. |

[23] | Huang, M.F., Chan, C.M., Kwok, K.C.S. and Lou, W.J. (2009) A Peak Factor for Predicting Non-Gaussian Peak Resultant Response of Wind-Excited Tall Buildings. Proceedings of the 7th Asia-Pacific Conference on Wind Engineering, Taipei, 8-12 November 2009, 423-426. |

[24] | Huang, M.F., Chan, C.M., Lou, W.J. and Kwok, K.C.S. (2012) Statistical Extremes and Peak Factors in Wind Induced Vibration of Tall Buildings. Journal of Zhejiang University, Science A (Applied Physics and Engineering), 13, 18-32. |

Journals Menu

Contact us

+1 323-425-8868 | |

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2024 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.