#### Table of Content

18 September 2019, Volume 51 Issue 5
Research Review
PATHOGENETIC MECHANISMS OF GLAUCOMA------RESEARCH PROCESS ON THE DEFORMATION OF LAMINA CRIBROSA 1)
Zhang Ting, Li Long, Song Fan
2019, 51(5):  1273-1284.  DOI: 10.6052/0459-1879-18-321
Asbtract ( 798 )   HTML ( 81)   PDF (21642KB) ( 460 )
Figures and Tables | References | Related Articles | Metrics

Glaucoma is the first cause of irreversible blinding eye disease in the world. Glaucomatous optic nerve damage is directly associated with the intraocular pressure, and tight control of intraocular pressure is still the only therapeutic approach available for the treatment of glaucoma, while the pathogenesis of glaucoma remains unknown. It has now been confirmed that the primary site of glaucoma is the lamina cribrosa: the pressure difference between the intraocular pressure and intracranial pressure respectively exerted on the anterior and posterior surfaces of lamina cribrosa can cause the change in the structure and morphology of lamina cribrosa, then the deformation of lamina cribrosa squeezes the optic nerves passing through the lamina cribrosa to make their damages; and finally, the damages produce irreversible visual loss. As a result, the pathogenesis of glaucoma is essentially associated with the mechanical properties of lamina cribrosa and mechanical environment surrounding the lamina cribrosa. Since lamina cribrosa was identified as the primary site of glaucomatous optic nerve damage, it has become the hot spot of glaucomatous optic nerve damage research. As an effective method, we can study the deformation of lamina cribrosa under the effect of intraocular pressure and intracranial pressure by developing mechanical model of lamina cribrosa, and analyze the effect of the deformation of lamina cribrosa on the optic nerve damage. This method has helped us to reveal the mechanism of glaucomatous optic nerve damage and the pathogenesis of glaucoma to some extent. This review will introduce the research progress and the present existing problems on the deformation of lamina cribrosa during glaucoma from the related experimental, theoretical, computational and clinical aspects.

Fluid Mechanics
BOUNCING BEHAVIORS OF A BUOYANCY-DRIVEN BUBBLE ON A HORIZONTAL SOLID WALL 1)
Zhang Yang, Chen Ke, You Yunxiang, Sheng Li
2019, 51(5):  1285-1295.  DOI: 10.6052/0459-1879-19-071
Asbtract ( 554 )   HTML ( 41)   PDF (6112KB) ( 353 )
Figures and Tables | References | Related Articles | Metrics

It is not only interesting but also complex that buoyancy-driven bubbles rising in viscous liquids. In particular, the interactions between bubbles and boundaries (e.g., solid walls) is relevant in practical applications and these interactions may have a significant effect in the global behaviors of the multi-phase fluids. In this work, the rising, collision and bouncing to a horizontal solid wall for a single bubble are studied by axisymmetric computations. The incompressible, density-variable Navier-Stokes equations with surface tension are used to describe the gas-liquid flow and are solved by a tree-based finite volume method (FVM). The evolution of bubble shape is implemented by using a volume of fluid (VOF) approach that combines a balanced surface tension force calculation and a height-function curvature estimation. To finely resolve the local but fast topological evolutions of bubble, the technique of adaptive mesh refinement (AMR) is used. Starting with the basic phenomenon of bubble impacting and bouncing, we explore the effects of Galilei number Ga and approach velocity $U_{\rm a}$. To study the bubble behaviors under different conditions both qualitatively and quantitatively, the evolution of the velocity vector field and a lot of parameters such as bouncing height $H$, bouncing period $T,$rising velocity $U$, axial coordinate $z$and coefficient of restitution $C_{\rm r }$are analyzed. Based on the results, we find that the bubble bouncing behaviors are pretty sensitive to the Galilei number. The increase of Ga promotes the bouncing and signifies the deformation of bubbles, increasing the collisions of bubbles, bouncing parameters and the coefficient of restitution. However, the value of $C_{\rm r}$is nearly unaffected by the variation of the approach velocity, indicating that $U_{\rm a}$is not a governing parameter for the bubble bouncing motion.

EXPERIMENTAL INVESTIGATION OF CAVITATION CHARACTERISTICS AND DYNAMICS IN COMPRESSIBLE TURBULENT CAVITATING FLOWS 1)
Wang Changchang, Wang Guoyu, Huang Biao, Zhang Mindi
2019, 51(5):  1296-1309.  DOI: 10.6052/0459-1879-19-128
Asbtract ( 430 )   HTML ( 17)   PDF (28608KB) ( 201 )
Figures and Tables | References | Related Articles | Metrics

Cavitation is the strongly compressible flows, where the cavitation compressibility effects could be closely associated with the cavity instabilities and dynamics. To investigate the cavity structure evolution and dynamics under both the re-entrant jet and shock wave mechanisms in compressible cloud cavitating flows, experiments were conducted in the divergent section with 10° in a venturi, using the simultaneous sampling technique to synchronize the observations of transient cavity behaviors with the wall-pressure signals measurements in cloud cavitating flows. A novel image processing algorithm is developed to analyze the transient cavity structures in detail. Based on the developed image processing method, the attached-type sheet cavity and shedding cloud-like cavity structures can be studied quantitatively. Results showed that unsteady cavity behaviors under re-entrant jet mechanism (RJM) can be divided into three stages: 1) the attached cavity growth, 2) the re-entrant jet formation and movement, and 3) the attached cavity breakup and cavity cloud shedding, collapse. The transient cavity structures under shock wave mechanism (SWM) present periodic behaviors: 1) the attached cavity growing, 2) the shock wave generating and propagating, and 3) the attached cavity collapsing and shedding. The period of re-entrant jet development is larger than that during the shock wave propagation. The interactions between shock wave and cavity cause significant local void fraction reduction, which induces complex cavity dynamics. During the shock wave propagation, both large pressure fluctuations and pressure peaks are captured, while in the process of the re-entrant jet movement, relatively smooth pressure fluctuations are measured, with small pressure fluctuations increase at the re-entrant jet head. Furthermore, the behaviors of attached-type sheet cavity, and shedding cloud-like cavity are totally different under re-entrant jet and shock wave mechanisms,indicating different mass transfer processes in different cavity shedding mechanisms.

