Prediction of mechanical-hysteresis behavior and complex moduli using the phase field crystal method with modified pressure controlled dynamic equation

Damping materials have been used in numerous engineering applications. The important property that plays a key role in this type of material is a damping capacity which is related to mechanical-hysteresis behavior of viscoelasticity. During the last decade, the phase-field-crystal (PFC) model has emerged as a robust tool to predict various material phenomena. This density-functional-type model has the advantage over a conventional phase-field model in terms of atomistic resolution and molecular dynamics in terms of computational expense. In this work, we propose a method to predict mechanical-hysteresis behavior and its related complex moduli parameters using a modified-pressure-controlled dynamics (MPCD) equation incorporating with the PFC model, denoted as PFC-MPCD method, in a three-dimensional solid structure. We modify the previously proposed pressure-controlled dynamics (PCD) equation by introducing the pressure-time derivative term which allows us to produce the results consistent with the standard linear solid model (SLS) at the broader frequency range. We apply sinusoidal pressure oscillation and compare the results predicted by both models. The results demonstrate that mechanical-hysteresis behavior and complex moduli parameters predicted by PFC-MPCD method are in good agreement with those of SLS model and consistent with numerous experimental observations whereas the results produced by original PCD equation tends to exhibit inaccurate results at the very frequency. We expect that this new PFC-MPCD method can be extended to a more complex system and still be capable to exhibit the accurate mechanical-hysteresis behavior and complex moduli parameters which results in predicting more reliable damping capacity parameter in damping material design.


Introduction
Damping materials have been used in many engineering applications and, in particular, in energy absorption such as vibration isolation and impact absorber. One of the most important properties that plays a key role in energy absorption is the damping capacity of material. This parameter indicates the ability of material in dissipating energy during the deformation which can be measured through mechanical-hysteresis curve of viscoelasticity. Viscoelasticity is the mechanical behavior which is characterized by three phenomena; viscoelastic-creep, stress relaxation and mechanical-hysteresis behavior. To select the appropriate materials in specific tasks, therefore, we need to predict these properties accurately and know the suitable working condition of these damping materials.
In the past few decades, computational materials science has emerged as an alternative choice to predict various material properties and phenomena. This new technique can help speeding up the design process and reducing the cost of the real experimental study. To predict this damping capacity accurately, we need a 2. Background

PFC-PCD method
The non-dimensional PFC free energy functional form employed in our study is given as [1,2]  where ò can be interpreted as the degree of undercooling which is inversely proportional to temperature and ρ is the density field variable. The evolution equation for ρ is the Cahn-Hilliard type which involves dissipative dynamics and mass conservation [2] r m m r r r r ¶ ¶ =  = - where L μ is the mobility coefficient and m d dr =  is the diffusion potential. This equation can be used to simulate evolution of ρ under specified (or fixed) temperature (from the value of ò), mass (from the value of r ), and volume (from value of grid spacing); typical simulations under this condition are solidification and microstructural evolution. One can also vary the grid spacing to simulate the evolution of the density profile under specified time-dependent deformation [16]. In the context of the viscoelastic behavior, this technique would be appropriate for the stress-relaxation calculation; however for the viscoelastic-creep and mechanicalhysteresis studies, where the system is subjected to specified external stress, additional dynamic equation is needed.
The PCD equation was first introduced by Kocher and Provatas [17]. The advantage of this equation is that it enables us to control internal pressure  int by specifying the external pressure. In their work, this equation was employed to investigate the full spectrum of solid-liquid-vapor transition problem. This PCD equation which refers to the original PCD method can be written as below where V is system volume which changes with time,  ext is the externally applied pressure, ¢ L def is the mobility coefficient and the internal pressure  int can be defined by This PCD equation when incorporates CH model equation (2) which adopts PFC free energy functional equation (1) will refer to as PFC-PCD method . We employed this method to investigate viscoelastic-creep behavior for one-dimensional crystal in our previous studies [18,19] and extended to three dimensional solid structure [20]. According to the results, although the PFC-PCD method predicted a reliable viscoelastic-creep behavior [18][19][20], however, it could not produce a mechanical-hysteresis behavior correctly, in particular, at high pressure oscillation frequency range [20] which might cause a problem in predicting a suitable working condition for damping materials.

