Applied Mathematics and Mechanics was founded in 1980 by CHIEN Wei-zang, a celebrated Chinese scientist in mechanics and mathematics. The current editor in chief is Professor LU Tianjian from Nanjing University of Aeronautics and Astronautics. The Journal was a quarterly in the beginning, a bimonthly the next year, and then a monthly ever since 1985. It carries original research papers on mechanics, mathematical methods in mechanics and interdisciplinary mechanics based on artificial intelligence mathematics. It also strengthens attention to mechanical issues in interdisciplinary fields such as mechanics and information networks, system control, life sciences, ecological sciences, new energy, and new materials, making due contributions to promoting the development of new productive forces.
More>Articles in press have been peer-reviewed and accepted, which are not yet assigned to volumes /issues, but are citable by Digital Object Identifier (DOI).
Display Method:
2026, Volume 47, Issue 8
publish date：August 01 2026
Display Method:
2026, 47(8): 959-977.
doi: 10.21656/1000-0887.472064
Abstract:
Previous editorials of this journal have discussed academic lineage, intelligent tools, criteria for evaluating research, disciplinary direction, and human development in an age of increasingly powerful tools. A more prior question must now be asked: what gives a study the right to begin, and under what conditions does its underlying question genuinely stand?
Previous editorials of this journal have discussed academic lineage, intelligent tools, criteria for evaluating research, disciplinary direction, and human development in an age of increasingly powerful tools. A more prior question must now be asked: what gives a study the right to begin, and under what conditions does its underlying question genuinely stand?
2026, 47(8): 978-989.
doi: 10.21656/1000-0887.460204
Abstract:
A unified semi-analytical framework combining the numerical equivalent inclusion method (NEIM) and the distributed dislocation technique (DDT) was employed to investigate the interaction between surface cracks and arbitrarily shaped inhomogeneities in a half-plane. Several representative examples were simulated to systematically examine the effects of crack angles, inhomogeneity stiffnesses, and the inhomogeneity distances from the free surface on the stress intensity factors and local stress fields. The results show that, the crack angle leads to notable changes in the crack-tip field and induces pronounced mode coupling. Stiff inhomogeneities decrease Mode Ⅰ stress intensity factor but increase Mode Ⅱ, whereas soft inhomogeneities exhibit the opposite trend. As the inhomogeneity moves farther from the free surface, its perturbation to the crack-tip field becomes significantly weaker. The semi-analytical results agree well with the finite element solutions, with discrepancies below 4.5%, demonstrating the accuracy and efficiency of the present framework. The numerical examples provide clear insight into the dominant mechanisms governing crack-inhomogeneity interaction under multiple influencing parameters.
A unified semi-analytical framework combining the numerical equivalent inclusion method (NEIM) and the distributed dislocation technique (DDT) was employed to investigate the interaction between surface cracks and arbitrarily shaped inhomogeneities in a half-plane. Several representative examples were simulated to systematically examine the effects of crack angles, inhomogeneity stiffnesses, and the inhomogeneity distances from the free surface on the stress intensity factors and local stress fields. The results show that, the crack angle leads to notable changes in the crack-tip field and induces pronounced mode coupling. Stiff inhomogeneities decrease Mode Ⅰ stress intensity factor but increase Mode Ⅱ, whereas soft inhomogeneities exhibit the opposite trend. As the inhomogeneity moves farther from the free surface, its perturbation to the crack-tip field becomes significantly weaker. The semi-analytical results agree well with the finite element solutions, with discrepancies below 4.5%, demonstrating the accuracy and efficiency of the present framework. The numerical examples provide clear insight into the dominant mechanisms governing crack-inhomogeneity interaction under multiple influencing parameters.
2026, 47(8): 990-998.
doi: 10.21656/1000-0887.460118
Abstract:
To address the challenges of mechanical property prediction and structural optimization design for auxetic honeycomb materials, a machine learning approach integrating the particle swarm optimization (PSO) and the long short-term memory (LSTM) networks was proposed. A training dataset consisting of 400 groups of geometric parameters (including straight wall lengths, cell heights, wall thicknesses, and cell angles) and their corresponding mechanical properties (including energy absorption, Young’s moduli, and Poisson’s ratios) was constructed through finite element simulation. The PSO algorithm was employed to globally optimize the hyperparameters of the LSTM model, to establish a multi-objective mechanical property prediction model. Furthermore, an inverse design framework based on the PSO-LSTM was developed, enabling the optimization of geometric parameters driven by target mechanical properties. Experimental results show that, the optimized PSO-LSTM model achieves R2 values of 0.983 4, 0.974 6, and 0.970 4 for energy absorption, Young’s moduli, and Poisson’s ratios, respectively, with mean squared errors (MSE) below 0.001 2. The relative errors of total energy absorption, Young’s moduli, and Poisson’s ratios for the model obtained through inverse design are 0.432%, 1.05%, and 0.327%, respectively. The proposed method provides theoretical support for the intelligent design and engineering application of auxetic materials.
To address the challenges of mechanical property prediction and structural optimization design for auxetic honeycomb materials, a machine learning approach integrating the particle swarm optimization (PSO) and the long short-term memory (LSTM) networks was proposed. A training dataset consisting of 400 groups of geometric parameters (including straight wall lengths, cell heights, wall thicknesses, and cell angles) and their corresponding mechanical properties (including energy absorption, Young’s moduli, and Poisson’s ratios) was constructed through finite element simulation. The PSO algorithm was employed to globally optimize the hyperparameters of the LSTM model, to establish a multi-objective mechanical property prediction model. Furthermore, an inverse design framework based on the PSO-LSTM was developed, enabling the optimization of geometric parameters driven by target mechanical properties. Experimental results show that, the optimized PSO-LSTM model achieves R2 values of 0.983 4, 0.974 6, and 0.970 4 for energy absorption, Young’s moduli, and Poisson’s ratios, respectively, with mean squared errors (MSE) below 0.001 2. The relative errors of total energy absorption, Young’s moduli, and Poisson’s ratios for the model obtained through inverse design are 0.432%, 1.05%, and 0.327%, respectively. The proposed method provides theoretical support for the intelligent design and engineering application of auxetic materials.
2026, 47(8): 999-1008.
doi: 10.21656/1000-0887.460178
Abstract:
The trippingin speed affects the wellbore effective pressure, leading to changes in formation fracture widths and posing a risk of mud loss. The trippingin speed and drilling fluid compressibility, along with the well depth, the bottomThe trippingin speed affects the wellbore effective pressure, leading to changes in formation fracture widths and posing a risk of mud loss. The trippingin speed and drilling fluid compressibility, along with the well depth, the bottomhole assembly, the formation petrophysical parameters, and the drilling fluid properties were incorporated to build a mathematical model coupling wellbore pressure transients and fracture deformation. The model was solved with the finite elementfinite volume coupling method and validated with data from the hardbrittle shale formation of the Longmaxi formation in the Zi X well of the Weiyuan Block in the Sichuan Basin. The results show that, i. when the trippingin speed increase from 0.5 m/s to 2.0 m/s, the wellbore pressure will increase from 86.0 MPa to 106.6 MPa, while the fracture width will rise from 0.478 mm to 0.881 mm. ii. Bigger fracture widths go with greater trippingin depths and higher drill string annular ratios. At a trippingin speed of 0.5 m/s, when the tripping depth increases from 1 500 m to 5 500 m, the fracture width will rise from 0.432 mm to 0.478 mm; when the drill string annular ratio increases from 0.59 to 0.65, the fracture width will rise from 0.463 mm to 0.487 mm. iii. As the trippingin speed increases from 0 to 1.5 m/s, at a tripping depth of 1 500 m, the fracture width will rise from 0.4 mm to 0.474 mm; at a tripping depth of 5 500 m, the fracture width will rise from 0.5 mm to 0.624 mm. Bigger fracture widths go with higher trippingin speeds, and the rising trend will be more pronounced once the trippingin speed exceeds 1.0 m/s; further, this rise becomes more pronounced with greater tripping depths. This study provides a theoretical basis for optimizing trippingin operation parameters and offers an important guidance for preventing mud loss incidents during drilling.hole assembly, the formation petrophysical parameters, and the drilling fluid properties were incorporated to build a mathematical model coupling wellbore pressure transients and fracture deformation. The model was solved with the finite elementfinite volume coupling method and validated with data from the hardbrittle shale formation of the Longmaxi formation in the Zi X well of the Weiyuan Block in the Sichuan Basin. The results show that, i. when the trippingin speed increase from 0.5 m/s to 2.0 m/s, the wellbore pressure will increase from 86.0 MPa to 106.6 MPa, while the fracture width will rise from 0.478 mm to 0.881 mm. ii. Bigger fracture widths go with greater trippingin depths and higher drill string annular ratios. At a trippingin speed of 0.5 m/s, when the tripping depth increases from 1 500 m to 5 500 m, the fracture width will rise from 0.432 mm to 0.478 mm; when the drill string annular ratio increases from 0.59 to 0.65, the fracture width will rise from 0.463 mm to 0.487 mm. iii. As the trippingin speed increases from 0 to 1.5 m/s, at a tripping depth of 1 500 m, the fracture width will rise from 0.4 mm to 0.474 mm; at a tripping depth of 5 500 m, the fracture width will rise from 0.5 mm to 0.624 mm. Bigger fracture widths go with higher trippingin speeds, and the rising trend will be more pronounced once the trippingin speed exceeds 1.0 m/s; further, this rise becomes more pronounced with greater tripping depths. This study provides a theoretical basis for optimizing trippingin operation parameters and offers an important guidance for preventing mud loss incidents during drilling.
The trippingin speed affects the wellbore effective pressure, leading to changes in formation fracture widths and posing a risk of mud loss. The trippingin speed and drilling fluid compressibility, along with the well depth, the bottomThe trippingin speed affects the wellbore effective pressure, leading to changes in formation fracture widths and posing a risk of mud loss. The trippingin speed and drilling fluid compressibility, along with the well depth, the bottomhole assembly, the formation petrophysical parameters, and the drilling fluid properties were incorporated to build a mathematical model coupling wellbore pressure transients and fracture deformation. The model was solved with the finite elementfinite volume coupling method and validated with data from the hardbrittle shale formation of the Longmaxi formation in the Zi X well of the Weiyuan Block in the Sichuan Basin. The results show that, i. when the trippingin speed increase from 0.5 m/s to 2.0 m/s, the wellbore pressure will increase from 86.0 MPa to 106.6 MPa, while the fracture width will rise from 0.478 mm to 0.881 mm. ii. Bigger fracture widths go with greater trippingin depths and higher drill string annular ratios. At a trippingin speed of 0.5 m/s, when the tripping depth increases from 1 500 m to 5 500 m, the fracture width will rise from 0.432 mm to 0.478 mm; when the drill string annular ratio increases from 0.59 to 0.65, the fracture width will rise from 0.463 mm to 0.487 mm. iii. As the trippingin speed increases from 0 to 1.5 m/s, at a tripping depth of 1 500 m, the fracture width will rise from 0.4 mm to 0.474 mm; at a tripping depth of 5 500 m, the fracture width will rise from 0.5 mm to 0.624 mm. Bigger fracture widths go with higher trippingin speeds, and the rising trend will be more pronounced once the trippingin speed exceeds 1.0 m/s; further, this rise becomes more pronounced with greater tripping depths. This study provides a theoretical basis for optimizing trippingin operation parameters and offers an important guidance for preventing mud loss incidents during drilling.hole assembly, the formation petrophysical parameters, and the drilling fluid properties were incorporated to build a mathematical model coupling wellbore pressure transients and fracture deformation. The model was solved with the finite elementfinite volume coupling method and validated with data from the hardbrittle shale formation of the Longmaxi formation in the Zi X well of the Weiyuan Block in the Sichuan Basin. The results show that, i. when the trippingin speed increase from 0.5 m/s to 2.0 m/s, the wellbore pressure will increase from 86.0 MPa to 106.6 MPa, while the fracture width will rise from 0.478 mm to 0.881 mm. ii. Bigger fracture widths go with greater trippingin depths and higher drill string annular ratios. At a trippingin speed of 0.5 m/s, when the tripping depth increases from 1 500 m to 5 500 m, the fracture width will rise from 0.432 mm to 0.478 mm; when the drill string annular ratio increases from 0.59 to 0.65, the fracture width will rise from 0.463 mm to 0.487 mm. iii. As the trippingin speed increases from 0 to 1.5 m/s, at a tripping depth of 1 500 m, the fracture width will rise from 0.4 mm to 0.474 mm; at a tripping depth of 5 500 m, the fracture width will rise from 0.5 mm to 0.624 mm. Bigger fracture widths go with higher trippingin speeds, and the rising trend will be more pronounced once the trippingin speed exceeds 1.0 m/s; further, this rise becomes more pronounced with greater tripping depths. This study provides a theoretical basis for optimizing trippingin operation parameters and offers an important guidance for preventing mud loss incidents during drilling.
2026, 47(8): 1009-1018.
doi: 10.21656/1000-0887.460152
Abstract:
The reflection problem of elastic waves at the elastic interface in couplestress solids was studied under the influence of an external magnetic field. Firstly, based on the couple stress theory and Maxwell’s electromagnetic theory, the governing equations for the propagation of elastic waves and 3 types of waves (the P wave, the SV wave and the SS wave) were derived. It was also found that the external magnetic field affects the propagation of elastic waves through the Lorentz force, but no new wave patterns are generated. Secondly, based on the elastic interface conditions, the corresponding linear algebraic equations were obtained, and the amplitude ratios of various reflected waves relative to the incident wave in the case of incident P waves and SV waves were derived. Then, the reflection coefficient was defined with the energy flow ratio, and the effects of the external magnetic field, 3 elastic constants, couple stress parameters, and the angular frequency of the incident wave on the elastic wave reflection coefficient, were discussed for the incident P wave. Finally, the accuracy of the numerical results was verified through the conservation of normal energy for each wave.
The reflection problem of elastic waves at the elastic interface in couplestress solids was studied under the influence of an external magnetic field. Firstly, based on the couple stress theory and Maxwell’s electromagnetic theory, the governing equations for the propagation of elastic waves and 3 types of waves (the P wave, the SV wave and the SS wave) were derived. It was also found that the external magnetic field affects the propagation of elastic waves through the Lorentz force, but no new wave patterns are generated. Secondly, based on the elastic interface conditions, the corresponding linear algebraic equations were obtained, and the amplitude ratios of various reflected waves relative to the incident wave in the case of incident P waves and SV waves were derived. Then, the reflection coefficient was defined with the energy flow ratio, and the effects of the external magnetic field, 3 elastic constants, couple stress parameters, and the angular frequency of the incident wave on the elastic wave reflection coefficient, were discussed for the incident P wave. Finally, the accuracy of the numerical results was verified through the conservation of normal energy for each wave.
2026, 47(8): 1019-1034.
doi: 10.21656/1000-0887.460127
Abstract:
The damage mechanisms and ballistic performance of ultrahigh molecular weight polyethylene (UHMWPE) fiber reinforced composite/metal armor plates were investigated under combined blast and fragment impact. The combined loading experiments were conducted with composite projectiles. The typical failure modes of the armor plates at varying impact velocities were systematically analyzed. On this basis, the finite element method was employed to validate the experimental results. Furthermore, the key influence of the panel structure layout on its dynamic response was further explored. The results demonstrate that, in the configuration with UHMWPE at the front, the rigid constraint imposed by the steel plate restricts the large deformation capacity of the fiber layers. However, it enhances the stress diffusion effect. Consequently, this configuration exhibits superior blast resistance. Conversely, in the configuration with UHMWPE at the back, the UHMWPE layer fully exerts its potential of viscoelastic deformation. It dissipates the fragment kinetic energy primarily through large deformation, resulting in better penetration resistance. The difference in protective efficacy between these 2 configurations originates from the distinct constraint effects from the material stacking sequence. This constraint effect directly governs the energy distribution mechanisms and the evolution of failure modes within the composite structure.
The damage mechanisms and ballistic performance of ultrahigh molecular weight polyethylene (UHMWPE) fiber reinforced composite/metal armor plates were investigated under combined blast and fragment impact. The combined loading experiments were conducted with composite projectiles. The typical failure modes of the armor plates at varying impact velocities were systematically analyzed. On this basis, the finite element method was employed to validate the experimental results. Furthermore, the key influence of the panel structure layout on its dynamic response was further explored. The results demonstrate that, in the configuration with UHMWPE at the front, the rigid constraint imposed by the steel plate restricts the large deformation capacity of the fiber layers. However, it enhances the stress diffusion effect. Consequently, this configuration exhibits superior blast resistance. Conversely, in the configuration with UHMWPE at the back, the UHMWPE layer fully exerts its potential of viscoelastic deformation. It dissipates the fragment kinetic energy primarily through large deformation, resulting in better penetration resistance. The difference in protective efficacy between these 2 configurations originates from the distinct constraint effects from the material stacking sequence. This constraint effect directly governs the energy distribution mechanisms and the evolution of failure modes within the composite structure.
2026, 47(8): 1035-1044.
doi: 10.21656/1000-0887.460177
Abstract:
The in-service supercritical circulating fluidized bed boiler was investigated. Based on actual operational data, a numerical simulation model for the temperature field and the stress field was established. The effects of boundary conditions, including the external heat flux density of the tube wall and the convective heat transfer coefficient, on the maximum wall temperature and maximum thermal stress of the water-wall tube, were discussed. In addition, the influences of structural parameters, such as the fin thickness and the tube pitch, on the maximum temperature difference and maximum thermal stress between adjacent water-wall tubes, were analyzed. The simulation results indicate that, under a constant convective heat transfer coefficient, both the maximum wall temperature and the maximum thermal stress of the water-wall tube will increase with the external heat flux density. Conversely, when the external heat flux density is kept constant, the maximum wall temperature and maximum thermal stress will decrease as the convective heat transfer coefficient between the tube wall and the working medium increases. Notably, when the convective heat transfer coefficient decreases to 0.5 kW/((m2·K), the maximum thermal stress will increase sharply. Furthermore, the maximum temperature difference and maximum thermal stress between adjacent water-wall tubes will increase with the decrease of the fin thickness, but increase with the tube pitch. Specifically, when the fin thickness is reduced by 2 mm, the maximum temperature difference between adjacent water-wall tubes will increase by approximately 5℃, accompanied by an increase in the maximum thermal stress from 213 MPa to 219 MPa. When the tube pitch increases by 2 mm, the maximum temperature difference will rise by about 7 ℃, and the corresponding maximum thermal stress will rise from 213.8 MPa to 215.1 MPa.
The in-service supercritical circulating fluidized bed boiler was investigated. Based on actual operational data, a numerical simulation model for the temperature field and the stress field was established. The effects of boundary conditions, including the external heat flux density of the tube wall and the convective heat transfer coefficient, on the maximum wall temperature and maximum thermal stress of the water-wall tube, were discussed. In addition, the influences of structural parameters, such as the fin thickness and the tube pitch, on the maximum temperature difference and maximum thermal stress between adjacent water-wall tubes, were analyzed. The simulation results indicate that, under a constant convective heat transfer coefficient, both the maximum wall temperature and the maximum thermal stress of the water-wall tube will increase with the external heat flux density. Conversely, when the external heat flux density is kept constant, the maximum wall temperature and maximum thermal stress will decrease as the convective heat transfer coefficient between the tube wall and the working medium increases. Notably, when the convective heat transfer coefficient decreases to 0.5 kW/((m2·K), the maximum thermal stress will increase sharply. Furthermore, the maximum temperature difference and maximum thermal stress between adjacent water-wall tubes will increase with the decrease of the fin thickness, but increase with the tube pitch. Specifically, when the fin thickness is reduced by 2 mm, the maximum temperature difference between adjacent water-wall tubes will increase by approximately 5℃, accompanied by an increase in the maximum thermal stress from 213 MPa to 219 MPa. When the tube pitch increases by 2 mm, the maximum temperature difference will rise by about 7 ℃, and the corresponding maximum thermal stress will rise from 213.8 MPa to 215.1 MPa.
2026, 47(8): 1045-1059.
doi: 10.21656/1000-0887.460114
Abstract:
The pipeline for marine transportation of high-temperature oil and gas is affected by the convective heat transfer of the fluid inside the pipe and the internal heat source. With the rise of the pipe body temperature, the heat will also transfer from the pipe wall surface to the external low-temperature seawater, to form a temperature gradient in the pipeline radial and axial direction. The 2D transient heat transfer in an axisymmetric single-layer pipe with non-homogeneous boundary conditions and an internal heat source was analyzed with the integral transform technique. First, the 2D transient temperature was expressed as a combination of axisymmetric 1D steady-state temperature with a heat source, a filtered 2D steady-state temperature, and a filtered homogeneous 2D transient temperature, and was non-dimensionalized. Next, based on the method of separating variables to determine the integral transform pairs in radial and axial directions, the heat conduction control equations and boundary conditions plus initial conditions were subjected to the integral transform process to separate the time and space dependence of transient temperatures, and the 1st-order linear ordinary differential equations with respect to time were obtained. Finally, the theoretical solution of the 2D transient temperature distribution was obtained through the inverse transformation of the integral transform pair, and the 2D transient temperature distribution was verified through comparison of the results with the steady-state temperature distribution. On this basis, the effects of different combinations of inner and outer Biot numbers and different heat source strengths G on the temperature distributions of pipes were investigated. The results show that, the temperature gradients and temperature values in the spatial and temporal dimensions of the pipe are larger for the same combination of internal and external Biot numbers and for the high intensity of the heat source. There are differences in the temperature gradient and overall temperature distribution in the spatial and temporal dimensions when the intensity of the heat source is the same but the combinations of the internal and external Biot numbers are different. In the radial direction and the temporal dimension, the temperature distribution curves of different combinations of inner and outer Biot numbers have an intersection point where the temperature gradient and the relative magnitude of the temperature values change, but this phenomenon does not occur in the axial direction.
The pipeline for marine transportation of high-temperature oil and gas is affected by the convective heat transfer of the fluid inside the pipe and the internal heat source. With the rise of the pipe body temperature, the heat will also transfer from the pipe wall surface to the external low-temperature seawater, to form a temperature gradient in the pipeline radial and axial direction. The 2D transient heat transfer in an axisymmetric single-layer pipe with non-homogeneous boundary conditions and an internal heat source was analyzed with the integral transform technique. First, the 2D transient temperature was expressed as a combination of axisymmetric 1D steady-state temperature with a heat source, a filtered 2D steady-state temperature, and a filtered homogeneous 2D transient temperature, and was non-dimensionalized. Next, based on the method of separating variables to determine the integral transform pairs in radial and axial directions, the heat conduction control equations and boundary conditions plus initial conditions were subjected to the integral transform process to separate the time and space dependence of transient temperatures, and the 1st-order linear ordinary differential equations with respect to time were obtained. Finally, the theoretical solution of the 2D transient temperature distribution was obtained through the inverse transformation of the integral transform pair, and the 2D transient temperature distribution was verified through comparison of the results with the steady-state temperature distribution. On this basis, the effects of different combinations of inner and outer Biot numbers and different heat source strengths G on the temperature distributions of pipes were investigated. The results show that, the temperature gradients and temperature values in the spatial and temporal dimensions of the pipe are larger for the same combination of internal and external Biot numbers and for the high intensity of the heat source. There are differences in the temperature gradient and overall temperature distribution in the spatial and temporal dimensions when the intensity of the heat source is the same but the combinations of the internal and external Biot numbers are different. In the radial direction and the temporal dimension, the temperature distribution curves of different combinations of inner and outer Biot numbers have an intersection point where the temperature gradient and the relative magnitude of the temperature values change, but this phenomenon does not occur in the axial direction.
2026, 47(8): 1060-1067.
doi: 10.21656/1000-0887.460166
Abstract:
The bouncing ball model has garnered significant attention since the 1 940 s, with scholars extensively investigating the dynamics of a point mass undergoing collisions with a periodically oscillating surface under gravitational influence. As the chaos theory and nonsmooth dynamics evolved, this model emerged as a quintessential testbed for the application of these theoretical frameworks. Herein, the inquiry was extended to the bouncing ellipsoid model. The equations of motion were meticulously derived for an ellipsoid experiencing perfectly elastic collisions. The analysis yields the conserved quantities of the system and elucidates the stability conditions of its trivial periodic orbits. Furthermore, the numerical simulations corroborate the theoretical findings, thereby substantiating the robustness of the proposed method.
The bouncing ball model has garnered significant attention since the 1 940 s, with scholars extensively investigating the dynamics of a point mass undergoing collisions with a periodically oscillating surface under gravitational influence. As the chaos theory and nonsmooth dynamics evolved, this model emerged as a quintessential testbed for the application of these theoretical frameworks. Herein, the inquiry was extended to the bouncing ellipsoid model. The equations of motion were meticulously derived for an ellipsoid experiencing perfectly elastic collisions. The analysis yields the conserved quantities of the system and elucidates the stability conditions of its trivial periodic orbits. Furthermore, the numerical simulations corroborate the theoretical findings, thereby substantiating the robustness of the proposed method.
2026, 47(8): 1068-1085.
doi: 10.21656/1000-0887.450230
Abstract:
The dynamics equations for multi-body systems with closed-loop constraints are a set of differential algebraic mixed equations. The main difficulties in their numerical solution are as follow: it is difficult to satisfy the displacement, velocity and acceleration constraint equations with high precision at the same time, which will cause large defaults with the accumulation of errors; there are often redundant constraints and singular configuration problems, which seriously affect the numerical behavior of the solution process and the accuracy of the results. Herein, a numerical method for solving singular configuration closed-loop multibody systems was presented, without the traditional constraint independence assumption. The correction of position and velocity were done to satisfy the constraint equations before the motion equation formulation, rather than the correction of the constraint default at the end of each integration step. This method works with any standard ODE solver. According to the geometric characteristics of constraints, the components of velocity and acceleration in the hypersurface constraint space were all determined by constraints. The orthogonal basis vectors of the constrained tangent space were obtained through QR decomposition of the column permutation. The relationship between the generalized velocity of the system and the independent generalized velocity was established, and the default values of the generalized velocity and the generalized acceleration were filtered out in the dynamic calculation process, and the differential algebraic mixed equation was transformed into a pure differential equation satisfying the velocity and acceleration constraints ( with the number of equations equal to the generalized coordinate number of the system). Then, the orthogonal projection default correction method was applied to the system, and the default magnitude of the system position and speed was effectively suppressed by several iterations. At the same time, according to the principle of virtual power equivalence, the physical significance of the Laplace multiplier in solving the constraint reaction torque of the cut hinge was further revealed. Numerical examples show that, the precision of displacement, velocity and acceleration constraints of the proposed method is higher than that of the traditional augmented method and the penalty function method. The proposed method can deal with redundant constraints and singular configuration problems, and can be programmed to solve multi-body systems with closed-loop constraints with high precision. The proposed method can identify the independent constraints dynamically and modify the motion equation accordingly. Numerical examples demonstrate the effectiveness of the proposed method.
The dynamics equations for multi-body systems with closed-loop constraints are a set of differential algebraic mixed equations. The main difficulties in their numerical solution are as follow: it is difficult to satisfy the displacement, velocity and acceleration constraint equations with high precision at the same time, which will cause large defaults with the accumulation of errors; there are often redundant constraints and singular configuration problems, which seriously affect the numerical behavior of the solution process and the accuracy of the results. Herein, a numerical method for solving singular configuration closed-loop multibody systems was presented, without the traditional constraint independence assumption. The correction of position and velocity were done to satisfy the constraint equations before the motion equation formulation, rather than the correction of the constraint default at the end of each integration step. This method works with any standard ODE solver. According to the geometric characteristics of constraints, the components of velocity and acceleration in the hypersurface constraint space were all determined by constraints. The orthogonal basis vectors of the constrained tangent space were obtained through QR decomposition of the column permutation. The relationship between the generalized velocity of the system and the independent generalized velocity was established, and the default values of the generalized velocity and the generalized acceleration were filtered out in the dynamic calculation process, and the differential algebraic mixed equation was transformed into a pure differential equation satisfying the velocity and acceleration constraints ( with the number of equations equal to the generalized coordinate number of the system). Then, the orthogonal projection default correction method was applied to the system, and the default magnitude of the system position and speed was effectively suppressed by several iterations. At the same time, according to the principle of virtual power equivalence, the physical significance of the Laplace multiplier in solving the constraint reaction torque of the cut hinge was further revealed. Numerical examples show that, the precision of displacement, velocity and acceleration constraints of the proposed method is higher than that of the traditional augmented method and the penalty function method. The proposed method can deal with redundant constraints and singular configuration problems, and can be programmed to solve multi-body systems with closed-loop constraints with high precision. The proposed method can identify the independent constraints dynamically and modify the motion equation accordingly. Numerical examples demonstrate the effectiveness of the proposed method.
2026, 47(8): 1086-1092.
doi: 10.21656/1000-0887.460039
Abstract:
The rotational inertia caused by the shear effect was considered, and a modified axial motion Timoshenko beam model was established under Newton’s 2nd law in the material coordinate system. The natural frequencies of the system were solved, and the laws of the natural frequencies changing with parameters such as the beam length, the beam height and the axial velocity, were analyzed. The results obtained by this model were compared with those by the traditional axially moving Euler-Bernoulli beam theoretical model, the axially moving Rayleigh beam theoretical model and the axially moving Timoshenko beam theoretical model. The results show that, when the length-to-height ratio of the axially moving beam is greater than 10, the results obtained by the 4 models are consistent. When the length-to-height ratio is relatively small, the present results are basically consistent with those by the classical axially moving Timoshenko beam model. In addition, the rotational inertia caused by the shear effect has a greater impact on higher-order frequencies. Therefore, the influence of the moment of inertia due to the shear effect shall be considered in the dynamic analysis of axially moving beams.
The rotational inertia caused by the shear effect was considered, and a modified axial motion Timoshenko beam model was established under Newton’s 2nd law in the material coordinate system. The natural frequencies of the system were solved, and the laws of the natural frequencies changing with parameters such as the beam length, the beam height and the axial velocity, were analyzed. The results obtained by this model were compared with those by the traditional axially moving Euler-Bernoulli beam theoretical model, the axially moving Rayleigh beam theoretical model and the axially moving Timoshenko beam theoretical model. The results show that, when the length-to-height ratio of the axially moving beam is greater than 10, the results obtained by the 4 models are consistent. When the length-to-height ratio is relatively small, the present results are basically consistent with those by the classical axially moving Timoshenko beam model. In addition, the rotational inertia caused by the shear effect has a greater impact on higher-order frequencies. Therefore, the influence of the moment of inertia due to the shear effect shall be considered in the dynamic analysis of axially moving beams.
Download Center
More>
CopyRight
Founded in 1980 ( monthly )
Supervisor：Chongqing Jiaotong University
Founder：CHIEN Weizang
Editor-in-Chief：LU Tianjian
ISSN 1000-0887
CN 50-1060/O3