INFLUENCE OF THE LENGTH OF HIGH-SPEED TRAIN ON THE FAR-FIELD AEROACOUSTICS CHARACTERISTICS 1)
Mo Huangrui, An Yi, Liu Qingquan
2019, 51(5):  1310-1320.  DOI: 10.6052/0459-1879-19-079
Asbtract ( 453 )   HTML ( 32)   PDF (11317KB) ( 222 )
Figures and Tables | References | Related Articles | Metrics

The estimation of the aeroacoustic characteristics of the high-speed train is of great importance in train design. Due to the extreme slender shape of the full marshaling of the train (generally 8 or 16 cars), the calculation of the far filed acoustics involves massive consumption of the CPU time which even makes the optimization unpractical. This paper numerically studies the influence of train length on its aeroacoustics characteristics by considering the acoustic contribution car by car. The nonlinear acoustic solver (NLAS) is used to calculate the acoustic source at the acoustic surface, and the FW-H analogy method is used to integrate on the acoustic surface to obtain the far-field results. Three different marshaling of the trains, i.e. 3 cars, 4 cars, and 6 cars, are studied and compared. The far-field acoustic level and frequency profile along the train is obtained. The results show that for different middle cars, the far-field acoustic profile along the length is very similar in both quantity and shape, except the offset in position. For the cars in the same position of different marshaling, both the sound pressure level distribution and the frequency profile are also very close to each other. Thus the key aeroacoustics characteristics of a long marshaling train could be estimated with a much smaller marshaling such as 3 cars. With the superposition of the acoustic surface data from short marshaling simulation, the aeroacoustics characteristics of long marshaling could be obtained. The comparison of the superposed results and the results calculated directly from long marshaling simulation is close enough for engineering use. This demonstrates that the proposed novel approach for estimating aeroacoustics of long marshaling could not only reduce the computational cost significantly but also be as accurate as direct simulation. This paper might provide a handy tool for engineering practice in this region.

EFFECT OF WAKE INDUCED-VIBRATION RESPONSES OF A SQUARE CYLINDER BEHIND THE STATIONARY SQUARE CYLINDER 1)
Tu Jiahuang, Tan Xiaoling, Yang Zhilong, Deng Xuhui, Guo Xiaogang, Zhang Ping
2019, 51(5):  1321-1335.  DOI: 10.6052/0459-1879-19-051
Asbtract ( 282 )   HTML ( 15)   PDF (25866KB) ( 128 )
Figures and Tables | References | Related Articles | Metrics

Wake-induced vibration (WIV) of a two-degree-of-freedom downstream square cylinder behind a stationary upstream square one are numerically investigated at low Reynolds numbers by using semi-implicit Characteristics-based split(CBS) finite element algorithm in this study. Firstly, the results are compared with the existing data in the literature to verify the accurateness of the method. Then, the influence of Reynolds number ($Re$) and reduced velocity ($U_{\rm r})$on the wake-induced vibration response of the downstream square cylinder is analyzed. Meanwhile, the results are also compared with those of the single square cylinder case. The numerical results show that the $Re$and $U_{\rm r}$have great influence on the dynamic response characteristics, such as amplitude, vibration frequency and motion trajectory of the downstream square cylinder. With the increasing of $Re$, the interaction between both cylinders changes from the vortex-induced interference to the wake one, resulting in the change of the frequency characteristics and the intensifying dynamic response of the downstream circular cylinder. The trajectories resemble figure "8" in the single square cylinder case. However, in the case of both cylinders arranged in-line, the trajectories of the downstream cylinder become complicated with the increasing of $Re$. The motion orbits are basically "8" figure at $Re=40$, 80. When the $Re$are 120, 160 and 200, the "dual-8" figure can be observed. Meanwhile, the wake field characteristics of the downstream square cylinder show 2S, 2S*, 2P, 2T, P+S and steady mode. Finally, the wake-induced vibration mechanisms of two square cylinders in the tandem arrangement is revealed according to the analysis of the characteristic of the flow field.

EXPERIMENT STUDY OF EFFECT OF NONEQUILIBRIUM PLASMA ON METHANE-OXYGEN DIFFUSIVE FLAME 1)
Zhou Siyin, Nie Wansheng, Che Xueke, Tong Yiheng, Zheng Tikai
2019, 51(5):  1336-1349.  DOI: 10.6052/0459-1879-19-149
Asbtract ( 305 )   HTML ( 6)   PDF (30526KB) ( 75 )
Figures and Tables | References | Related Articles | Metrics

Based on the self-designed plasma injector, a nonequilibrium oxygen plasma is generated by dielectric barrier discharge to study the effect of pure oxygen plasma on methane-oxygen diffusive flame. The discharge characteristics, thermal and gas dynamic effects of the plasma, together with the flame characteristic parameters, the discharge power, and the cost-benefit ratio are all experimentally analyzed by utilizing various diagnostic methods, such as schlieren imaging, thermocouple, infrared thermometer, broadband and CH* chemiluminescence imaging. Results show that the discharge current increases notably due to combustion. The current increases with the discharge voltage rises, yet it decreases with the flow rate of oxygen increases. The discharge power of oxygen plasma is higher than that of air plasma, yet the light emission intensity is obviously weaker under a certain flow rate and voltage. The heating of discharge plasma is mainly restricted within the injector. Since the maximum temperature increment led by the plasma is only 30.6 K, it is assumed that the thermal effect of oxygen plasma is too weak to influence the flame combustion. Moreover, according to the experimental results of air discharge plasma, the heat release process of the discharge reactions should mainly be determined by species, which contain oxygen. An induced jet, which has three velocity components, is generated by the discharge. The angle of the original oxygen jet is enlarged due to this induced jet, especially for a higher voltage. The relative capability of the plasma on enlarging the jet angle is stronger under a lower flow rate. The flame becomes more stable and its heat release enhances under certain conditions mainly due to the plasma gas dynamic effect, which changes the mixing between fuel and oxidizer. The lowest cost-benefit ratio of the plasma is only 2.2％ in this study, and the plasma performs better under a high flow rate or a small mixing ratio.

