Abstract
The elastic material properties which change momentarily and locally under the high deformation rate due to the movement of a wavefront are presented. The work contains mathematical formulation, semianalytical results, numerical formulations, and simulation results demonstrating the effectiveness of modifying the rheological properties of the elastic material upon shock load or contact with a rigid obstacle. While the semianalytical solutions can be obtained in a narrow time interval, numerical solutions allow us to track the process of wavefront reflections from edges. The effectiveness of reducing the physical quantities significant for impact in the presented examples reaches 30–70% of forces or accelerations, depending on the adopted criteria.
1. Introduction
Impact protection materials are a problem of great practical importance. Improving the parameters of materials responsible for protecting the human head, human body, or broadly understood engineering structures is highly desirable. Therefore, new solutions, mainly using intelligent materials, are constantly sought.
The aim of the work is to present the concept of an intelligent material reducing the effects of impact loads. The material used can be both rigid and effectively reduce effects of impacts. The improvement in efficiency is possible due to the local and temporary change in stiffness induced by the moving mechanical wave. It should be emphasized that the designed material should have a high usable stiffness and only in a critical moment should locally and temporarily change its properties. In this way, the usability of an element or a structure with the use of intelligent material will remain unchanged throughout the entire range of use, including the moment of impacts.
The crucial significance of mechanical properties in practical applications has made it one of the most fundamental and extensively researched areas in materials science. Sudden impacts cause significant harm to lives and devices. Thus, the development of crash-resistant devices plays an essential role in engineering applications such as aircraft, vessels, and automobiles. Hence, the demand for the safety of structures against severe loading is increasing day by day as structure safety is generally associated with the materials used for construction [1–3].
In practical applications, materials that absorb impact energy using the phenomenon of plastic deformations are most often selected [4–6]. The use of aluminium or steel enables achieving desired properties such as lightweight, strength, and the ability for controlled deformation, which is crucial for the effective absorption of kinetic energy in the case of impacts. During plastic deformation, the material experiences tensile, shear, compressive, or torsional stresses, depending on the type of loading conditions [7–9]. Consequently, the kinetic energy is dissipated internally through mechanisms such as grain boundary sliding, dislocation motion, or fracture propagation within the material’s microstructure. This energy absorption process mitigates accelerations and forces exerted on the structure, thereby minimizing the risk of injuries or structural damage.
Another engineering material that finds application in reducing accelerations is the carbon fiber composite. These composites consist of a polymer resin reinforced with carbon or glass fibers. They are significantly lighter than steel while exhibiting high strength and stiffness. Carbon fiber-reinforced polymer utilizes brittle fracturing to reduce the effects of impact. Carbon fibers are incredibly strong and stiff but they are also brittle. When subjected to dynamic loading, such as during an impact, carbon fibers experience sudden fracturing, resulting in the absorption and dispersion of kinetic energy [10–12]. This process reduces accelerations and forces acting on the structure, providing protection against damage.
There are also attempts to use the so-called smart materials in acceleration reduction, especially magnetorheological materials. Magnetorheological (MR) materials are used to reduce the effects of impact by leveraging their ability to change viscosity or hardness in response to a magnetic field. Under dynamic loading, MR materials can adapt their rheological properties, allowing for effective absorption of kinetic energy and reducing accelerations and forces acting on the structure. These are interesting solutions, however, requiring additional control devices that activate the material. This is a major limitation in practical applications [13–16].
Also, cellular materials, such as honeycombs and foams, have been studied for impact energy absorption. The undoubted advantages of this type of materials are low density and high stiffness. Due to the geometric complexity of the material’s structure, experimental research is carried out mainly. Examples include hybrid nickel-strengthened aluminium foams. Hybrid metallic foams are highly strain-rate sensitive, so they are able to dissipate more energy under dynamic loads [17]. In [18], various recipes for energy-absorbing granular materials were experimentally tested. The obtained physicomechanical properties are presented, namely, bulk density, flexural, compressive, and impact strength. Based on these results, the rate and mode of energy absorption can be predicted. There are also attempts to find cellular structures inspired by nature [19]. The authors presented how the mechanical properties of metal foam can be improved by changing its structure at different hierarchical levels, inspired by features and principles important for the impact resistance of biological patterns. An interesting example of cellular material is the metallic hollow sphere material. The mechanical behaviors were experimentally characterized [20–22]. Functionally graded porous materials are also being developed. It is characterized by a gradual change in the composition and structure of the volume, which results in corresponding changes in the material properties. The specific functions are attempted to be modeled using the finite element approach [23, 24].
A promising area of research is mechanical metamaterials, which, thanks to their nontrivial structure, have interesting properties [25]. These are innovative materials that are used to reduce the effects of impact. They possess a unique structure composed of repeating units with different mechanical properties. These units, known as metamaterial cells, are optimized to absorb and dissipate kinetic energy, thereby reducing accelerations and forces acting on the structure. By precisely manipulating mechanical properties such as stiffness, elasticity, and damping, mechanical metamaterials offer the ability to design custom structures that minimize the impact effects [26–31]. The numerical search for appropriate structures requires time-consuming simulations of 3D objects with a large number of degrees of freedom. In addition, in the case of structural dynamics computations, the simulation time is significantly longer. Therefore, the direct search (optimization) of metamaterial structures is extremely expensive and practically impossible. Recent progress in the development of energy-absorbing metamaterials and their mechanical properties under dynamic loading is presented in [32].
The paper proposes a simplified analysis of the behavior of a hypothetical metamaterial. The material is elastic, and the mechanical wave traveling in it causes a local change in the stiffness of the material. The influence of local reduction or increase of Young’s modulus on the response of the structure was investigated. The information obtained indicates the way to search for the structure of the appropriate metamaterial that will be able to implement the assumed strategy of operation. A simple 1D model enables deep exploration of the problem due to its simplicity and relatively short computation time.
The remaining sections of this paper are as follows. The mathematical model is presented in Section 2. To model the traveling stiffness zone, we use the Heaviside functions. This leads to a nonconstant longitudinal force in a rod. The problem definition and equation of motion are given. The semianalytical solution to the model is presented in Section 3 based on a separation of variables method and integral transformation. A numerical finite element method scheme for the problem is constructed in Section 4. The obtained results confirm the correctness of the developed numerical model. The influence of the number of local material changes on the system’s response is depicted in Section 5. The paper concludes in Section 6.
2. Mathematical model
2.1. Problem Definition
Let us consider a homogeneous straight clamped-free bar of length and a constant cross section clamped on the left side (Figure 1).
The (x, t) function defines longitudinal displacements. So, we can write the boundary conditions as follows:
Assume that the material of the rod is linearly elastic and thus follows Hook’s law. Moreover, the bar has a density and a global stiffness . The following zero initial displacement and the initial velocity are obtained:
In addition, we assume that along with the wave moving in the rod, its stiffness changes locally.where the wave speed is given by the following formula:
We will consider the following two cases: local strengthening and local weakening of stiffness. The difference of the Heaviside functions defines the traveling stiffness zone of width . At this point, it should be noted that we solve semianalytically only the stage of the wave reaching the support .
2.2. Equation of Motion
In order to obtain the equation of free vibrations of a rod, we will examine the motion of an infinitesimally small element created by cutting the rod with two very close cross sections. The longitudinal force balance equation takes the following form:
Due to its elasticity, the longitudinal force in the rod can be determined depending on the deformation, which in turn depends on the displacement.
According to (3), the longitudinal force (6) is not constant, so the increase in longitudinal force can be written in the following form:
Finally, after rearrangements of equation (5), taking into account (7), and dividing by , it takes the form of the equation of motion.
It is a partial differential equation with variable coefficients. Since the variables are Dirac deltas and Heaviside steps, the potential solution is a solution in the distribution sense.
3. Semianalytical Solution
In order to solve equation (8) along with the boundary conditions (1) and the initial conditions (2), we will use the separation of variables method as follows:
Substitution of (1) into formula (9) leads to boundary conditions for the function . To solve the boundary problem, we define eigenfunctions and eigenvalues. The assumed boundary conditions are satisfied by the following series of sines:
So, we get an infinite sequence of eigenvalues. Each eigenvalue corresponds to an eigenfunction. The function is determined from the following formula:
The integral (11) is called the Fourier transformation. Taking into account (11), initial conditions (2) take the following form:
According to (11) and (10), the equation of motion (8) after rearrangements can be written in the following form:where the integrals are presented for the interval in the following form:
The second integral differs by the integration interval when . We can see that the solution (14) is only possible for . In the case of , we have to use l’Hôpital’s rule. The formula is given in the following form:
We limit the series to a finite number of terms. Finally, the system of equation (13) can be written in matrix form as follows:
The final solution is obtained by substituting the components of the vector into the series (10). The presented system of differential equation (16) will be solved using the numerical Newmark integration method with initial conditions (12). The self-prepared computer program was written in the Fortran programming language. The advantage of the obtained semianalytical solution is the possibility of using it to verify the solution using the finite element method. The disadvantage of the semianalytical solution is that the mathematical model is limited only to the first phase of the process, i.e., to the wave reaching the support . This is due to the mathematical structure of the function , i.e., the moving stiffness (3), which does not provide for wave reflection from the fixed end of the rod; therefore, further observation of the solution is unphysical. Taking into account the wave reflection from the support in the function would be extremely complicated mathematically. However, it is sufficient to correctly verify the developed numerical model, which will be used later in the paper. In addition, it should be remembered that the semianalytical solution applies only to selected boundary conditions (1).
4. Finite Element Model
In further analysis of the material with a local stiffness change zone, we will use the finite element method. The rod will be divided into finite elements of length . The considered moving zone of length may include only 1 element but also several elements. Figure 2 shows the case of the two involved elements.
We are dealing with the wave equation, so we will use linear shape functions describing the distribution of nodal displacements in a finite element as follows:and similarly, virtual displacements as follows:
Classical characteristic matrices of the rod can be written in the following form:
Moving components require our attention. According to (8), virtual energy of the traveling zone can be presented in the following form:
The abovementioned integrals can be computed considering (17) and (18) as follows:
In the case of , the same methodology was used.
As a result of the minimization of the energy (21), we obtain the following vector:where is a parameter that describes the current position of the moving component in a finite element.
The vector (23) describes the wave direction from right to left. When the zone of changed stiffness moves from left to right, the vector takes the following form:
According to (19) and (20), we build a global system of differential equations. Finally, we obtain a matrix diagram of the rod dynamics with a moving stiffness zone as follows:
This discrete problem can be solved, for example, by the Newmark method. The abovementioned numerical scheme was hand coded in the Fortran programming language. The developed computer program has an undeniable advantage. It takes into account the moving stiffness in both directions according to the vectors (23) and (25). In the following, a test example that illustrates the correctness of the numerical model is presented. It should be noted that the solution applies only to the arrival of the wave at the support . Figures 3 and 4 present the obtained results compared with the semianalytical solution from Section 3.
(a)
(b)
(a)
(b)
The dataset for a rod used in the simulations was as follows:(i)Young’s modulus: (ii)Moving Young’s modulus: (iii)Density of a rod: (iv)Length of a rod: (v)Width of a moving zone: (vi)Initial velocity: .
The semianalytical solution presented above shows the correctness of the results obtained by the numerical method. Successively increasing the number of terms of the series (10), which is the solution of the semianalytical solution, converges to the exact solution. This suggests that the semianalytic solution is correct. Unfortunately, direct mathematical proof of the correctness of the solution is not possible due to its final nature.
The differences in the semianalytical solution and the finite element method result from their nature. In the finite element solution, we use a mesh of elements. In simulations of structural dynamics concerning wave equations such as a rod, during the process investigation, waves are reflected between the mesh nodes. This causes additional oscillations to appear in the solution. We will never obtain an identical solution to the same solution with both methods, especially when we study displacement velocities (Figure 4). While in the semianalytical solution, the scope of considerations is limited to a particular case, numerical solutions offer much wider possibilities, allowing the observation of reflections of moving stiffness.
5. Results
In this work, the aim was to address the question of the effectiveness of locally modifying material properties under impact conditions. Commercial computational packages were not suitable for this task because they do not allow for accurate modeling of the problem. In the numerical model, it was necessary to consider the fact that a moving narrow zone of material softening/hardening, influenced by a propagating wave, simultaneously changes the wave propagation velocity in that region. Consequently, it alters the spatial location of the area where material parameters change. Dealing with variable coefficients, the problem becomes nonlinear. Therefore, a custom computer program was developed for simulations and its key components are described. Another aspect that had to be taken into account in the numerical simulations was the need to incorporate certain imperfections of the physical models in which local material weakening occurs in a gradual rather than abrupt manner. For this reason, the stated objective could only be achieved using a custom computational program.
Another important goal of numerical calculations was to assess the influence of the spatial width of the weakening zone and the value of the elasticity modulus of the weakened zone relative to the original modulus of elasticity. The own computer program was used for optimizing these parameters while maximizing the reduction of instantaneous accelerations as the objective function.
In the following, we will discuss the results obtained numerically using the matrices described in Section 4 in the case of a bar with one end fixed and the free end loaded with a force impulse. We will consider the acceleration of the free end of the rod and the axial force appearing in the finite element at the free end. The influence of the width of the hardening/softening zone will be exhibited as well as the ratio of the contribution of to the initial elasticity modulus .
Figure 5 presents the force and acceleration at the subjected end for the softened and stiffened traveling zone. The following data were assumed: the ratio = 0.9 and the width of the traveling gap equal to 1% of the length of the bar. In the case of the softening material, the axial force appearing in the edge element just after the impact is about 15% of the value related to the force in the unmodified material. In the same example, the hardening material exhibits about 60% higher force. The acceleration registered at the end of the rod with the temporarily weakening material is 25% less than the reference acceleration and 30% greater for the hardening material.
(a)
(b)
Figure 6 shows the influence of the amount of local material weakening on the response of the system. The range of weakening of the modulus of elasticity was assumed in the range from 0 to 95% of the value of the basic modulus of elasticity in the rod . In practice, the physical modulus in the moving zone changes from to 0.05 . It can be seen that the reduction of the internal force at the loaded end member is proportional to the ratio.
Figure 7 shows displacements in time for various relative additional elasticity module . Displacements are related to the displacement of the free end with unmodified elasticity module . In the analysed example, = 1. A slower but deeper increase in time of the resultant displacement value is visible in the case of softened material. This allows one to reduce the recorded accelerations and thus makes the impacts gentler.
The diagram in Figure 8 displays the internal force value measured at the free end over time. The horizontal sections of the diagram correspond to the increasing width of the modified gap zone. The figure illustrates that the width of the moving softened gap can be quite narrow. Although the wider gap slightly improves the results in relation to the narrow gap and reduces the acting forces and accelerations, it happens to a slight extent. This is theoretically explicable because the wavefront in the elastic material, in the absence of parabolization of the governing equation and blurring of the response due to damping, is relatively sharply outlined. Therefore, only a narrow zone is involved in the process described in the work.
Figure 9 depicts the amplitude in terms of the width of the modified moving zone. The increasing amplitude along with the lengthening of the first period of vibration results in a reduced value of the acceleration at the free end of the rod. The increase in the period of free-end oscillation is directly proportional to the width of the softening gap. This relationship is also evident in Figure 7, which depicts the impact of the softening modulus range. The extension of the oscillation period, even with increasing amplitude, leads to a reduction in the acceleration amplitudes. Figure 9 shows a threefold increase in displacement amplitudes, with a tenfold increase in the time it takes to return the observed point to its equilibrium position. This is the reason for the significant decrease in accelerations and forces. It should be emphasized again that the narrow, softened moving zone somehow distributes the impulse energy over the area of the object.
6. Conclusion
The paper shows that the locally and temporarily weakening material significantly reduces critical utility values, such as peak accelerations and peak force values acting at the edge of the tested object. The material weakens in the zone of high deformation velocities, i.e., in the location of the traveling wave front resulting from the impact. Only a narrow layer of the material is weakened, while the rest does not change its mechanical and functional properties.
Mathematical formulation and semianalytical solution of the task allowed us to verify the developed numerical model of finite elements. In the case of numerical modeling, the correctness of the description of the passage of the front of the weakened zone through the area of successive finite elements is important. The same applies to the end-of-zone transition of a weakened zone, where part of the finite element area is weakened and parts are not. In the numerical model, the direction of the elastic wave is important. Semianalytical and numerical results were consistent within the scope possible to be assessed with the appearance of known artifacts.
The simulation results show that material with locally weakening stiffness allows avoiding large and short-term acceleration peaks more effectively and can be successfully used in protection against accidents. The shown examples show that the reduction of forces or accelerations can range from 30 to 70% of the value of models with ordinary material. In further research, it is necessary to correlate the mathematical model with a group of real materials and identify directions for developing material compositions that meet theoretical expectations. In addition, it is essential to verify the effectiveness of impact energy absorption in two- and three-dimensional systems. It is worth investigating the potential effectiveness of multilayer materials in which individual layers would serve different mechanical functions. Properly selected relationships between the gap width and the relative additional elasticity module significantly reduce the forces acting during the impact. The obtained results are an important clue in the search for a metamaterial structure that implements the strategy of local stiffness change.
Abbreviations
| : | Density of the rod |
| : | Rod cross section area |
| : | Length of the rod |
| : | Global stiffness of the rod (Young’s modulus) |
| : | Local moving stiffness |
| : | Width of the moving stiffness zone |
| : | Wave speed in the rod |
| : | Initial velocity |
| : | Longitudinal displacements of the rod |
| : | Static longitudinal displacements of the rod |
| : | Heaviside function |
| : | Dirac delta function |
| : | Time |
| : | Longitudinal force in the rod |
| : | Fourier finite integral transformation of the function |
| : | Vector of subsequent integral transformation |
| : | Inertial matrix in semianalytical solution |
| : | Stiffness matrix in semianalytical solution |
| : | Finite element length |
| : | Displacement in the node |
| : | Virtual displacements |
| : | Virtual energy of the traveling zone |
| : | Current position of the moving component |
| : | Inertial matrix in the finite element approach |
| : | Stiffness matrix in the finite element approach |
| : | Vector of the moving zone |
| : | Nodal displacement vector |
| ., ..: | Dots over the letter denote derivatives with respect to time. |
Data Availability
The data used to support the findings of this study are available on request from the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported within the projects (UMO-2017/26/E/ST8/00532) and (UMO-2019/33/B/ST8/02686) funded by the Polish National Science Centre.