Standard linear solid model
The standard linear solid (SLS) is a model of modeling viscoelastic behavior of viscoelastic material known as Zener model [22][23][24]. It exhibits more realistic behavior involving single exponent in both viscoelastic-creep and stress-relaxation. Furthermore, this model predicts the finite storage moduli at high frequency range and providing peak locations for both loss factor parameters. The simplified governing equation of the SLS model is given as below where e is strain response predicted by the model. ¢  is external applied pressure. K′ is the elastic coefficient which represents the elastic part of the model. τ′ is the shape factor coefficient which enables us to adjust the peak magnitude of all loss factors. To exhibit mechanical-hysteresis behavior, the sinusoidal pressure oscillation is required. So, the ¢  can be specify as where ¢  A is the amplitude of sinusoidal pressure oscillation. The ω=2πf is angular frequency and f is the frequency of sinusoidal pressure oscillation input. The t variable is time domain. Substitute equation (6) in equation (5), it becomes Then, equation (7) yields the solution below where E A is the strain oscillation amplitude. The δ parameter is phase angle which defines time-delayed response due to viscoelastic effect. Both parameters are functions of frequency as shown below , a r c t a n 1 9 The strain responses e(t) corresponding to pressure input ¢  t ( ) at various sinusoidal pressure oscillation frequencies are illustrated in figures 1(a)-(c). The figures demonstrates that the amplitude and phase angle are affected by pressure-oscillation frequency. The governing equations which describe these changes are given by equation (9) and illustrated in figures 2(a) and (b). Another important aspect is mechanical-hysteresis behavior curve showed in figure 2(c). It is established according to e(t) versus ¢  t ( ) based on equations (6) and (8). This curve provides us with the information about the energy dissipation during each oscillation cycle at specific frequency. Furthermore, this curve also contributes us the viscoelastic-measure-parameter called complex moduli parameters. These parameters consist of storage moduli (E′), loss moduli (E″) and loss tangent factor (tan(δ)). Each parameter is frequency dependent function governed by equation (10) (Derivation of complex moduli parameters is provided in appendix).
The illustration of these parameters are given in figure 2(d). The E′ parameter provides the measure of stored energy which represents the elastic part of the material. Whereas, the E″ parameter provides the energy dissipation or heat loss during the deformation. This parameter depicts the viscous portion of the material and d tan( ) parameter is defined by the ratio between loss modulus and storage modulus; , which provides a measure of damping in the material. If this parameter is greater than unity, the system behaves liquidlike material while the opposite case the system exhibits solid-like material.

.1. MPCD equation
Mechanical-hysteresis is one of the behaviors which characterizes viscoelasticity. This behavior requires sinusoidal pressure oscillation as an input to establish mechanical-hysteresis curve. According to our previous study [20], we found the existing PFC-PCD method predicted mechanical-hysteresis curve inaccurately, in particular, at high pressure oscillation frequency when compared to the SLS model. With regard to the Taylor expansion [20], we observed that the PCD equation closely resembled Kelvin-Voigt model which can only exhibit viscoelastic-creep but not mechanical-hysteresis in solid material, in particular, at high frequency range. This inaccuracy arose from the missing pressure time derivative term in the existing PCD equation. Therefore, we have introduced this pressure time derivative term to the existing equation. We refer to this new modified PCD version as MPCD equation as shown below is the time pressure derivative term which is introduced to the existing PCD equation. τ is the shape factor coefficient that enables us to adjust the magnitude of peak in the loss factor parameters. This additional term helps MPCD equation in producing mechanical-hysteresis behavior consistent with SLS model, in particular, at high pressure oscillation frequency where the existing version is failed. It should be noted that the MPCD equation is established according to the SLS model (Zener model) which restricts the model to describing viscoelastic-creep and relaxation in exponential manner [22]. Therefore, this limitation also applies to the MPCD equation.