THE NUMERICAL INVESTIFATION ON HYDRODYNAMIC AND STRUCTURAL STRENGTH OF A COMPOSITE HYDROFOIL 1)
Chen Qian, Zhang Hanzhe, Wu Qin, Fu Xiaoying, Zhang Jing, Wang Guoyu
2019, 51(5):  1350-1362.  DOI: 10.6052/0459-1879-19-107
Asbtract ( 263 )   HTML ( 6)   PDF (24814KB) ( 100 )
Figures and Tables | References | Related Articles | Metrics

This paper investigates systematically the structural deformation characteristics of composite hydrofoils with own unique properties by numerical method to solve the fluid structure interaction problem on the composite hydrofoils. It establishes a numerical model of composite hydrofoils with fluid structure interaction in this study, which verifies the correctness of the numerical model by comparing the calculated results to the experimental results of Zarruk, and the numerical results show that the tip torsion angle of the composite hydrofoils increases with the increase of the Reynolds number. This paper investigates and analyses systematically the influence of the ply angle on the hydrodynamic performance and strength performance of composite hydrofoils in view of the numerical method. The results illustrate that the tip torsion angle of composite hydrofoils decreases with increasing the ply angle, while the tip displacement of composite hydrofoils decrease firstly and then increase with increase of ply angle. The study proposes two dimensionless, the dimensionless torsion angle and the dimensionless displacement, to eliminate the effect of engineering constants of composite on the structure deformation of composite hydrofoils, which is used to further investigate the influence of the bending-torsion coupling effect on the deformation characteristics of the composite hydrofoils. The strength of the composite hydrofoils is judged and studied by using the Tsai-Wu failure criterion. The numerical results show that the Tsai-Wu coefficient of composite hydrofoils decreases firstly and then increases with the increase of the ply angle, where the lowest Tsai-Wu coefficient and the largest Tsai-Wu coefficient appears in the 0$^\circ$layered composite hydrofoil and the 50$^\circ$layered composite hydrofoil, respectively.

Soild Mechanics
PROBABILISTIC CONTROL VOLUME METHOD FOR THE SIZE EFFECT OF SPECIMEN FATIGUE PERFORMANCE1)
Li Yabo, Song Qingyuan, Yang kai, Chen Yiping, Sun Chengqi, Hong Youshi
2019, 51(5):  1363-1371.  DOI: 10.6052/0459-1879-19-118
Asbtract ( 393 )   HTML ( 4)   PDF (3758KB) ( 226 )
Figures and Tables | References | Related Articles | Metrics

CONTINUUM DAMAGE MECHANICS-BASED CRITICAL SINGULARITY EXPONENT AND FAILURE TIME PREDICTION1)
Zhou Sunji, Cheng Lei, Wang Liwei, Wang Ding, Hao Shengwang
2019, 51(5):  1372-1380.  DOI: 10.6052/0459-1879-19-120
Asbtract ( 366 )   HTML ( 8)   PDF (4005KB) ( 166 )
Figures and Tables | References | Related Articles | Metrics

The accelerating increase of response quantities, such as strain and acoustic emission signals in the vicinity of failure time, has been revealed as a precursor of catastrophic failure. This critical precursor has been widely validated as a valid way to predict failure time by the retrospective prediction of volcanic eruptions, landslides and laboratory experiments of rock failures. But the scatter of exponents in the critical power law singularity relationship that describes acceleration in precursory signals leads to a debate on the actual value of critical exponent and the uncertainty of failure time prediction. Consequently, the uncertainty resulting from the scatter of exponents is a key difficulty in using such methods for prediction of the failure time through the use of acceleration precursors. Thus understanding the underlying mechanisms for the magnitude and variation of critical power-law exponents becomes a central problem for understanding the process of failure and failure time prediction. This paper presents a multi-scale damage mechanic model describing the accelerating process of time dependent failure. Theoretic derivations and demonstrations of critical power law relationship are presented for two typical load process, i.e. brittle creep failure under constant nominal stress and the load process of linearly increasing the nominal stress with time. It is found that values of the singularity exponents have a relationship with the parameter that defines the nonlinear level of damage evolution rate depending on to the local true stress. The physical expressions of critical parameters are deduced and the physical meanings of these critical parameters are explained. It is declared that the observed variation of the critical power law exponents not only due to the fluctuation in the measurement data, but has its intrinsic physical controls. Then a method is suggested to predict the failure time when the critical singularity exponent is unknown. This proposed methodology is validated through granite creep failure experiments in laboratory and the challenges for practical applications are demonstrated.

COUPLING EFFECTS OF WRINKLES AND GRAIN BOUNDARY ON THE FRACTURE OF GRAPHENE1)
Ren Yunpeng, Cao Guoxin
2019, 51(5):  1381-1392.  DOI: 10.6052/0459-1879-19-181
Asbtract ( 288 )   HTML ( 4)   PDF (15708KB) ( 61 )
Figures and Tables | References | Related Articles | Metrics

Graphene fabricated via chemical vapor deposition (CVD) is typically polycrystalline and also includes many wrinkles. The fracture of the polycrystalline graphene with wrinkles under an uniaxial tensile load is investigated via molecular dynamics simulations. With a tensile load perpendicular to the grain boundary, wrinkles can significantly increase the failure stress of bi-crystalline graphene with a small tilt angle, and the increase in the failure stress is up to around 50%. The wrinkle effect on the failure stress decreases with the increase of the tilt angle, resulting in that the failure stress of bi-crystalline graphene is insensitive to the tilt angle and slightly lower than that of pristine graphene, which agrees with the experimental results very well. With a tensile load along the grain boundary, the failure stress is insensitive to the wrinkle. In addition, wrinkles can significantly increase the failure strain, up to 100%. The influence mechanism can be described as follows: wrinkles will cause the out-of-plane deformation in graphene, resulting a partial release of the tensile pre-stress induced via the 5-7 rings of the grain boundary and consequently an increase in the failure stress of bi-crystalline graphene; the interaction between 5-7 rings of the grain boundary is eliminated, resulting that the failure stress is insensitive to the tilt angle; the flattening of wrinkles can significantly decrease the stretching ratio of C---C bonds, resulting an obvious increase in the failure strain. The present study provides a useful help to understand the fracture of graphene.

