力学学报

ADVANCES IN DYNAMICS AND VIBRATION CONTROL OF LARGE-SCALE FLEXIBLE SPACECRAFT

Cao Dengqing, Bai Kunchao, Ding Hu, Zhou Xubin, Pan Zhongwen, Chen Liqun, Zhan Shige

With the gradual implementation of the major projects in aerospace engineering, the spacecraft design is heading towards the direction of ultra-high speed, large scale and multi-function, and its launching and operating environment is likely to worsen. The problems on vibration and active and/or passive control in spacecraft launching process, dynamic modeling and response analysis of large flexible spacecraft in orbit, and the hybrid control of structural vibration and attitude maneuver of the spacecraft are getting more and more complicated and difficult to deal with. The enlarged scale and increased flexibility of the spacecraft structure (suchas the large aperture antenna and solar panels) present a challenge to the ground test and semi-physical simulation. The dynamics and control problems involved in the large-scale flexible spacecraft such as the whole-spacecraft vibration isolation in the spacecraft lunching process, the dynamic modeling and vibration response analysis of large-scale flexible structures, and the coupling dynamics and hybrid control of structural vibration and attitude and or orbital maneuver of the large-scale flexible spacecraft are presented. The key scientific issues seriously in the fields of spacecraft dynamics and control could then be extracted as follows: the dynamic modeling and order reduction of multi-rigid flexible body systems (including the dynamic modeling of the flexible structure with large deformation, the collaborate simulations with multi-solvers, model reductions, the analytical approach for the dynamic modeling of composite structures, etc.), the construction of state space model of complicated structures and its controllability investigation (including the theoretical and experiment methods of the state space model formulation, the observability and controllability of the control system for complex structures, etc.), and the design of hybrid control law of structural vibration and attitude maneuver for the large-scale flexible spacecraft (concerning the robust hybrid control of attitude maneuver and structural vibration, the collaborative control of actuating mechanism and piezoelectric actuator, etc.

2019, 51(1): 1-13. doi: 10.6052/0459-1879-18-054

HTML PDF

2019, 51(1): 14-15.

HTML PDF

EVOLUTION PATH FOR THE PARTICLE BREAKAGE OF GRANULAR MATERIALS: A MICROMECHANICAL AND THERMODYNAMIC INSIGHT

Shen Chaomin, Liu Sihong

Particle breakage of granular materials is ubiquitous in nature and engineering practices and often takes place under high stress levels. The phenomenon of particle breakage may not only influence the mechanical response of granular materials, resulting the contraction of the volume of the material and reduction of the shearing strength, but is also closely associated to a variety of engineering problems. The existing research is mainly focused on depicting the evolution of the particle breakage and uses a quantifiable parameter to relate the particle breakage to the subsequent mechanical response. However, little attention has been paid to exploring the underlying physics of the driving force that initiates and attenuates the particle breakage. In this study, we present the formulation of an elastic-breakage model for the isotropic compression of frictionless spheres in the framework of thermodynamics. In the model, both the elastic strain energy and the dissipation due to particle breakage are formulated using the micro-macro averaging procedure, which is often used in micromechanics of granular materials. The evolution path of the particle breakage is determined using the maximum energy dissipation hypothesis. As the modelling does not involve any empirical results, all the model parameters have concrete physical meanings. Comparison of the model prediction with the experimental data in the literature showed that the initial gradation has different effects on the elastic bulk modulus and the breakage stress: the bulk modulus increase initially and then decrease with the fractal dimension of the gradation, which implies that there is a peak bulk modulus for a certain value of the fractal dimension; while the breakage stress increases monotonically with the increase of the fractal dimension. In addition, both the bulk modulus and the breakage stress increase monotonically with the increase of the polydispersity of the particle sizes. The evolution path of the gradation due to particle breakage is found to indeed satisfy the maximum dissipation hypothesis. Both experimental results and model prediction show that the compression curve of granular materials can be divided into three stages: the elastic compression stage under low compressive stress, particle breakage stage and the pseudoelastic compression after sufficiently large amount of particle breakage.

2019, 51(1): 16-25. doi: 10.6052/0459-1879-18-340

HTML PDF

CONTACT FABRIC CHARACTERISTICS OF GRANULAR MATERIALS UNDER THREE DIMENSIONAL STRESS PATHS

Liu Jiaying, Zhou Wei, Ma Gang, Li Yiao, Liu Qiwen

Macroscopic mechanical characteristics of granular materials are closely related to the microscopic contact force and fabric. Generally speaking, the strong contact system contributes to the force transmission of internal granular system, and then its corresponding fabric tensor has important influence on the macroscopic stress. Microscopic numerical methods, such as discrete element method, can reproduce the laboratory tests with reasonable macroscopic responses and extract macro- and micro-data conveniently for investigating the underlying mechanism of the granular system. Based on discrete element method (DEM), a series of true triaxial tests for granular materials under constant $p$ and constant $b$ stress paths are carried out, and the evolutions of macro- and micro-mechanical parameters of granular materials, the multiple relationship between three-dimensional fabric tensor and stress tensor and the macro-stress characteristics reflected by strong contact system are studied. The results demonstrate that some macro- and microscopic parameters at the stress peak and critical state in the granular system are independent on the loading path. Non-coaxiality between fabric tensor and stress tensor is observed under three-dimensional stress path, but the evolution of the joint invariant of the two tensors is independent on the 3D loading path. Compared to fabric tensor of weak contact system, the fabric tensor of strong contact system reflects better the characteristics of macroscopic stress tensor. Fabric tensors of strong and weak contact systems contribute differently to the granular macroscopic response. To divide the strong and weak contact system, there is a range for the threshold, however adopting the average contact force is relatively simple and reasonable.

2019, 51(1): 26-35. doi: 10.6052/0459-1879-18-338

HTML PDF

STUDY ON LOW-SPEED WATER ENTRY OF SUPER-HYDROPHOBIC SMALL SPHERES

Huang Chao, Wen Xi, Liu Moubin

Water-entry ofobject is a type of fluid-structure interaction problem withcomplex physics, which has a wide range of engineeringapplications. During the water-entry process, as the object passesthe water-air interface, the possible air entrapment may formcavity under specific conditions and the dynamic movement ofcavity can generate a jet pointing towards the object. Both thecavity and jet can greatly influence the force on the object andtherefore influence its movement. Moreover, the super-hydrophobicsurface can form multi-scale fluid-solid coupling phenomena in theprocess of water-entry of object, which in turn affects the motionof the object and the macroscopic fluid flow. For small scalewater-entry problem, the surface or interfacial force rather thanbody force is usually dominant. The objective of present study isto obtain the cavity behaviors and motion of characteristics forthe problem of water-entry of small super-hydrophobic spheres in awider parameter space. Based on such purposes, the high-speedphotography experiment method is used to study the low-speedwater-entry and cavity dynamics of super-hydrophobic spheres witha radius of 0.175$\sim $10 mm. Five dynamic modes are obtainedincluding sphere floating oscillation, quasi-static impact cavity,shallow seal impact cavity, deep seal impact cavity and surfaceseal impact cavity. The relationship between these dynamic modeswith the Weber number and Bond number is discussed, and thedimensionless relationship between the floating oscillation andsinking behaviors of the sphere is derived. The results show thatthe cavity behaviors of water-entry of super-hydrophobic spheresare mainly related to the Weber number and the Bond number. In therange of $BoO(1)$, the shallow seal state also does not happen. Thecritical relationship between the floating oscillation and thesinking of sphere can be described by the scaling laws. A formulais developed to categorize the floating oscillation andquasi-static impact cavity phenomena.

2019, 51(1): 36-45. doi: 10.6052/0459-1879-18-310

HTML PDF

THE COMPOSITE-MAPPING HYBRID ALGORITHM AND ITS APPLICATIONS OF VIBRATION PLATE BURIED IN PARTICLES

Pan Jingwu, Li Jian, Hong Guangyang, Li Hongying

The interaction between structure and granular materials exists widely in various engineering fields, and the research method of this continuum-discrete coupling problems faces numerous challenges. The composite-mapping hybrid algorithm is presented to research dynamics of continuum-discrete coupling problems. The continuum model is divided into the inner region and the border region of particle contact. In the border region, the composite spheres method is applied to construct the profile of continuum efficiently in order to facilitate fast contact detection between the continuum and particles. In the inner region, the finite element mapping method is introduced to precisely calculate the internal force and deformation of continuum, and the method also contains Rayleigh damping mapping processes. The program with the composite-mapping hybrid algorithm is developed based on the compute cluster and GPU parallel computing technique. The numerical simulation of the square vibration plate which supported at four fixed edges and buried in particles is done to study continuum-discrete coupling dynamics problems. The results show that the proposed composite-mapping hybrid algorithm is appropriate for realization of the compute cluster and GPU paralleled computing technique and improvement of computational efficiency. In analysis on buried plate problems, motion and deformation of the plate can be easily and accurately measured by means of the algorithm. Simultaneously, contact detection can be achieved rapidly in the interface between continuum and discrete, and mechanical parameters of displacement, deformation and vibration modes can also be calculated. The influence of excitation frequency and amplitude on square plate's nonlinear vibration has been studied through excitation with constant amplitude-changing frequency and with constant frequency-changing amplitude, and the period-doubling has been found. Meanwhile, the energy dissipation of granular media in this continuum-discrete coupling system is provided.

2019, 51(1): 46-55. doi: 10.6052/0459-1879-18-343

HTML PDF

DETERMINATION OF DISCRETE ELEMENT MODEL CONTACT PARAMETERS OF NYLON POWDER AT SLS PREHEATING TEMPERATURE AND ITS FLOW CHARATERISTICS

Tan Yuanqiang, Xiao iangwu, Zhang Jiangtao, Jiang Shengqiang

Nylon powder is a commonly used powder material in Additive Manufacturing whose fluidity is closely related to temperature. Exploring powder fluidity at preheating temperature in Additive Manufacturing is the basis for studying the fluidity and spreading properties of powder in selective laser sintering (SLS) process. Choosing nylon powder in SLS technology as a raw material and the flow behavior of nylon powder is studied by discrete element method (DEM), which is a hot topic of numerical simulation and powder spreading process optimization in Additive Manufacturing. Based on Hertz-Mindlin model, Hamaker theory model and Coulomb's law, Van der Waals and electrostatic force are introduced to describe the contact dynamics of nylon powder at preheating temperature. The DEM model of nylon powder at preheating temperature was established based on the mechanical parameters and the rationality of the model was verified by comparing with the experimental results. The flow process of nylon powder in a heated rotating roller was simulated by DEM which checked the correctness of the model. The effects of particle size and particle size distribution on the flow characteristics of nylon powder were studied. The results show that the adhesion force of nylon powder is the result of the interaction of electrostatic force and van der Waals force. With the increase of particle size, the collapse angle of nylon powder decreases and the fluidity of nylon powder increases. And the nylon powder fluidity with uniform particle size distribution is stronger than that of Gaussian particle size distribution. The results can guide the optimization of powder spreading process in SLS.

2019, 51(1): 56-63. doi: 10.6052/0459-1879-18-341

HTML PDF

SEGREGATION BEHAVIOR OF SINTER IN VERTICALLY ARRANGED COOLER WITH HIGE PERFORMANCE GPU SIMULATION

Zhang Xuekuan, Xu Jiy, Sun Junjie, Zhang Yongjie, Zhang Zhenghao, Ge Weiy

Energy-saving and emission-reduction technologies are increasingly required in the iron and steel industry, leading to urgent demanding for very efficient methods of the waste heat recovery and dust emission reduction. The vertically arranged sinter cooler is a new and an efficient apparatus to recover the sensible heat and reduce the dust pollution in the sintering process, which attracts much more attention in recent years. However, the segregation phenomenon is very severe in current design due to the wide diameter distribution of the sinter particles, leading to great reduction of the heat recovery. In order to solve this problem, the structure of the vertically arranged sinter cooler and the operating conditions should be optimized. However, it is very hard to obtain the detailed information of the distribution of the sinter particles in an industrial-scale apparatus. Along with the development in the computer science, the discrete element method (DEM) could provide more and more power for the study of particulate systems, which obtains detailed information of the particles. Thus, DEM is adopted to study the segregation of sinter particles in the vertically arranged sinter cooler. To alleviate the problem of huge computing load, the graphics processing unit (GPU) is adopted to accelerate the DEM simulation. It is found that the inlet tube has significant influence on the distribution of the sinter particles, so that three types of feeding tube structures are designed and tested. Although the sinter particles of different diameters are evenly mixed originally, severe particle segregation occurs in the sinter cooler, where the small and large sinter particles are mostly located at the center and in the periphery regions, respectively. It is obvious that the level of segregation changes with the structure of the feeding tube, which shows that both the number and inclined angle of the feeding tubes will affect the final segregation. The results show that four inlet tubes with small inclined angles are better for tailoring the size distribution of the sinters. So optimization of the structure of the inlet tubes could reduce the segregation of the sinter particles and the efficiency of sensible heat recovery will be improved accordingly.

2019, 51(1): 64-73. doi: 10.6052/0459-1879-18-339

HTML PDF

INFLUENCE OF WATER LEVEL ON ICE LOAD ON UPWARD-DOWNWARD CONICAL STRUCTURE BASED ON DEM ANALYSIS

Long Xue, Liu Shewen, Ji Shunying

During the interaction between the sea ice and the conical structure, the ice load is affected by the position of the sea ice acting on the conical structure when the tidal water level changes. In this study, the discrete element method (DEM) with the bond and failure model is adopted to simulate the breaking process of sea ice acting on the conical structure. In DEM simulations, the influence of the ice temperature at the top and bottom surface of ice cover on the ice strength is considered. The ice load and the failure mode of ice cover simulated with DEM are compared well with the field data in the Bohai Sea. The DEM results indicate that the bending failure of ice cover occurs when it acts on the upward or downward cone. For both of the upward and downward cone, the ice load increases with the increase of the cone diameter at water level, but the ice load on upward cone is larger than that of downward cone. Meanwhile, the breaking length of ice cover on upward cone is smaller than that on the downward cone. The reasons for the difference of ice load and breaking length between upward cone and downward cone are analyzed based on the DEM results and the filed observation in the Bohai Sea. Moreover, when the sea ice interacts with the interface of upward and downward cone, bending damage generally occurs. When the height of the ice cover center is very close to the height of the upward-downward cone interface, the local ice crushing occurs, but the ice load does not increase significantly. Therefore, the downward conical structure has better performance of ice resistance with effectively reducing the ice load. The DEM can be adopted to understand the mechanism of ice cover failure process, and provide reference for the anti-ice design of structures in cold ocean engineering.

2019, 51(1): 74-84. doi: 10.6052/0459-1879-18-342

HTML PDF

SHOCK-TUBE EXPERIMENTAL STUDY AND KINETIC MODELING OF JP-10 PYROLYSIS

Xiong Zhuang, Wang Su, Zhang Can, Yu Hongru

Pyrolysis of JP-10 hydrocarbon fuel was studied in a single-pulse shock tube over temperature range from 1150 K to 1300 K. The main decomposition products were identified by gas chromatography as ethylene, acetylene, propylene, $n$-butene, 1,3-butadiene, cyclopentadiene, cyclopentene, benzene, toluene, and a small amount of methane, ethane, xylene and 1-methylcyclopentene. By summation of all product concentrations in each run, the rate coefficient of JP-10 pyrolysis was experimentally determined. Comparative rate measurements were used to eliminate the effects of shock's non-ideality and boundary layer. A small amount of the internal standard compound, whose rate expression for decomposition is well established, was added in the test gas mixtures, and the reaction temperatures were determined according to the decomposition extents of the internal standard compound under the same experimental conditions in a shock tube. The reaction temperatures determined from the decomposition extents of the internal standard compound are usually less than those at the region 5 behind reflected shock calculated by shock velocity measurements. The temperatures determined by two methods are consistent between 1150 K and 1300 K, the difference is within 20 K, and the difference increases with the temperature increase. Based on the experimental study, kinetic modeling of JP-10 pyrolysis was carried out according to San Diego Mechanism. The yields of three main products, ethylene, acetylene and 1,3-butadiene, have a good agreement between the experimental and the simulation results,while the experimental results of cyclopentene yield are much higher than simulation,indicating that both fully and partially ring-opening reactions in JP-10 pyrolysis are important decomposition reaction pathways.

2019, 51(1): 85-93. doi: 10.6052/0459-1879-18-102

HTML PDF

NUMERICAL STUDY ON THE AERODYNAMIC PERFORMANCE OF THEFLEXIBLE AND CORRUGATED FOREWING OF DRAGONFLY IN GILDINGFLIGHT

Liu Huixiang, He Guoyi, Wang Qi

Dragonflies are capable of carrying out dramatic flight manoeuvres, gliding flight is a common mode of flight for dragonfly, and dragonfly wings are the source of dragonflies dramatic flight manoeuvres. Unlike typical engineered airfoil, dragonfly wings are not smooth, wing cross-section are highly corrugated. It has been shown that corrugations could enhance the spanwise stiffness in the wings and influence aerodynamic performance of the dragonfly wings. Flexibility is another characteristic of the dragonfly's wings, which is mainly manifested as the flexible deformation of the wings during the flight. To explore corrugations and flexibility effect on aerodynamic performance of the dragonfly forewings in gliding flight, a computational fluid dynamics (CFD) model and a computational structural mechanics (CSD) model of the corrugated dragonfly forewing are established based on current research, and the modal analysis verified that the model has sufficient accuracy. The corrugated rigid and flexible dragonfly forewing are acquired by using CFD method and CFD/CSD coupling method respectively. The simulation indicated that flexible and corrugated forewings is subjected to aerodynamic load, which only produces bending deformation without torsion deformation in gliding flight, and the aerodynamic response time is short. Compared with the aerodynamic performance of the rigid and corrugated forewings, the result showed that veins and cuticular membrane of flexible forewing are deformed which caused the lift coefficient and drag coefficient decrease, the leading edge vortex of flexible forewing is much higher than rigid forewings at large angle of attack because of deformed vein. The aerodynamic performance of rigid forewings is better below 10 degree angle of attack,and the aerodynamic performance of flexible forewings is better at large angle of attack as result.

2019, 51(1): 94-102. doi: 10.6052/0459-1879-18-157

HTML PDF

PIV EXPERIMENTAL INVESTIGATION ON THE BEHAVIOR OF PARTICLES IN THE TURBULENT BOUNDARY LAYER

Gao Ti, a, Sun Jiao, Fan Ying, Chen Wenyi

The particle image velocimetry (PIV) is used to conduct experimental research in the solid-liquid two-phase plane turbulent boundary layer. The turbulence statistics such as the average velocity profile, turbulence intensity and Reynolds stress of the particle phase and single-phase clean water are compared to analyze the behavior of the particles in the turbulent boundary layer. The concept of multi-scale spatial locally-averaged vortices is utilized to extract the spatial topologies of the spanwise vortex head and the statistic of the prograde vortex is acquired. From that, the spatial topologies of the fluctuating velocity and streamlines around the prograde vortex at different normal positions can be obtained. The degree of development of the prograde vortex and the surrounding turbulence coherence structure can be compared and analyzed. The results show that compared with the clean water conditions, the buffer layer of the turbulent boundary layer of the particle phase becomes thinner, the logarithmic region moves downward, the turbulence intensity is enhanced, and the Reynolds stress in the logarithmic law region is increased. The fluctuating velocity of the particle phase is different from the clear water condition around the vortex, and the particles can be effectively transferred by the burst process around the spanwise vortex. The particle-laden flow has a large prograde vortex core and develops as the normal position rises. The vortex and the band stretch longer in the flow direction. At the same time, it is found that there is always a retrograde vortex in the lower left of the prograde vortex under both conditions, and the formation of retrograde vortex in the particle-laden flow is weaker than that of single-phase fluids. The number of prograde vortices in both conditions decreases with the increase of the normal position, and finally gradually stabilizes.

2019, 51(1): 103-110. doi: 10.6052/0459-1879-18-211

HTML PDF

EVOLUTION CHARACTERISTIC AND WORKING MECHANISM ANALYSIS OF ROTATING THIN-WALLED STRUCTURES IN POST-CRITICAL TURBULENT INTERVAL BASED ON FIELD MEASUREMENT

Wang Hao, Ke Shitang

Recent study found that the time-varying characteristic of the load may have a significant effect on the vibrational strength and energy mechanism. The most important structures in fire/nuclear power plants (such as cooling towers, chimneys, etc.) are all typical rotating thin-walled structures. To reveal the vibration evolution characteristic and working mechanism of thin-walled structures in post-critical turbulent interval, the vibration responses of eight typical rotating thin-walled structures of high Reynolds number flow ($Re \ge $3.5$\times $10$^6$) are measured. Firstly, non-stationary identifications of signals with different time intervals are performed after depressing and filtering noise. The time-varying mean and extreme estimation of response are studied based on non-stationary analysis model. Besides, the frequency domain evolution characteristics are studied based on evolution spectrum method. On this basis, proportion of resonance component in wind-induced response and its effect are discussed. Then, self-resonant frequency and damping ratio of the structures are identified, and the damping mechanism of different rotating thin-walled structures is studied. The evolution characteristic and working mechanism are revealed as follows. (1) The wind-induced vibration response of the rotating thin-shell structure in post-critical turbulent interval is characterized by stable frequency evolution characteristics and non-stationary evolution characteristics in intensity aspect; (2) The wind-induced vibration problem of rotating thin-walled structures in post-critical turbulent interval should be studied as quasi-static and resonance excitation points separately. The vibration energy distributions of resonance excitation points at different regions of the cooling tower were similar, but the PSD functions of quasi-static points were dramatically different from each other; (3) Vibration energy distribution of the resonant excitation points showed a phased trend, and the proposed resonance spectral expression takes three variation stages of responses into account and achieves high prediction accuracy; (4) With the concept of equivalent damping ratio proposed in this paper, the damping ratio prediction formula is proposed. More importantly, these analysis results show that resonance effects and non-stationary effects on wind-induced effects of rotating thin-walled structures in post-critical turbulent interval are generally notable, and the irrationality of 5% damping ratio value commonly used in the current project for this type of rotating thin-walled structure has been demonstrated.

2019, 51(1): 111-123. doi: 10.6052/0459-1879-18-125

HTML PDF

A LATTICE BOLTZMANN SIMULATION OF FLUID FLOW IN POROUS MEDIA USING A MODIFIED BOUNDARY CONDITION

Cheng Zhilin, Ning Zhengfu, Zeng Yan, Wang Qing, Sui Weibo, Zhang Wentong, Ye Hongtao, Chen Zhili

The lattice Boltzmann method has been considered as an effective method for the simulation of hydrodynamic flows. Handling the boundary condition accurately in simulation is extremely essential for a reliable study. In this paper, a multiple relaxation time lattice Boltzmann model with different boundary conditions was applied to mimic the flows in periodically symmetric and irregular structures. The scope of application and accuracy for different boundary conditions in various geometries was investigated. In addition, a hybrid boundary treatment method was introduced to simulate the non-Darcy flow in porous media, the simulation results of which were also compared to the results obtained using pressure boundary condition. The results show that for the symmetric and periodic flow simulation, both the body force and the pressure driven boundary treatments are perfectly equivalent and both can accurately capture the flow characteristics. While for the fluid flow in irregular structures, the body force and pressure boundary conditions are not equivalent, and the body force one has limited use and can only be applied to periodic structures. This implies that one must be cautious of the reliability of modeling when conducting model validation with simple structures. It seems that the regular structures could be inadequate to validate the modeling, which depends on the research issues, i.e., the flow patterns in what kinds of structures. Furthermore, the generalized periodic boundary condition proposed by previous authors combines periodic density momentum with a pressure gradient in one dimension is also not appropriate to conduct flow simulation in irregular models since this method ignores the effect of asymmetric obstacles in the direction perpendicular to the main streamlines. Moreover, the hybrid boundary condition can be used to perform flow simulations not only in periodic structures but also the irregular ones. In particular, for the inertial flow of fluids in porous media, the relatively high Reynolds number can be achieved readily with the hybrid boundary condition. For the pressure driven boundary condition, the pressure gradient comes from the density difference between the inlet and outlet. To provide a higher Reynolds number, it is necessary to implement a great density contrast in inlet and outlet nodes. However, this approach is inconsistent with physical situation and causes undesirable errors in simulation. All in all, the hybrid boundary condition has greater advantages over the pressure boundary condition.

2019, 51(1): 124-134. doi: 10.6052/0459-1879-18-179

HTML PDF

CIRCULATION CONTROL ON WIND TURBINE AIRFOIL WITH BLUNT TRAILING EDGE

Qiao Chenliang, Xu Heyong, Ye Zhengyin

The blunt trailing edge wind turbine airfoil has the advantages of high structural strength and insensitive to surface contamination, but its larger drag coefficient makes the aerodynamic characteristics of the blunt trailing edge wind turbine airfoil unsatisfactory. The circulation control method is implemented on a blunt trailing edge wind turbine airfoil in order to improve the aerodynamic characteristics of the blunt trailing edge wind turbine airfoil, and weaken the strength of the shed vortex. The circulation control on a blunt trailing edge wind turbine airfoil is investigated using numerical simulation methods. The effectiveness of the circulation control method in increasing lift and reducing drag is researched. The variation of lift and drag coefficients with the jet momentum coefficient is studied, and the aerodynamic figure of merit and the control efficiency of the circulation control with different jet momentum coefficients is analyzed. The results show that the circulation control can significantly enhance the lift coefficient of the blunt trailing edge airfoil, and effectively reduce the drag coefficient of the airfoil; the lift coefficient of the airfoil increases with the jet momentum coefficient, and two typical control regimes, namely, separation control and super-circulation control, are observed; the power coefficient of jet increases with increasing jet momentum coefficient, and growth rate gradually increases; the output power of blade also increases with increasing jet momentum coefficient, but the growth rate decreases gradually. It is demonstrated that the circulation control method can significantly improve the aerodynamic characteristics and power output characteristics of blunt trailing edge wind turbine airfoils, and it has a good application prospect in large wind turbine flow control.

2019, 51(1): 135-145. doi: 10.6052/0459-1879-18-164

HTML PDF

FREE VIBRATION ANALYSIS OF THIN-WALLED AXISYMMETRIC STRUCTURES WITH BOUNDARY ELEMENT METHOD

Zhou Qi, Chen Yongqiang

The dual reciprocity method(DRM) is extended to study the eigenvalue and eigenmode of thin-walled axisymmetric structures. First the displacement in the domain integral can be approximated by a set of radial basis functions and the domain integral can be converted to the boundary using DRM. Then the displacement and the traction can be expanded as Fourier series and integrate along the circumferential direction. The obtained boundary integral equation can be used for analysis of elastostatics of axisymmetric structure distributed body force and elastodynamics subject to asymmetric loading. The special case of the source point on the axis of symmetry is discussed in detail. New schemes are suggested for dealing with singular matrices for cases with and without body force respectively according to the degenerate form of the fundamental solution and the particular solution. For the thin walled structure, the sinh transformation is applied to improve the accuracy of evaluation of the nearly singular integrals. The developed project has been used to analyze elastostatics with body force and the free vibration of the thin axisymmetric structures. Numerical results indicate that the proposed method for dealing with singular matrices can effectively deal with the situation where the source point is on the axis of symmetry. and when the thickness ratio reaches $10^{-3}$, the relative error of the results can approach $10^{-3}$, which is better than those of FEM.

2019, 51(1): 146-158. doi: 10.6052/0459-1879-18-251

HTML PDF

SPHERICAL INDENTATION METHOD TO DETERMINE STRESS-STRAIN RELATIONS AND TENSILE STRENGTH OF METALLIC MATERIALS

Zhang Zhijie, Cai Lixun, Chen Hui, Bao Chen, Liu Xiaokun

It is significant to obtain the uniaxial stress-strain relations and tensile strength of materials by indentation method for evaluating the integrity of structures in service. Based on energy equivalent assumptions, i.e., the Von Mises equivalence about a RVE (representative volume element) and the median energy equivalence in the effective deformation region, a semi-analytical spherical indentation(SSI) model is proposed to describe the relation among indented load, depth, the diameter of spherical indenter and Hollomon-law parameters if indented materials are uniform, continuous, isotropic and power-law hardening. The stress-strain relations and tensile strengths of the materials were obtained by the tested load-depth curves under spherical indenter loading. Considering the damage effect in the indentation process, a correction model for determining elastic modulus of the metallic materials due to the spherical indentation tests is proposed . For eleven kinds of tested ductile metallic materials with spherical indentation, the Young's modulus, stress-strain relations and tensile strength predicted by SSI model are in good agreement with the uniaxial tensile results.

2019, 51(1): 159-169. doi: 10.6052/0459-1879-18-200

HTML PDF

ANALYTICAL STRESS SOLUTION RESEARCH ON AN INFINTE PLATE CONTIANING TWO NON-CIRCULAR HOLES

Zeng Xiangtai

The analytical stress solution for an infinite plate containing a single hole of arbitrary shape can be obtained by complex variable method. As to the doubly-connected domain problem that an infinite plate contains two round holes or two elliptical holes, it can also be solved using a variety of methods, such as the bi-polar coordinate method, the stress function method, the complex variable method, and the Schwarz alternating method. The complex variable method combined with conformal mapping is one of importance methods which can be used to obtain analytical stress solution, but it is not yet used to solve the problem of an infinite plate containing a square hole and an elliptical hole. Taking advantage of the conformal mapping method, the problem that an infinite plate contains two specific holes, which far-field uniform stress is applied at infinity and the boundaries of the two holes are subjected to uniform vertical compression, can be solved. The key step of this method is to find the corresponding mapping function with which the considered region can be mapped onto a ring in the image plane. Based on the Riemann mapping theorem, we propose a general form of the mapping function and figure out the concrete mapping function for the specific problem using optimization method. The basic equation set for solving the two analytical functions is established through the stress boundary condition of the two holes. Then the analytical stress solution can be obtained according to the two analytical functions. The analytical stress solution is compared with numerical stress solution of ANSYS finite element method. Effects of separation distance, size of elliptical hole, and the orientation of holes on tangential stress of the boundary is investigated using the newly derived solution. The stress distributions on the line that connects the centers of the two holes under different loads are presented.

2019, 51(1): 170-181. doi: 10.6052/0459-1879-18-197

HTML PDF

A LARGE DEFORMATION ELASTOPLASTIC MODEL BASED ON THE INTERMEDIATE CONFIGURATION

Meng Chunyu, Tang Zhengjun, Chen Mingxiang

In the large deformation elastoplastic constitutive theory, a basic problem is the decomposition of elastic deformation and plastic deformation. In the usual case, two decomposition methods are adopted. One method is to decompose deformation rate (or strain rate) into elastic and plastic parts. Among them, the elastic deformation rate and the objective rate of Kirchhoff stress are linked by elastic tensors, and by using this way, a sub-elastic model can be established. In the mean time,the plastic deformation rate is related to Kirchhoff stress by using flow law. Another method is to decompose the deformation gradient tensor based in the intermediate configuration. It is supposed that by considering the virtual unloading process, an unstressed intermediate configuration can be obtained and the so-called hyperelastic-plastic model can be established. In this paper, a large number of properties of a large deformation elastoplastic model which is based on the deformation gradient multiplicative decomposition are studied, and the model is built in the intermediate configuration. These properties include: in different configurations, the existence of the plastic spin rate; the symmetry of back stress; the orthogonality of plastic deformation rate and yield surface; and the relationships between plastic spin rate, back stress, plastic deformation rate and yield surface. First of all, by using the tensor function representation theorem which is in the appendix, some special properties of the isotropic function are obtained and some formula relationships are established, and some simple relationships of tensor value function between the intermediate configuration and the current configuration are derived. Secondly, based on these properties and relationships, and in combination with the laws of thermodynamics, the mathematical expression of the model in different configurations is established, which consists of objective rate representation and continuous tangential stiffness. Thus, some properties of the large deformation elastoplastic model based on the intermediate configuration are obtained. Finally, the model is compared and analyzed with the four models.

2019, 51(1): 182-191. doi: 10.6052/0459-1879-18-138

HTML PDF

PERFORMANCE OF COLLINEAR THERMOELECTRIC GENERATOR CONSIDERING THE HEAT DISSIPATION IN THE SIDE SURFACE

Wang Xueqiang, Ju Chengjian, Dui Guansuo

Thermoelectric material is an environment-friendly function material, which can~convert energy between heat and electricity. And it holds extensive~application~potentiality in power generation and refrigeration. The traditional thermoelectric generator is a $\pi $-type structure, which requires the length of the thermoelectric legs to be equal. In some cases, the structure is not conducive to the optimal design of the thermoelectric generator. Intense thermal stress and even stress~concentration will be induced in the thermoelectric generator due to high-temperature working condition, leading to shortening its working life. In addition, since the operating temperature of the thermoelectric generator is higher than ambient temperature, part of the heat will inevitably be dissipated to the environment, which will affect the thermoelectric performance and mechanical performance of the thermoelectric generator. Therefore, the heat dissipation cannot be neglected when analyzing this kind of problem. For these phenomena, in this work, a novel collinear-type thermoelectric generator model is proposed considering the heat dissipation in the side surface. And the legs of the proposed thermoelectric model can be optimized independently. Then, based on the finite element method, performance of collinear thermoelectric generator considering the heat dissipation in the side surface is simulated. And the thermoelectric performance and mechanical performance under the Dirichlet boundary condition is analyzed. Simultaneously, the temperature field, electric potential field and stress field in the thermoelectric generator are obtained. The influence of various convective heat transfer coefficient on thermoelectric performance and mechanical performance of the thermoelectric generator is investigated. The results demonstrate that thermal convection can decrease the energy conversion efficiency of the thermoelectric generator. When the convective heat transfer coefficient reaches 100 W/(m$^{2}\cdot$\textcelsius), the efficiency is 0.047 9 which is 28% lower than the conversion efficiency of 0.066 7 in adiabatic state. Though heat loss from the side surface is increased due to heat convection, thermal stress is reduced. In practical application, proper design and improvement of the thermal insulation system should be carried out to improve the efficiency of energy conversion.

2019, 51(1): 192-197. doi: 10.6052/0459-1879-18-255

HTML PDF

TRAJECTORY SIMULATION OF SELF-PROPELLED ELASTIC RODS IN FLUID

Wang Jun, Zhu Yongning, Xu Jian

The research for flexible structures coupling with fluid simulates and promotes the development of soft robotics. A fast and accurate numerical method is significant for a real-time simulation of robots. This research provides theoretical and critical information for experiments to reduce the possibility of failure, by anticipating the path of the underwater soft robots and the possible requiring parameters for materials. This paper researches for bio-swimming elastic rods coupling with two-dimensional incompressible Newton fluid. First, we discretize elastic rods to extensive springs and rotational springs and stablish the kinetic equations based on energy reflecting the influence of internal force on swimming rods, and solve the governing equation by leap-frog algorithm. Second, semi-Lagrangian method is used to establish a fluid solver. Finally, the simplified coupling algorithm based on immersed-boundary method is raised, calling immersed-boundary method for momentum equations. These equations update the velocity of grid near the coupling interface directly to replace the function source force play in immersed-boundary method. Combining the fluid solver, rods solver and the coupling algorithm, a integral program solving the problem of the dynamics for immersed rods numerically. Simulating the rods swimming with a referenced pose of sine curvature. Comparing the simulating result to existed experiments of filament swimming and to the swimming trajectory of soft-body fish, we find the numerical result conform to experiment result, leading to the conclusion that this algorithm and dynamic model can simulate the trajectory of discrete elastic rods swimming underwater smoothly under the influence of both internal elastic force and coupling soft body-fluid interaction. Using this program, we test several critical parameters relating to swimming performance of rods, including iteration numbers, frequency and initial phases, finding that changing the initial phase of rods will alter onward direction of elastic rods. These results prove the feasibility and possibility of this algorithm and program being used for guiding development of real soft swimmers.

2019, 51(1): 198-208. doi: 10.6052/0459-1879-18-230

HTML PDF

A NUMERICAL METHOD FOR DYNAMICS OF PLANAR MULTI- RIGID-BODY SYSTEM WITH FRICTIONAL TRANSLATIONAL JOINTS BASED ON LUGRE FRICTION MODEL

Wang Xiaojun, Lü Jing, Wang Qi

A numerical method for the dynamics of the planar multi-rigid-body system with frictional translational joints is presented in this paper. The multibody system consists of the several rigid bodies which are linked with ideal revolute joints and imperfect translational joints. The frictional forces on the slider in the imperfect translational joint are modeled by the LuGre friction law which can effectively describe stick-slip motion in the mechanical system. The sizes of the clearances and the impacts between the guide and the slider in the translational joints can be neglected when the clearance sizes are very small, so the geometric constraints of the translational joints are treated as bilateral constraints. In this work, firstly, the complementarity conditions and formulations about the normal forces on the slider in the translational joint are given. Secondly, the dynamical equations of the multibody system are obtained by the Lagrange's equations of the first kind and the Baumgarte stabilization method for the constrained multibody systems in order to reduce the constraint drift in the numerical simulation of the multibody systems. Thirdly, the problems of determining contact situations of the slider in the translational joint and solving normal forces on the slider in all contact situations are formulated and solved as a linear complementarity problem (LCP). Finally, two numerical examples of the planar multi-rigid-body system with a frictional translational joint are given to illustrate their dynamical behaviors such as stick-slip motion and several contact situations of the slider in the translational joint. The LuGre friction model and the Coulomb friction model are used in two numerical examples to compare the dynamical behaviors of two mechanical systems. The numerical results obtained by our method are compared with that obtained by other method. The numerical results of the examples show the availability of the method presented in this paper.

2019, 51(1): 209-217. doi: 10.6052/0459-1879-18-222

HTML PDF

NONLINEAR VIBRATIONS OF AXIALLY MOVING BEAMS WITH NONHOMOGENEOUS BOUNDARY CONDITIONS

Zhang Dengbo, Tang Youqi, Chen Liqun

The transverse nonlinear vibration of the axial moving system has become one of the hot subjects at home and abroad. At present, most of the related studies are on the homogeneous boundary conditions. However, the nonhomogeneous boundary conditions are more common in the engineering practice. There are relatively few publications on the axial moving system with homogeneous boundary conditions. In order to study the effect of nonhomogeneous boundary conditions on the transverse nonlinear vibration of an axial moving beam, the nonlinear parametric vibrations of an axially accelerating viscoelastic Euler beam under nonhomogeneous boundary conditions are studied in this paper. The variable tension caused by the axial acceleration is introduced. A nonlinear integro-partial-differential equation and corresponding nonhomogeneous boundary conditions of an axially accelerating viscoelastic beam are presented. The effects of nonhomogeneous boundary conditions are highlighted. The method of multiple scales is used to establish the solvability conditions. The steady-state response of the beam was obtained from the solvability condition. According to the Routh-Hurvitz criterion, the stability of the response was determined. Some numerical examples are introduced to demonstrate the effect of the viscoelastic coefficient, mean speed, axial speed fluctuation amplitude, nonlinear coefficient, and the nonhomogeneous boundary conditions on the steady-state response. The larger viscoelastic coefficient leads to the smaller instability interval of the trivial solutions and the smaller stable steady-state response amplitude especially for the second mode. The instability interval and the stable steady-state response amplitude with the nonhomogeneous boundary conditions are larger than that with the homogeneous boundary conditions. The differential quadrature scheme is introduced to confirm the approximate analytical results. The numerical results show reasonable agreement with the approximate analytical results.

2019, 51(1): 218-227. doi: 10.6052/0459-1879-18-189

HTML PDF

SERIES-MODE PITCHFORK-HYSTERESIS BURSTING OSCILLATIONS AND THEIR DYNAMICAL MECHANISMS

Zhang Yi, Han Xiujing, Bi Qinsheng

Bursting oscillations is a spontaneous physical phenomenon existing in natural science, which has various patterns according to their dynamical regimes. For instance, bursting of point-point type means bursting patterns related to transition behaviors among different equilibrium attractors. Pitchfork-hysteresis bursting, induced by delayed pitchfork bifurcation, is a kind of point-point type bursting pattern showing simple dynamical characteristics. The present paper takes the Duffing system with multiple-frequency parametric excitations as an example in order to reveal bursting patterns, related to delayed pitchfork bifurcation, showing complex characteristics, i.e., the series-mode pitchfork-hysteresis bursting oscillations. We considered the case when one excitation frequency is an integer multiple of the other, obtained the fast subsystem and the slow variable of the Duffing system by frequency-transformation fast-slow analysis, and analyzed bifurcation behaviors of the fast subsystem. Our study shows that two or multiple pitchfork bifurcation points can be observed in the fast subsystem, and thus two or multiple pitchfork-hysteresis bursting patterns are created when the slow variable passes through these points. In particular, the pitchfork-hysteresis bursting patterns are connected in series, and as a result, the so-called series-mode pitchfork-hysteresis bursting oscillations are generated. Besides, the effects of parameters on the series-mode pitchfork-hysteresis bursting oscillations are analyzed. It is found that the damping of the system and the maximum excitation amplitude show no qualitative impact on corresponding dynamical mechanisms, while the smaller one may lead to vanish of busting oscillations. Our findings reveal the road from simple dynamical characteristics of point-point type bursting oscillation related to complex one, thereout, a complement and expansion for nowadays bursting dynamics is obtained.

2019, 51(1): 228-236. doi: 10.6052/0459-1879-18-223

HTML PDF

STRUCTURAL DYNAMIC AND STABILITY ANALYSIS OF A STRIPPED SOLAR SAIL

Zhang Junhui, Cui Yangyang, Tong An

Propelled by solar pressure, solar sail is thought as the most promising interstellar exploration technology. Two kinds of solar sail architecture have been proposed: boom-supporting solar sails and rotating solar sails, among them, stripped solar sail whose membrane is divided into separate narrow membrane strips is the ideal architecture for boom-supporting solar sails. How to precisely calculate structural dynamic characteristics of a stripped sails is worth studying. In this paper, the dynamic characteristics and stability of a stripped solar sail is studied. The whole solar sail structure is regarded as an assembly consisting of several sequentially connected boom-strip components, and one boom-strip component includes four boom segments and four membrane strips. In the stripped solar sail, the booms are under multiple-step axial loads and their vibrations are coupled with the vibrations of stripped membranes. Considering the coupling effects between booms and membrane strips, a closed form vibration model of the stripped solar sail is established by the distributed transfer function method. Based on this model, vibration characteristics and buckling loads of the stripped solar sail can be determined accurately and efficiently. Numerical results indicate that stripped solar sail architecture benefits to enhancing structural stiffness and stability, the more membrane strips will provide the higher stiffness and buckling loads, increasing membrane pre-stress will decrease base frequency of the stripped solar sail and worsen the structural stability, free vibration frequencies and buckling loads will increase with stiffness of the supporting booms. This presented analytical method is more efficiency and accurate than the numerical methods, which provides an efficient analysis tool for structure design and attitude control of stripped solar sails.

2019, 51(1): 237-244. doi: 10.6052/0459-1879-18-275

HTML PDF

WIND TUNNEL INVESTIGATION OF THE AERODYNAMIC CHARACTERISTICS OF PURPLE WISTERIA COMPOUND LEAVES

Yu Kejie, Shao Chuanping

The capability of reconfiguration of tree leaves is of significance in the designs of solar panel, aerofoil, bionic antenna, and wind power generation tree. The leaf was clamped at the base end of the rachis and vertically suspended in the center of wind tunnel test section, and tested with its front and back surface facing on-coming stream respectively at step-by-step increasing wind speed from 0 to 25 m/s. Results show that the changing process of the leaf can be divided into three stages: earlier steady, intermediate transition, and later steady, and critical wind speeds are observed. In earlier stage, the downstream bending curvature of the rachis increases rapidly with wind speed, and multi-layer wing steady and multi-layer multi-shape steady states exist. In intermediate stage, large amplitude low frequency rachis vibration, and small amplitude high frequency lobules vibration are observed. In later steady stage, two-layer structure or single streamlined body of conic, or wedge or U-shape cross section can be found. As wind speed increasing, the number of lobule layers and the width of the compound leaf $b$ decrease, until the single streamlined body formed. As $Re$ increasing, the drag coefficient of the leaf decreases rapidly at first, then slowly approaching to a constant. The absolute value of the negative Vogel component $\vert \alpha \vert$ decreases with the increase of lobule number of the leaf. $\vert \alpha \vert $ of the leaf with its back surface is larger than that with its front surface facing wind, but they tend to converge with the increase of lobule number. Rachis vibration occurs if the frequency of vortex shedding from the leaf is close to the natural frequency of the rachis. The second critical wind speed $V_2/\sqrt{E/\rho}$, at which the rachis vibration begins, is shown to be the function of $b/l$ and $d/l$, where $E, \rho , d, $and $l$ are respectively elastic module, mass density, diameter and length of the rachis, $b$ is width of the deformed leaf, and a figure about this function is drawn using experimental data.

2019, 51(1): 245-262. doi: 10.6052/0459-1879-18-140

HTML PDF

A GFEM WITH C$^1$ CONTINUITY

Tian Rong

The function $f$ is said to be of class C$^{1}$ if the first order derivatives of $f$ exist and are continuous. A C$^{1}$ approximate can be applicable, totally up to users' choice, to solve the weak or the strong forms of PDEs, which provides an opportunity on designing a better-fit numerical method. The Partition of Unity Finite Element Method (PUFEM, Babuska and Melenk (1997)) gains broad attention due to a strikingly advantageous feature: A user-tailorablly high order approximation while without complicating numerical implementations in a standard FE code. However, the smoothness of the global approximate function of PUFEM is inherent to that of the partition of unity function that is usually taken as the standard finite element shape function. How to construct a PUFEM of class C$^{1}$ based on the C$^{0}$ finite element shape functions is still a pending problem. Based on the recently developed extra-dof free partition of unity approximation, we develop in this paper a C$^{1}$ continuous generalized finite element approximation using only a C$^{0}$ finite element shape function constructed on a Cartesian grid. The approximation is applied to discretize the Poisson's equation in both strong forms and weak forms. Numerical tests show that the approximation can be applicable to numerical solution to both the strong and the weak form of PDEs, and it is able to deliver a high order of accuracy and convergence without necessarily altering grid topology and increasing nodes. The necessary condition for using the C$^{1}$ approximate is a domain discretization based on a Cartesian grid (not needed to be uniform). The new approximate can be used to both fluids (like FDM) and solids (like material point method). The difference from FDM is that the field function and its derivatives at an arbitrary point of the domain of interest can be computed directly using the "shape function" of the new approximate. When working with the material point method, the new approximate can be expected to reduce or eliminate quadrature errors and cross-grid oscillations.

2019, 51(1): 263-277. doi: 10.6052/0459-1879-18-188

HTML PDF

IMPROVED EXTENDED SCALED BOUNDARY FINITE ELEMENT METHODS

Jiang Shouyan, Li Yun, Du Chengbin

Combining the main advantages of the extended finite element methods (XFEM) and the scaled boundary finite element methods (SBFEM), improved extended scaled boundary finite element methods ($i$XSBFEM) are proposed. The proposed methods can provide a new way for the simulation of fracture problems. Similar to XFEM, two orthogonal level set functions are used to characterize the internal crack surface in materials, and how the element is partitioned by a crack can be judged by level set functions. These elements partitioned by the crack are treated as a subdomain of SBFEM, and then the element stiffness matrix of these discontinuous elements can be directly solved by SBFEM, thus avoiding the need for further element subdivision to the solution of the discontinuous element stiffness matrix in XFEM. At the same time, with the help of the main idea of XFEM, the real displacement of the intersection point between the crack and the element boundary is considered as the additional degrees of freedom of the element nodes, thereby it gives explicit physical meaning of additional degrees of freedom. For the element containing the crack tip, several layers of elements around the crack tip are selected as a super element, and the super element is used as a subdomain of SBFE to solve the stiffness matrix. The node displacement inside the super element can be obtained by the SBFE displacement approximation. The stress intensity factor can be directly obtained based on the singular displacement or stress at the crack tip, without the need of other numerical methods. Finally, several numerical examples are given to verify the effectiveness of the proposed $i$XSBFEM. Compared with the standard XFEM, the relative error convergence of the $i$XSBFEM based on the displacement norm is better, and the stress intensity factors computed by stress based and displacement based method in $i$XSBFEM both are in good agreement with the analytical solution.

2019, 51(1): 278-288. doi: 10.6052/0459-1879-18-218

HTML PDF

2019, 51(1): 289-291.

HTML PDF

REVIEW OF THE 2018 SYMPOSIUM ON APPLICATION OF FLUID-STRUCTURE INTERACTION IN NAVAL ARCHITECTURE AND OFFSHORE RENEWABLE ENERGY

Li Ye, Wang Benlong, Zhan Shige

"The 2018 Symposium on Application of Fluid-Structure Interaction in Naval Architecture and Offshore Renewable Energy" held in January 2018 in Shanghai Jiao Tong University is briefly introduced. The academic reports in the symposium by the invited scientists are reviewed by category. The viewpoint and suggested future focused research directions on the application of fluid-structure interaction in naval architecture and offshore renewable energy in the open discussion session are summarized.

2019, 51(1): 292-297. doi: 10.6052/0459-1879-18-445

HTML PDF

REVIEW OF THE THIRD NATIONAL SYMPOSIUM ON BIOMECHANICS\FOR YOUNG SCHOLARS

Lü Yonggang, Zhan Shige

In this paper, The Third National Symposium on Biomechanics for Young Scholars was briefly introduced, and the academic reports of experts and speakers were summarized and reviewed.

2019, 51(1): 298-303. doi: 10.6052/0459-1879-18-245

HTML PDF

2019 Vol. 51, No. 1

NoticeMORE+

Exquisite ArticlesMORE+

Download Center

Author Center