Approximated dynamics behavior of MPCD equation
In this section, we aim to demonstrate that the approximation behavior of MPCD equation is consistent with SLS model. We begin with the definition of where V(t) is current volume at time t. L a is the length of unit cell in the reference state. L(t) is the length of unit cell at time t. e(t) is strain at time t. It is obvious that we can alter volume time derivative term to strain time derivative term by chain's rule. where dV/dt=3L 2 , dL/de=L a . Substitute equation (13) into equation (11), we obtain Since  int changes its value with e(t), therefore,  int can be written as a function of e(t). However, if the small deformation is postulated around the undeformed state, e 0 =0,  int can be written in Taylor expansion form [32] as ( )only at first derivative expansion, equation (16) yields Substitute equation (17) into equation (14), we obtain Regarding the binomial series of Ref. [32]. Consider only the first term of the expansion; (1+Δe) −2 ≈1−2Δe. Substitute this relation into equation (18). We obtain Equation (20) has the structure like SLS model equation (5) where the pressure time derivative term is included on the right hand side of the equation. Nonetheless, the viscoelastic behavior stems from the MPCD equation not from PFC free energy. According to equation (20), it is obvious that the PFC free energy only influences the  K variable, the second term on the left-hand side of equation (20), but not the structure of the equation. Therefore, the alternative choices of PFC free energy equation (1) can be applied and the PFC-MPCD method should still produce similar viscoelastic behavior correspond to those PFC free energy choices.

Method of simulation
For simplicity in all numerical experiments, the PFC free energy functional used in this study is given by equation (1); it should be noted that MPCD equation can be used with any free energy functional forms which is not limited to equation (1). This free energy functional formulation is minimized by periodic density field of a body-centered-cubic (BCC) crystal structure.
where q a =2π/L a . The value of q a can be determined from minimization of the free energy density [33] leading to = q 1 2 a and, consequently, p = L 2 2 a . One can also approximate the ranges of ò and r where the BCC solid is stable. This is achieved by calculating the free energy density of the solids using the one-mode approximations for stripes, hexagonal, and BCC structures and the free energy of the liquid using a uniform density. Then the phase diagram can be determined from the common tangent construction and the ranges of ò and r where the BCC solid is stable can be obtained [34].
In all simulations, we construct a BCC unit-cell with a periodic N×N×N grid where N=16 and the lattice parameter This results in the reference grid spacing, p D = r 2 8 0 in all dimensions, and the reference volume The initial density profile is constructed with the one-mode approximation in equation (21), where the reference density, r 0 , is specified. The profile is then relaxed using the evolution equation (equation (2)), while fixing Δr=Δr 0 and r r = 0¯, until the equilibrium density profile is reached. At this state, the reference internal pressure,  0 , can be calculated using equation (4). Also, the numerical method employed to solve the evolution equation (equation (2)), PCD equation (equation (3)), and the MPCD equation (equation (11)) is the Fourier pseudo-spectral method where Fourier transform is used to calculate spatial derivatives and the forward finite difference scheme is used to discretize the time derivative. For the viscoelastic simulations, both PCD and MPCD equations are employed where pressure sinusoidal oscillation profile of  t ext ( ) is specified and the change in the volume is used to update the grid spacing through the relation We note that the grid spacing Δr(t) is identical in all directions due to the assumption of crystal under hydrostatic stress. During the simulation, the mass conservation is imposed by changing the average density through Also, the mechanical equilibrium is maintained throughout the deformation process; this is accomplished by relaxing the system using the evolution equation at every time step of the time evolution from the PCD and MPCD equations. The system response is characterized by the strain defined by which is identical in all directions. The profiles of e(t) and  t ext ( ) is then used to quantify the viscoelastic behavior predicted form the PFC-PCD and PFC-MPCD methods. In hysteresis simulation, the applied pressure is set to where  Ā is the pressure amplitude. The pressure deviation from the reference pressure will then be With this type of pressure function, the system is then allowed to evolve temporally until the steady-stead strainresponse is reached. The strain response profile at steady state will exhibit the sinusoidal behavior similar to . We note that in all hysteresis simulations,  Ā is set to  0.02 0 . Finally, all parameters obtained from numerical study according to PFC-PCD and PFC-MPCD methods have been thoroughly compared and discussed. It should be noted that the values of r 0 and ò in all simulations are restricted at −0.260 and 0.175, respectively. These values of r 0 and ò ensure the stability of the BCC crystal as shown in the phase diagram provided by Wu et al [34].