STUDY ON TENSILE FRACTURE BEHAVIOR AND MECHANICAL PROPERTIES OF GO BASED ON MOLECULAR DYNAMICS METHOD1)
Li Dongbo, Liu Qinlong, Zhang Hongchi, Lei Pengbo, Zhao Dong
2019, 51(5):  1393-1402.  DOI: 10.6052/0459-1879-19-175
Asbtract ( 591 )   HTML ( 15)   PDF (29211KB) ( 402 )
Figures and Tables | References | Related Articles | Metrics

In comparison with graphene, the graphene oxide (GO) has better hydrophilicity, dispersion performance and reaction activity, which makes it easier to interact with other materials to form composites with excellent properties. But on the other hand, due to the complexity of the electronic structure of the GO, there are some differences in the research of mechanics of GO at present. For this purpose, in the present paper, a random distribution model of GO containing the hydroxyl groups,epoxy groups and carboxyl groups was established based on the molecular dynamics method. And then, the tensile fracture behavior was analyzed by uniaxial tensile numerical simulation. The results illustrated that the epoxy groups which are far away from hydroxyl and carboxyl groups had induce effect on the fracture of GO. In addition, the mechanisms were explained from three aspects,which are the chemical bonding, system energy and stress distribution. Additionally, the influence of the coverage of oxygen-containing functional groups of hydroxyl groups, epoxy groups and carboxyl groups on the stress-strain curve, ultimate strength and ultimate strain of GO were further studied. The results showed that the ultimate strength and ultimate strain of GO decreased with the increase of the coverage of oxygen-containing functional groups of hydroxyl groups,epoxy groups and carboxyl groups. According to the analysis, the main reason is that the oxygen-containing functional groups of hydroxyl groups, epoxy groups and carboxyl groups destroy the sp$^{2}$hybrid form in the original graphene surface, thus weakening the bonding energy between atoms. Thereby, The greater the coverage of the oxygen-containing functional groups of hydroxyl groups, epoxy groups and carboxyl groups, the greater the amount and degree of weakened bonding energy, and the lower the ultimate strength of GO. The research results have certain reference value and significance for engineering application and basic research of GO.

THEORETICAL PREDICTION AND EXPERIMENTAL VERIFICATION OF WRINKLE AMPLITUDE IN A SQUARE MEMBRANE SUBJECTED TO DIAGONAL TENSION1)
Cao Jinjun, Zhang Huiting, Zhang Liang, Peng Fujun, Yun Weidong
2019, 51(5):  1403-1410.  DOI: 10.6052/0459-1879-19-109
Asbtract ( 285 )   HTML ( 6)   PDF (5970KB) ( 106 )
Figures and Tables | References | Related Articles | Metrics

Flexible membrane structures are widely used in the key parts of aerospace vehicles. Surface flatness is one of the main factors affecting the performances of membrane structures. Wrinkle amplitude is an important factor for evaluating the surface flatness of membrane reflector antennas. Wrinkle amplitude is strongly related with the transverse strain, perpendicular to wrinkling direction. Based on the stability theory of thin plates, a theoretical model is proposed to predict the wrinkle amplitude in a square membrane subjected to diagonal tension. The effect of transverse tensile force on the membrane deformation is taken into account. The displacement perpendicular to wrinkling direction is decomposed into three parts: the transverse displacement induced by Poisson's effect, the wrinkling displacement induced by out-of-plane deformation, and the tensile displacement induced by transverse tensile force. The formulation of wrinkle amplitude is reworked. Based on digital image correlation (DIC) technology, speckle experiment is carried out for a square membrane subjected to diagonal tension. The three-dimensional displacement of square membrane is measured by the binocular vision three-dimensional measurement system. The three-dimensional deformed shapes and wrinkle waveforms of membrane are obtained. We study the nonlinear relationship between wrinkle amplitudes and tensile loads. Compared with an existing model, our model greatly improves the accuracy of prediction to wrinkle amplitude, which is in good agreement with experimental results. The theoretical research presented in this paper can provide valuable guidance for the establishment of a fine numerical model and the implementation of algorithm.

MULTI-SCALE METHOD OF PLAIN WOVEN COMPOSITES SUBJECTED TO LOW VELOCITY IMPACT BASED ON ASYMPTOTIC HOMOGENIZATION1)
Zhang Jiehao, Duan Yuechen, Hou Yuliang, Tie Ying, Li Cheng
2019, 51(5):  1411-1423.  DOI: 10.6052/0459-1879-19-133
Asbtract ( 326 )   HTML ( 14)   PDF (40088KB) ( 142 )
Figures and Tables | References | Related Articles | Metrics

A multi-scale approach was presented to analyze low velocity impact response and damage of plain woven composites. Firstly, by using the maximum principal stress failure criterion and direct stiffness degradation model to characterize the damage initiation and damage evolution of fiber and matrix, micro-scale unit cell under the periodical boundary condition was established to predict the elastic and strength properties of fiber bundles, which were substituted into the meso-scale unit cell. After that, the progressive damage simulation of meso-scale unit cell under six boundary conditions was carried out based on the mixed failure criteria of Hashin and Hou, and continuum damage model. Then the effective properties of 0$^\circ$and 90$^\circ$subcell were predicted based on the asymptotic homogenization method by using meso-scale unit cell as the media, and the subcell model of plain woven composites was established. The subcell model was then extended into a macro-scale low velocity impact model. Based on the above methods, the mechanical response and damage characteristics of plain woven composites under low velocity impact were studied. The results show that macro-scale impact simulation results agree well with experimental results, which verifies the correctness of multi-scale approach. The maximum contact force, absorbed energy and delamination area increase with the increasing impact energy, and the delamination damage morphology gradually transforms from ellipse to circle.The long axis direction of matrix tensile damage and matrix compressive damage are orthogonal and consistent with the material principal direction of subcell respectively, and the damage area of the former is much larger than that of the latter.

HYBRID PERTURBATION-GALERKIN METHOD FOR GEOMETRICAL NONLINEAR ANALYSIS OF TRUSS STRUCTURES WITH RANDOM PARAMETERS1)
Huang Bin, He Zhiyun, Zhang Heng
2019, 51(5):  1424-1436.  DOI: 10.6052/0459-1879-19-099
Asbtract ( 304 )   HTML ( 7)   PDF (8701KB) ( 97 )
Figures and Tables | References | Related Articles | Metrics

The hybrid perturbation-Galerkin stochastic finite element method is used to solve the geometrical nonlinear truss structures with random parameters. The power polynomial expansions are adopted to express the random secant elastic modulus and random responses with respect to displacement terms, respectively. Using the high-order perturbation method, the coefficients of the power polynomial expansions can be obtained, so that the expression of the geometrical nonlinear displacement can be determined. The coefficients terms of the power polynomial expansions obtained are used as the Galerkin trial functions, and the Galerkin projection technique is employed to determine the coefficients of these trial functions. Since the trial functions come from the linear combination of the perturbation solutions, the trial functions selected are self-adaptive to the nonlinear problem. The numerical example about multi random variables with different probability distributions show that since no probability conversion is required for the proposed method so that the conversion errors in the calculation process are avoidable, which results in that the accuracy of the proposed method is higher than that of generalized polynomial chaos method (GPC method). Meanwhile, when the results are equally accurate, the nonlinear algebraic equations about the coefficients of the trial functions obtained by the proposed method is more easy to be solved than that by the GPC method, and the calculation cost of the suggested method is less than that of the GPC method. When the fluctuation of the random variable becomes large, the statistical moments of the structural response calculated by the hybrid perturbation-Galerkin method are closer to the results of the Monte Carlo simulation method than that of the high-order perturbation method, which illustrates that the hybrid perturbation-Galerkin method is effective in solving the stochastic geometric nonlinear problems.

Dynamics, Vibration and Control
A NEW METHOD FOR THE PROBABILITY DENSITY OF MAXIMUM ABSOLUTE VALUE OF A MARKOV PROCESS1)
Chen Jianbing, Lü Mengze
2019, 51(5):  1437-1447.  DOI: 10.6052/0459-1879-19-104
Asbtract ( 257 )   HTML ( 8)   PDF (11863KB) ( 76 )
Figures and Tables | References | Related Articles | Metrics

The probability distribution of maximum absolute value of stochastic processes or responses of stochastic systems is a significant problem in various science and engineering fields. In the present paper, theoretical and numerical studies on the time-variant extreme value process and its probability distribution of Markov process are performed. An augmented vector process is constructed by combining the maximum absolute value process and its underlying Markov process. Thereby, the non-Markov maximum absolute value process is converted to a Markov vector process. The transitional probability density function of the augmented vector process is derived based on the relationship between the maximum absolute process and the underlying state process. Further, by incorporating the Chapman-Kolmogorov equation and the path integral solution method, the numerical methods for the probability density function of the maximum absolute value is proposed. Consequently, the time-variant probability density function of the maximum absolute of a Markov process can be obtained, and can be applied, e.g., to the evaluation of dynamic reliability of engineering structures. Numerical examples are illustrated, demonstrating the effectiveness of the proposed method. The proposed method is potentially to be extended for the estimation of extreme value distribution of more general stochastic systems.

STUDY ON THE RIGID-LIQUID COUPLING DYNAMICS OF A CYLINDRICAL-TANK SPACECRAFT IN TRANSLATION BASED ON MUTI-SCALE METHOD1)
Li Xiaoyu, Yue Baozeng
2019, 51(5):  1448-1454.  DOI: 10.6052/0459-1879-19-111
Asbtract ( 325 )   HTML ( 5)   PDF (1372KB) ( 84 )
Figures and Tables | References | Related Articles | Metrics