Strain response and mechanical hysteresis behavior
In this part, we demonstrate the strain response under sinusoidal pressure oscillation and its mechanicalhysteresis behavior exhibited by PFC-MPCD and PFC-PCD methods. The steady-state strain responses regarding the PFC-MPCD method are given figures 4(a)-(c) whereas the strain responses regarding the PFC-PCD method are shown in figures 4(d)-(f) for comparison. According the numerical results at oscillation frequency f=0.004 and f=0.1, the steady-state strain response predicted by both methods are in good agreement to each other. Both methods demonstrate the strain amplitude reduction and increased phase lag (phase angle) at elevated oscillation frequency. These results are qualitatively consistent with those of SLS model prediction in figures 1 (a)-(b). However, at oscillation frequency f=1, both methods seem to behave differently. The amplitude of strain response exhibited by PFC-MPCD method tends to decline to the finite value while phase lag between strain response and pressure oscillation predicted by this method seems to be inphase at this frequency. On the contrary, the amplitude of strain response exhibited by PFC-PCD method tends to vanish. Moreover, the phase lag between the strain response and pressure oscillation predicted by the PFC-PCD method becomes larger at this frequency which is contrast to SLS model prediction as shown figure 1(c). To obtain more insight into these phenomena, we study the impact of frequency f on amplitude ratio A r and phase lag (δ) exhibited by PFC-MPCD and PFC-PCD methods as indicated in figures 5(a)-(b). It is obvious that both methods exhibit amplitude ratio reduction with increasing frequency. However, this prediction according to the PFC-PCD method tends to vanishingly small at high-frequency range whereas the prediction regarding PFC-MPCD method remains at finite value. Besides, the phase lag predicted by both methods tends to increase at elevated frequency but until the certain value is reached. In particular, the result exhibited by PFC-MPCD method tends to reduce its magnitude with increasing frequency after its peak location, unlike PFC-PCD method prediction which tends to exhibit higher phase lag without the peak location. These results obtained from the numerical study indicate that the prediction demonstrated by both methods are consistent at lowfrequency range but tends to differ at high frequency.
The magnitude of the amplitude ratio can also be illustrated in the mechanical-hysteresis loop according to figures 6(a)-(b) where the results in figures 4(a)-(f) are plotted on a pressure-strain axis. As the frequency  increases, the hysteresis loops predicted by both methods rotate counter-clockwise, indicating smaller strain amplitude (or amplitude ratio). At frequency f=1, the hysteresis loop predicted by the PFC-PCD method almost forms a vertical line with the y-axis in figure 6(b) because the amplitude ratio is driven to an extremely small value or almost zero amplitude. This is in contrast to the hysteresis loop predicted by PFC-MPCD method which remains tilted away from the y-axis as indicated by figure 6(a); this figure shows that strain amplitude has not been driven to zero amplitude at the high-frequency range.
In summary, the results demonstrated by PFC-MPCD method as shown in figures 5(a)-(b) and 6(a) qualitatively agree with the SLS model as indicated in figures 2(a)-(c). While, all those results exhibited by PFC-PCD method as illustrated in figures 5(a)-(b) and 6(b) tend to differ from the SLS model, in particular, at the high-frequency range as indicated in figures 2(a)-(c).