Based on the engineering background of liquid-filled spacecraft, the nonlinear dynamic characteristics of rigid-liquid coupling dynamic system are studied by means of multi-scale method. By using the multi-dimensional modal method, the free boundary problem describing the nonlinear sloshing of the liquid in a cylindrical tank under horizontal excitation is transformed into a finite dimensional nonlinear ordinary differential equation system in which the coupled liquid sloshing modal coefficients are assumed as state variables. The analytical expressions of sloshing force and sloshing moment acting on the tank wall caused by liquid sloshing are derived, and then the coupled dynamic equations for the spacecraft translation and liquid fuel sloshing are obtained. The dynamic characteristics, especially the first-order primary resonance, of rigid-liquid coupling system are analyzed by multi-scale method. The natural frequency of rigid-liquid coupling system is solved by the characteristic equation of natural frequency, and the amplitude-frequency response equation of liquid steady-state solution is derived when the external excitation frequency is close to the first natural frequency of coupling system. Combined with the numerical method, the amplitude-frequency response curve and excitation-amplitude response curve of liquid steady state solution are studied. The results show that with the change of liquid filling ratio, the amplitude frequency response curve of liquid steady state solution will exhibit soft, hard spring features conversion phenomenon. In addition to that, a `jump' phenomenon was also observed during this kind of soft-hard spring feature conversion. Furthermore, it is shown that the soft and hard spring characteristic conversion point of amplitude frequency response curve will be affected by gravity acceleration and spring stiffness coefficient. The above results show that the dynamic characteristics of the rigid-liquid coupling system considering the nonlinear effect are essentially different from those shown by the traditional linear system model. The investigations presented have important reference value for further analysis of rigid-liquid coupling nonlinear dynamic characteristics of liquid-filled spacecraft.

RESEARCH ON ELASTIC LINE METHOD BASED ON ABSOLUTE NODAL COORDINATE METHOD1)
Fan Jihua, Zhang Dingguo, Shen Hong
2019, 51(5):  1455-1465.  DOI: 10.6052/0459-1879-19-076
Asbtract ( 555 )   HTML ( 8)   PDF (1864KB) ( 143 )
Figures and Tables | References | Related Articles | Metrics

Compared with the floating frame of reference formulation, the absolute nodal coordinate formulation (ANCF) has significant advantages in dealing with the nonlinear large deformation problem of flexible bodies. ANCF defines the nodal co-ordinates in a global co-ordinate system, uses the global slopes instead of angles to define the orientation of the elements, has a constant mass matrix, and does not have the Coriolis centrifugal force. However,the elastic force matrix is a nonlinear term, and the solution will be time consuming and occupied resources. According to this,the elastic line method is introduced for solving the elastic force, the Green-Lagrangian strain tensor is defined on the centerline, the curvature formula is used to define the bending strain, and the torsion angle formula is used to define the torsional strain. At the same time, the finite element method is used to describe the displacement field of the three-dimensional flexible beam, and the constant mass matrix, the generalized stiffness matrix and the generalized force matrix of the beam element are solved, and then the dynamic equations of the element are obtained. The dynamic equations of the three-dimensional beam are obtained by the transformation matrix. Then the differences between the continuum mechanics method and the elastic line method are pointed out theoretically, and the dynamic simulation software is compiled. Finally, the dynamics behaviors of a flexible pendulum and a tracked vehicle are numerical investigated with the continuum mechanics method and the elastic line method. The results show that the elastic line method can effectively improve the calculation efficiency under the premise of ensuring accuracy.

PERFORMANCE ANALYSIS OF GROUNDED THREE-ELEMENT DYNAMIC VIBRATION ABSORBER1)
Xing Zikang, Shen Yongjun, Li Xianghong
2019, 51(5):  1466-1475.  DOI: 10.6052/0459-1879-19-154
Asbtract ( 351 )   HTML ( 4)   PDF (17559KB) ( 61 )
Figures and Tables | References | Related Articles | Metrics

In grounded dynamic vibration absorbers (DVA), the changing tendencies of the fixed-point amplitude with the natural frequency ratio are not monotonous. Thus, the results obtained by optimizing this type DVA based on classical fixed-point theory may be the local optimum parameters. The primary system can obtain smaller amplitudes when selecting other parameters. The optimization of grounded type DVAs are worthy to further study. In addition, the damper of the DVA inevitably has some elasticity. Accordingly, a grounded three-element DVA is studied by analyzing the influence of system parameters on fixed-point positions and maximum amplitude in this paper. The local optimum parameters of the DVA are obtained and the performance is investigated. Firstly, the motion differential equation of the system is established, and the amplitude amplification factor of the primary system is obtained. It is found that there are three fixed-points independent of damping on the amplitude-frequency response curve. In most cases, by optimizing the damping ratio, the tendency of the larger of the fixed point changing with the system parameters can represent the tendency of the maximum amplitude changing with the system parameters. Therefore, the expressions of the fixed-point are obtained by using the Shengjin's formula. For more accuracy, the numerical algorithm is used to obtain the relationship between the maximum amplitude and the system parameters, and it is found that there are local optimum parameters in the system. Finally, in order to obtain the local optimum parameters the grounded three-element DVA is compared with the grounded DVA. The study shows that the local optimum parameters of the two DVAs are the same except the stiffness ratio. When the natural frequency is smaller than local optimum frequency ratio, the maximum amplitude of the primary system of the grounded three-element DVA model is much smaller than that of the grounded DVA.

VIBRATION EXPERIMENT AND NONLINEAR MODELLING RESEARCH ON THE FOLDING FIN WITH FREEPLAY
He Haonan, Yu Kaiping, Tang Hong, Li Jinze, Zhou Qiankun, Zhang Xiaolei
2019, 51(5):  1476-1488.  DOI: 10.6052/0459-1879-19-119
Asbtract ( 434 )   HTML ( 6)   PDF (19220KB) ( 152 )
Figures and Tables | References | Related Articles | Metrics

The influence of the freeplay on ground vibration responses of the folding fin and the nonlinear modelling of the freeplay, which exists in the joints connecting the inboard fin and the outboard fin, are studied in this paper. The first five modal parameters are obtained by hammer tests of the linear structure without freeplay. The sweep base excitation is then applied to the root-fixed folding fin with adjustable freeplay. The experimental results show that the existence of freeplay can lead to some nonlinear phenomena of the structure, such as a forward and backward sweep difference, a jump, multi harmonics, and a frequency shift. This kind of lumped nonlinearity has a huge influence on the first bending mode. For example, both the increase of the excitation magnitude and the decrease of the angle of freeplay will raise the fundamental frequency which gradually tends to the result of the linear structure without freeplay, but they have little effects on the second torsion mode compared with those on the first bending mode. The finite element model of the folding fin is established in the numerical simulation analysis and it is reduced by the developed model reduction method which is suitable for similar folding structures with lumped nonlinearity. The freeplay and contact of the joints are simulated by nonlinear torsional springs with the combination of linear and 3/2 order stiffness according to the Hertz theory. The linear part of the folding fin is verified firstly via comparing the first four frequencies and mode shapes computed from numerical analysis and the above hammer test. The dynamic response of the nonlinear structure under base excitation is calculated by Bathe's two-step implicit composite algorithm. The transfer functions obtained can simulate the frequency variation characteristics in the experiments, thus verifying the proposed nonlinear modelling method of the joint.

RESEARCH ON MULTI-STABILITY CHARACTERISTICS OF GEAR TRANSMISSION SYSTEM WITH TWO-SPACE COUPLING1)
Shi Jianfei, Gou Xiangfeng, Zhu Lingyun
2019, 51(5):  1489-1499.  DOI: 10.6052/0459-1879-19-093
Asbtract ( 278 )   HTML ( 8)   PDF (23437KB) ( 62 )
Figures and Tables | References | Related Articles | Metrics

By defining the system parameters as parameter variables and forming the parameter space, the nonlinear global dynamics of the gear transmission system under the coupling of parameter space and state space are studied in detail in this work. The correlative relationship between multiple parameters, multiple initial values and multiple stable behaviors is also obtained. Firstly, a method for calculating and identifying the multi-stable behavior of a nonlinear system under the coupling of two spaces is designed. Secondly, based on the designed method and combined with phase diagram, Poincaré map, bifurcation diagram, top Lyapunov exponent and basin of attraction, the existence and distribution of multi-stable behavior for the gear transmission system in different parameter planes are investigated numerically to better understand the motion mechanism of the system. In addition, the distribution characteristic of multi-stable behavior in the state plane is also studied on the base of the cell-to-cell mapping method. The multi-stable behavior and bifurcation that may be hidden in the parametric plane and the state plane are fully revealed. The formation mechanism of multi-stable behavior is analyzed as well. The results show that there are a large number of multiple stable behaviors which are banded distribution in the parametric planes of the gear system under the coupling of two spaces. Two different erosion phenomena, such as internal erosion and boundary erosion, are clearly observed in the state plane. The sensitivity of bifurcation points or bifurcation curves to initial values leads to the occurrence of multi-stable behavior. When the amplitude of backlash or error fluctuation changes in a small range, the global dynamic characteristics of the gear system are less affected by the backlash or error disturbance. However, the global dynamic characteristics are greatly affected by the meshing frequency. Global dynamic characteristics of the gear system become rich and complex under two-space coupling.

Biomechanics,Engineering and Interdiscipliary Mechanics
A RESONANCE FREQUENCY EXTRACTION METHOD FROM LOW Q-FACTOR MATERIALS BASED ON RESONANT ULTRASOUND SPECTROSCOPY1)
Zhang Qiang, Fan Fan, Wang FeiNiu, Shen Fei, Niu Haijun
2019, 51(5):  1500-1506.  DOI: 10.6052/0459-1879-19-049
Asbtract ( 326 )   HTML ( 2)   PDF (3928KB) ( 69 )
Figures and Tables | References | Related Articles | Metrics

Resonant ultrasound spectroscopy (RUS) allows identification of the elastic coefficients of solid materials vibrating under an ultrasonic excitation from the measurement of their inherent frequencies. Retrieving the resonant frequencies is therefore a key signal processing step in RUS. However, according to the attenuation characteristics of low $Q$-factor (quality factor) materials, the resonance spectrum obtained by the experiment is flat and the resonance frequencies can not be directly observed from the spectrum. Therefore, in order to retrieve more effective resonance frequencies than traditional approach from the low $Q$-factor materials, a new extraction method of resonance frequencies was proposed to solve the limitation in this paper. The empirical mode decomposition method was used to decompose the frequency response of the specimen into finite Intrinsic Mode Function (IMF) components with special oscillation characteristics. According to the prior information of resonant ultrasound spectroscopy (RUS),the relevant IMF component was selected to retrieve reliable resonance frequencies from the resonance spectrum. The short fiber filled epoxy (a kind of bone-like materials, $Q \approx$25) was adopted as the specimen to calculate the elastic coefficients and engineering moduli compared with the traditional linear prediction method. The experimental results show that the new method has high computational efficiency and is more sensitive to the weak excitation modes of low $Q$-factor materials. The number of effective resonance frequencies (26) are more than traditional linear prediction methods (21), which also satisfied 5 times estimation requirement of elastic constants. In addition, the optimized elastic moduli are closer to the standard values of the short fiber filled epoxy. In conclusion, the EMD-based method can retrieve a sufficient quantity and effective resonance frequencies from the flat spectrum of low $Q$-factor materials, which can not only improve the reliability of the estimation of mechanical parameters, but also extend the application range of resonant ultrasound spectroscopy.

AN INTEGRAL DIFFERENTIATION PROCEDURE FOR DYNAMIC TIME-HISTORY RESPONSE ANALYSIS OF STRUCTURES1)
Li Hongjing, Mei Yuchen, Ren Yongliang
2019, 51(5):  1507-1516.  DOI: 10.6052/0459-1879-19-105
Asbtract ( 295 )   HTML ( 2)   PDF (963KB) ( 137 )
Figures and Tables | References | Related Articles | Metrics

Traditionally, the response of displacement is selected to be the basic unknown and the responses of velocity and acceleration are usually expressed by linear weighted sum of the displacement when the differential quadrature (DQ) method is applied to the solution of dynamic problems. Either the linear equations or matrix equations (Sylvester equation) has to be processed in such procedure for dynamic solutions and the derived algorithm is conditionally stable in general. In this paper,the DQ principle is used in the inverse way to implement a high-accuracy explicit algorithm for the operation of convolution, and the algorithm is applied to dynamic analysis via the solution of Duhamel's integral. The dynamic response can be solved over a finite time interval according to this procedure, so that the total time-history of response could be obtained step by step. The inverse of Vandermonde matrix is required only once if the distribution of DQ nodes are completely consistent in each time interval and the response at several time instants during the interval can be obtained simultaneously. Hence, the procedure for dynamic solutions numerically achieves a high computational efficiency. It is proved that the spectral radius of the transfer matrix in the dynamic algorithm is always equal to 1, so the algorithm has unconditional stability and no numerical dissipation occurs during the calculation. The numerical accuracy of this algorithm depends on $N$, the number of DQ nodes within the analyzing time interval, and an algebraic accuracy with the order of $N-1$can be achieved. In practice, 10 and even more DQ nodes of $N$are suggested in order to gain high accuracy for dynamic problems.

A UNIFIED COMPUTATIONAL FRAMEWORK FOR FLUID-SOLID COUPLING IN MARINE EARTHQUAKE ENGINEERING: IRREGULAR INTERFACE CASE1)
Chen Shaolin, Cheng Shulin, Ke Xiaofei
2019, 51(5):  1517-1529.  DOI: 10.6052/0459-1879-19-156
Asbtract ( 319 )   HTML ( 5)   PDF (8787KB) ( 84 )
Figures and Tables | References | Related Articles | Metrics

The simulation of seismic wavefield at seafloor and ocean acoustic field, the influence of complex seabed media and seabed topography needs to be considered, involving the coupling between seawater, saturated seabed, elastic bedrock and structure. That means, we target simulation where several types of equations are involved such as fluid, solid and saturated porous media equation. The conventional method for this fluid-solid-saturated porous media interaction problem is to use exsisting solvers of different equations and coupling method, which needs data mapping, communication and coupling algorithm between different solvers. Here, we present an alternative method, in which the coulping between different solvers is avoided. In fact, when porosity equals to one and zero,the saturated porous media is reduced to fluid and solid respectively, so we can use the porous media equation to describe the ideal fluid and solid, and the coupling between porous media,solid and fluid turns to the coupling between porous media with different porosity. Based on this idea, firstly the Biot's equations are approximated by Galerkin scheme and the explicit lumped-mass FEM is chosen, that are well suited to parallel computation. Then considering the conditions of coupling on the irregular interface between porous media with different porosity,by solving the normal and tangential interface forces, the coupled algorithm is derived, which is proved to be suitable for the coupling between fluid,solid and saturated porous media. Thus, the coupling problem between fluid, solid and saturated porous media can be brought into a unified framework, in which only the solver of saturated porous media is used. The three-dimensional parallel code for this proposed method is programed. Considering the situation of sag terrain in water-bedrock, water-saturated seabed-bedrock system, the unified framework proposed in this paper is combined with the transmission boundary conditions to analyze the dynamic response of P wave incident, and the unified framework are verified by the results meeting the interface conditions.

INFLUENCE OF VORTEX CORE SIZE ON AERODYNAMIC CALCULATION OF WIND TURBINE IN FREE VORTEX WAKE METHOD1)
Xu Bofeng, Liu Bingbing, Feng Junheng, Zuo Lu
2019, 51(5):  1530-1537.  DOI: 10.6052/0459-1879-18-440
Asbtract ( 469 )   HTML ( 3)   PDF (4204KB) ( 103 )
Figures and Tables | References | Related Articles | Metrics

A MODEL OF NEUTRON IRRADIATION EMBRITTLEMENT FOR METALS1)
Ye Xiangping, Liu Cangli, Cai Lingcang, Hu Changming, Yu Yuying, Hu Ling
2019, 51(5):  1538-1544.  DOI: 10.6052/0459-1879-19-025
Asbtract ( 310 )   HTML ( 4)   PDF (16534KB) ( 73 )
Figures and Tables | References | Related Articles | Metrics

Irradiation embrittlement of metals is very important in the field of nuclear energy safety. In order to describe the irradiation embrittlement behavior of metals, a neutron irradiation embrittlement model for annealed metals was proposed based on Johnson-Cook constitutive model. The fracture true stress of the irradiated samples was taken as the same as the unirradiated sample. This model can predict the whole true stress-strain curve, and the fracture true strain of the irradiated annealed materials by using the yield strength only. The tensile true stress-strain curves, the fracture true stress, and the fracture true strain of high-purity aluminum with different doses were measured by quasi-static tensile tests. The results showed that a higher dose results in a higher yield strength and a lower fracture true strain. However, the fracture true stress is almost unchangeable. The size and number density of irradiation-induced defects for high-purity aluminum with different doses were mearsured by TEM microscope. The results showed that a higher dose results in a higher size and a higher number density of voids, but the size and the number density of dislocation loops were difficult to measure accurately due to their size and number density are too small. The parameters of this model were fitted by experimental data of high-purity aluminum, and the application effect of this model was checked. The predicted results of this model agree well with the experimental data, regardless wtether the yield was obtained by quasi-static tensile tests or by the size and number density of irradiation-induced defects. The critical irradiation dose of annealed high-purity aluminum predicted by this model also agrees well with literature.

AN EXTENDED STRENGTH AND YIELD CRITERION FOR GEOMATERIALS1)
Wan Zheng, Meng Da, Song Chenchen
2019, 51(5):  1545-1556.  DOI: 10.6052/0459-1879-19-074
Asbtract ( 264 )   HTML ( 5)   PDF (10832KB) ( 173 )
Figures and Tables | References | Related Articles | Metrics

Soil is a kind of typical friction material, while natural rock have certain cohesive properties and metal is a kind of completely cohesive material. Based on the analysis of three typical material strength criterion expressions, namely SMP, Lade-Duncan and Extended Von-Mises criterion, an extended criterion, VML criterion, is proposed by using the invariant expression of stress tensor, which can be degraded to the above three typical criteria respectively. In the deviatoric plane, the new criterion can be adopted to describe a variety of opening forms from curved triangle to circle, and in the meridian plane, the power function is adopted as the failure criterion formula, which can be used to describe the nonlinear property of the influence of hydrostatic pressure on the strength characteristics. As for the yield behavior of soil, geomaterials have typical compression-shear coupling characteristics. Therefore, in order to describe the volumetric coupling phenomenon under the two paths of shear and isotropic compression conditions, the droplet yield surface is adopted as the yield criterion. For the section shape on the deviatoric plane, the distribution and characteristics of the strength curve in the deviatoric plane with a fixed hydrostatic stress are discussed, and the influence of the stress Lode angle on the concavity and convexity of the strength curve in the deviatoric plane is discussed. Finally, the new criterion is used to verify the failure and yield test results of various materials. The comparison between prediction and test results have been employed. The comparison results have demonstrated the rationality of the proposed VML criterion.

Journal Comprehensive Information
ANALYSIS OF THE OUTPUT AND ACADEMIC INFLUENCE OF THE CHINESE JOURNAL OF THEORETICAL AND APPLIED MECHANICS1)
Liu Junli,Guan Cuizhong,Hua Fang
2019, 51(5):  1557-1564.  DOI: 10.6052/0459-1879-19-252
Asbtract ( 526 )   HTML ( 22)   PDF (11328KB) ( 126 )
Figures and Tables | References | Related Articles | Metrics

Under the leadership of Tsien Hsue-shen, the Chinese Journal of Theoretical and Applied Mechanics (CJTAM), which is the oldest comprehensive journal of mechanics founded in 1957, has witnessed the development track of Chinese mechanics discipline and the growth of several generations of Chinese mechanics talents since the 1950s. This paper utilized the bibliometric indicators to analyse the output and discipline influence of CJTAM from its inception to the present, and used the data to illustrate its contribution to the development of Chinese mechanics and the exchange of academic achievements. Firstly, we sorted out the amount of papers, the institutes and the authors published in CJTAM; Afterwards, its ranking among the journals of mechanics was analysed; Finally, we analysed the characteristics of the content of the paper, including the distribution of major disciplines and the research hotspots of recent publications.