Prediction of complex moduli parameters
In this part, we aim to demonstrate the difference in complex moduli parameters prediction between PFC-MPCD method and PFC-PCD method. All complex moduli parameters exhibited by the PFC-MPC method are illustrated in figure 7(a). while the results obtained from the PFC-PCD method are showed in figure 7(b) for comparison. According to the numerical study, the storage modulus (E′) predicted by PFC-MPCD tends to monotonically increased with increasing frequency and becomes constant at the certain frequency while loss  modulus (E″) and loss tangent factor d tan ( ( )) tends to increase in magnitude until their peak locations are reached. This is in good agreement with the SLS model prediction in figure 2(d) and consistent with real solid behavior [21]. In particular, all results predicted by PFC-MPCD method are agreeing with numerous experimental observations which can be found in [25,26] for shear and young modulus, which are proportional to storage modulus, in In-Sn and Pb-Sn as well as for loss angle in some alloys [27] which is related to loss tangent factor. Besides, this consistency is even more pronounced in rubber [28][29][30] and in elastomer [31]. This is in contrast to the results obtained from PFC-PCD prediction in figure 7(b) where all complex moduli parameters monotonically increased with increasing frequency without peaks location for loss modulus (E″) and loss tangent factor d tan ( ( )) as well as without the certain frequency value which allows storage modulus E′ to become constant.
To further investigate into the difference between the results predicted by PFC-MPCD and PFC-PCD methods, we compare each parameter exhibited by both methods in figures 8(a)-(c). It is obvious that the results obtained from both methods are in good agreement in a low-frequency range but tend to differ at very frequency. All complex moduli parameters obtained from PFC-PCD method prediction have increased as increasing frequency without the limit. This is unlike the results obtained from PFC-MPCD method prediction. Although all complex moduli predicted by the method increases as increasing frequency, they tend to become limited at a certain frequency value. In particular, all loss parameters, loss modulus (E″) and loss tangent factor d tan ( ( )), predicted by the PFC-MPCD method can demonstrate at least one peak location. These observations are in good agreement with the real solid behavior suggested by [21] and experimental evidence that can be found in [25][26][27][28][29][30][31]. The peak location in loss modulus (E″) parameter also indicates the information on the suitable working frequency where the energy is most dissipated in damping material design. On the contrary, the  PFC-PCD method cannot provide this information as indicated in figure 8(b). It should be noted that the particular form of the PFC free energy equation (1) is applicable to crystalline materials such as metals [34]; therefore, it might seem that the results from this work only restricts to these types of materials. However, since the viscoelastic behavior originates from the MPCD equation, not the PFC free energy, the alternative choices of PFC free energy equation (1) can be used and the PFC-MPCD method should still exhibit similar viscoelastic behavior correspond to those PFC free energy choices. The capability to generalize the PFC free energy enables us to model viscoelastic behavior in different types of materials such as polymeric materials and elastomer where viscoelasticity behavior is much more pronounced than that in metals.

Conclusion
In this work, we propose a new modified pressure-controlled dynamics (MPCD) equation incorporating with the phase-field-crystal (PFC) method which refers to as PFC-MPCD method. This PFC-MPCD method differs from the original PFC-PCD method in that the pressure time derivative term has been introduced to the new version which enables the method to predict mechanical-hysteresis behavior and complex moduli parameters correctly, in particular, at very frequency range where the original PFC-PCD method is incapable. This modification is inspired by the standard linear solid (SLS) model; generally used to predict various viscoelastic phenomena in a solid structure such as viscoelastic-creep, stress-relaxation and mechanical-hysteresis behavior; where the pressure time derivative term is available in its formulation. This time derivative term has played a key role in producing mechanical-hysteresis behavior correctly. The numerical study regarding PFC-MPCD method demonstrates a good agreement between the results predicted by PFC-MPCD method and SLS model in many aspects. In the system response aspect, the PFC-MPCD method shows that the magnitude of strain amplitude reduces with increased sinusoidal pressure oscillation frequency but then becomes finite value at very frequency. This is in contrast to the PFC-PCD method where the magnitude of strain amplitude vanishes at very frequency range. Furthermore, the magnitude of phase angel predicted by PFC-MPCD method is consistent with SLS model unlike PFC-PCD method exhibits extremely large magnitude of phase angle at the highfrequency range without the peak location like SLS and PFC-MPCD method prediction. In mechanicalhysteresis behavior aspect, the prediction of mechanical-hysteresis behavior obtained from both methods is in good agreement at early frequency. However, the mechanical-hysteresis curve predicted by PFC-PCD method tends to be vertical with y-axis which indicates that the system produces zero strain at a very frequency which is inconsistent with SLS and PFC-MPCD method which produce a limited strain value at a very frequency while all loss parameters predicted by the method pass through at least one peak location. These results generally agree with the real solid material behavior and experimental observations which is in contrast to PFC-PCD method prediction at the high-frequency range. In summary, the PFC-PCD method can generally be employed to predict mechanical-hysteresis behavior and complex moduli parameters, however, the application should be limited at the low-frequency range where the viscous effect is not dominated. Whereas, the PFC-MPCD method is more suitable to use in predicting mechanical-hysteresis behavior and complex moduli parameters at the broader pressure oscillation frequency. This is consistent with a real solid material behavior which allows us to predict a suitable working condition for damping materials. We expect that our method can be extended to a more complex system with microstructural features such as grain boundary, dislocation, surface or even with other free energy functional form where viscoelasticity behavior is much more pronounced and still be capable to exhibit the accurate mechanical-hysteresis behavior and complex moduli parameters.