Investigation of Viscoelastic-Creep and Mechanical-Hysteresis Behaviors of Hydrostatically Stressed Crystal Using the Phase Field Crystal Method

Em-Udom, J.; Pisutha-Arnond, N.

doi:https://doi.org/10.1155/2020/2821402

Advances in Mathematical Physics

On this page

Abstract Introduction Background Results Conclusion Appendix Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2020 | Article ID 2821402 | https://doi.org/10.1155/2020/2821402

Investigation of Viscoelastic-Creep and Mechanical-Hysteresis Behaviors of Hydrostatically Stressed Crystal Using the Phase Field Crystal Method

J. Em-Udom¹and N. Pisutha-Arnond¹

Academic Editor: Ivan Giorgio

Received24 Jul 2019

Accepted31 Oct 2019

Published01 Feb 2020

Abstract

The phase field crystal (PFC) method is a density-functional-type model with atomistic resolution and operating on diffusive time scales which has been proven to be an efficient tool for predicting numerous material phenomena. In this work, we first propose a method to predict viscoelastic-creep and mechanical-hysteresis behaviors in a body-centered-cubic (BCC) solid using a PFC method that is incorporated with a pressure-controlled dynamic equation which enables convenient control of deformation by specifying external pressure. To achieve our objective, we use constant pressure for the viscoelastic-creep study and sinusoidal pressure oscillation for the mechanical-hysteresis study. The parametric studies show that the relaxation time in the viscoelastic-creep phenomena is proportional to temperature. Also, mechanical-hysteresis behavior and the complex moduli predicted by the model are consistent with those of the standard linear solid model in a low-frequency pressure oscillation. Moreover, the impact of temperature on complex moduli is also investigated within the solid-stabilizing range. These results qualitatively agree with experimental and theoretical observations reported in the previous literature. We believe that our work should contribute to extending the capability of the PFC method to investigate the deformation problem when the externally applied pressure is required.

1. Introduction

The phase field crystal (PFC) method has emerged as a computational model with atomistic resolution and diffusive time scale. The method has an advantage over molecular dynamics (MD) in terms of the time scale that is not restricted by lattice vibration time; this is due to the specification of the order parameter as the local-time-average atomic number density and the evolution of the order parameter through dissipative dynamics. The method also has advantage over the phase field (PF) method in terms of atomic resolution and self-consistency; this benefit stems from the free energy functional that can be minimized by the periodic order-parameter field whereas the free energy functional of the PF method is formulated by a spatially uniform field which diminishes numerous physical features that occurs due to the periodicity of crystalline phases, e.g., spatial symmetry, elastic and plastic interaction, nucleation, multiple crystal orientations, and motion of dislocations [1–3]. The PFC method was first introduced by Elder et al. [1, 2] and has been extensively used to study material phenomena and behavior ranging from phase transformation [4, 5] and topological defect dynamics (vacancy [6], grain boundary [7–9], dislocation [3, 10, 11], and crack [12]) to elasticity and elastic constants [13–15] and plasticity [11, 16]. Moreover, the PFC method was also interpreted and derived according to the classical density functional theory (CDFT) of the freezing point of view [4]. This derivation provided an additional field variable which extended the capability of the PFC method to investigate the material phenomena in a more complex situation such as phase transformation [4, 17] and segregation-induced grain-boundary premelting [18] in binary alloys.

Although there were numerous successful PFC studies on the mechanical behavior of materials, such as elasticity, plasticity, and phase transformation in pure materials and binary alloys, the investigation of viscoelastic behavior and the related properties predicted by the PFC model is still limited and not thoroughly explored. Thus, it is of interest to study the capability of the PFC model to predict these phenomena. Viscoelasticity is the behavior with both viscous and elastic characteristics; thus, materials with viscoelastic behavior can return to their original configuration when unloaded but do so in a time-dependent manner. The viscous property also leads to the fact that material response depends on the rate at which it is deformed [19, 20]. Viscoelastic effect has played significant roles in many materials ranging from amorphous polymers [21], semicrystalline polymers [22], biopolymers [23], and bitumen materials [24] to metals at very high temperature [25]. Furthermore, viscoelasticity governs many engineering applications in various fields including damping material design for noise reduction [26, 27] and vibration control in engineering structure [28], shape-memory ceramics and polymers [29], piezoelectric materials [30], self-healing materials [31], viscoelastic gels in surgery application [32], low-loss materials [33], and nanoscale resonators [34, 35].

In quantifying the viscoelastic effect, there are three types of studies: viscoelastic creep, stress relaxation, and mechanical hysteresis. In this work, we not only focus on the viscoelastic creep but also extend our study to mechanical-hysteresis behavior and its related parameters in a body-centered-cubic (BCC) solid (3D case) predicted by the PFC method, unlike our previous works [36, 37] which were solely limited to viscoelastic creep in a one-dimensional (1D case) crystal study; moreover, the mechanical-hysteresis behavior and its related parameters were not addressed in those works. Since the studies of viscoelastic creep and mechanical hysteresis involve measuring the strain response from an applied stress, we also need a method to control the applied stress, and thus, we employ the PFC model that is incorporated with our previously proposed pressure controlled dynamic (PCD) equation [37]. The use of the PCD equation enables convenient control of deformation using external pressure which allows studies of experimentally observed phenomena involving a system under different types of applied pressure; we will refer to this particular PFC method as the PFC-PCD model. The PCD equation is similar to the dynamic equation introduced by Kocher and Provatas [38]. They used their PCD equation to control internal pressure by specifying external pressure to investigate the full spectrum of solid-liquid-vapor transitions; unfortunately, the viscoelastic-creep and mechanical-hysteresis behaviors which can be produced from that equation were not the focus area and, of course, not addressed in their study [38]. The difference between two PCD versions is due to the different derivation process in which our version is established regarding the thermodynamics of a hydrostatically stressed crystal solid proposed by Larché and Cahn [39] and classical irreversible thermodynamic frameworks [40–42]. Moreover, the deviation of viscoelastic-creep behavior in 1D crystal between two PCD versions was also investigated in our previous study [37]. The advantage of the PCD equation is that it enables the external pressure to be specified as an input to the simulation, which allows the system volume, or the strain response, to be measured; this is in contrast to the conventional PFC simulations where volume, or grid spacing, is an independent variable [43].

For the investigation of viscoelastic behavior, we simulate a hydrostatic stress via a BCC solid subjected to two types of applied pressure: constant pressure for the viscoelastic-creep study and sinusoidal pressure oscillation for the mechanical-hysteresis study. We find that, in both cases, the solid exhibits delayed strain responses from the applied pressure which is indicative of viscoelasticity, and the degree of viscoelasticity can be controlled by the parameter in the PCD equation. Furthermore, the strain responses from the simulations can be compared with functional forms of the solutions to the standard linear solid model to extract viscoelastic quantities and properties. These properties include relaxation time of viscoelastic creep and the complex moduli from the hysteresis. The parametric studies show that the relaxation time is proportional to temperature while the magnitude of deformation at steady-state condition is proportional to both temperature and atomic density. Also, the dynamic behavior, mechanical hysteresis, and complex moduli predicted by the PFC-PCD equation under sinusoidal pressure oscillation are consistent with those by a standard linear solid model at a certain frequency range. However, at high frequency range, the behavior exhibited by the PFC-PCD model differs from that by the standard linear solid model. The reason behind this difference is thoroughly analyzed by conducting Taylor’s expansion on PCD equation in Appendix C. Moreover, the impact of temperature on complex moduli is also investigated, and the results well agree with experimental and theoretical observation. We are positive that our work should provide the extension of capability for the PFC method to investigate the deformation problem in a real crystalline solid when externally applied pressure is necessary.

This paper is organized as follows. In Section II, we provide the background information on the PFC method, the pressure-control dynamic equation, and the standard linear solid model. In Section III, the method of simulation is presented. The results and discussions are provided in Section IV, and we conclude the work with the summary in Section V.

2. Background

2.1. PFC Method

The PFC method is characterized by its free energy functional and the dynamic equation. The simplest free energy functional was introduced by Elder and coworkers [1, 2]: is the free energy density; , , , and are fitting parameters; and is the atomic number density field. This PFC free energy functional can be minimized by a constant profile, representing a liquid phase, and a periodic profile, representing a crystalline solid phase. The density periodic field can be expressed in terms of Fourier expansion of the form where is a reciprocal lattice vector (RLV), is an amplitude of the density wave corresponding to , is an average atomic number density, and c.c. denotes a complex conjugate. The minimization of essentially results in the truncation of the summation in Equation (2) to terms with small values of as the terms with higher values of have vanishingly small . Consequently, this particular PFC model (Equation (1)) favors crystal structures that can be constructed from a single set of RLVs such as stripe, hexagonal, and BCC structures in one, two, and three dimensions, respectively. To this end, one can construct the so-called “one-mode” approximation of the BCC crystal by limiting the terms in Equation (2) to those corresponding to RLVs or , where is the lattice parameter. By assuming that of the terms with the same are equivalent, one obtains where and is the density-wave amplitude.

To simplify the expressions in Equation (1), one can nondimensionalize the variables by using the following substitutions [2]: where the quantities with tildes are nondimensional and is the dimensionality of the problem. The expressions in Equation (1) simplify to where can be interpreted as the degree of undercooling which is inversely proportional to temperature. The one-mode approximation of the BCC crystal becomes where . The value of can be determined from minimization of the free energy density [44] leading to and, consequently, .

One can also approximate the ranges of and where the BCC solid is stable. This is achieved by calculating the free energy density of the solids using the one-mode approximations for stripe, 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 where the BCC solid is stable can be obtained [45].

The evolution equation for is the Cahn-Hilliard type which involves dissipative dynamics and mass conservation [2]. where is the mobility coefficient and 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 ), 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 [43]. In the context of the viscoelastic behavior, this technique would be appropriate for the stress-relaxation calculation; however, for the viscoelastic-creep and mechanical-hysteresis studies, where the system is subjected to specified external stress, an additional dynamic equation is needed.

2.2. Pressure-Controlled Dynamic Equation

The study of viscoelastic creep and mechanical hysteresis involves specifying the external stress, or, in this study, external pressure, and investigating how the (average) internal pressure changes temporally; this requires the dynamic equation that contains external pressure as an independent variable. Such equation, referred to as the pressure-controlled dynamic (PCD) equation, was first introduced by Kocher and Provatas [38] expressed as where is the externally applied pressure and is the mobility coefficient. This equation can be used in conjunction with Equation (7) to simulate the system under specified . If the external pressure is higher (lower) than the internal pressure, or , the value of will be negative (positive) which indicates that the system will undergo compressive (tensile) deformation.

A slightly modified version of the PCD was recently proposed by Em-Udom and Pisutha-Arnond [37]. The equation is expressed as where is the mobility coefficient. This equation was developed using the principle of thermodynamics of a hydrostatically stressed solid [39, 46] and classical irreversible thermodynamic framework [40–42]. The comparison between results between Equations (9) and (10) shows that the steady-state behavior obtained from both equations is identical; however, the nonequilibrium behaviors are different. This difference is due to the extra volume term () on the right-hand side of Equation (10). This extra volume term results in a higher deformation rate in tension (due to increasing ) than the deformation rate in compression (due to decreasing ). Nevertheless, for small deformation, the numerical difference from Equations (9) and (10) is unlikely to change the qualitative interpretation of the results. Despite the fact that we employ Equation (9) in this work, the choice of the PCD equation is immaterial for this study. Hereafter, we will refer to the PFC method that incorporates the PCD equation as the PFC-PCD model. Lastly, the tilde notation will also be omitted for simplicity. More detail on the PFC-PCD model is provided in Appendix A.

2.3. Standard Linear Solid Model

We review the standard linear solid (SLS) model which is a method of modeling viscoelasticity using linear combinations of springs and dampers as shown in Figure 1. The strain response from this model will be compared with that from the PFC-PCD simulation in order to extract viscoelastic quantities and properties. The governing differential equation for the SLS model is [20, 47] where is the strain response from the applied stress . Equation (11) can be written in simple form as follows: where and . The variable is the scaled quantity defined as .

For the creep behavior, we consider the case where the system (initially at zero stress and strain) is subject to constant stress:

The corresponding strain is then where is the deformation at and is the relaxation time; the subscript “” indicates the creep phenomena. The quantity is the deformation at steady state while relates to the rate at which the system reaches the steady-state deformation. The plot of the strain response is shown as an example in Figure 2. For the hysteresis behavior, we consider the system induced by harmonic excitation given by where is the stress amplitude, is the angular frequency, and is the frequency. At steady state, the strain response yields the following solution (see Appendix B):

where is the strain amplitude. As a result, the amplitude ratio, , is

Also, the loss tangent factor, defines the phase angle due to time-delayed response. From the solution, the complex moduli (or dynamic modulus) can be defined from (all derivation details are given in Appendix B)

The quantity is a property of viscoelastic materials and can be written further as , where is the storage modulus and is the loss modulus. The quantity measures the stored energy and is an elastic response of material while measures energy dissipation and is a viscous response of a material. The expressions for and can be obtained from rearranging the last equality in Equation (20), yielding

The loss tangent, , is the ratio of storage and loss modulus or ; this quantity is a measure of damping in the material.

The behavior of the SLS model, a function of excitation frequency, is shown in Figure 3(a). At low frequencies, or is low, which indicates that the system behaves in an elastic manner. This is due to fact that the dashpot in Figure 1 has sufficient time to displace, resulting in the stress and strain profiles that are almost in phase. The value of is relatively small compared with that at higher frequencies due to large strain amplitude (Figure 3(b)). At high frequencies, the system also behaves elastically as seen from small ; at this condition, the change in stress is too rapid for the dashpot to operate, and thus, the strain response is in phase with the stress profile. Due to small strain amplitude, the value of is relatively large.

(a)

(b)

(c)

(d)

Figure 3

Hysteresis behavior regarding the SLS model. (a) Complex moduli, storage modulus (solid line on -axis), loss modulus (dashed line on -axis) and loss tangent factor (dense dashed line on -axis) established by the SLS model at conditions , , and and frequency, (-axis), ranges between 0.01 and 100: box area depicts Figure 12(b). (b) Amplitude ratio, (-axis), at frequency, (-axis), ranges between 0.004 and 10. (c) System response schematic demonstrates the phase lag between the stress or pressure input at time , (solid line), and the strain response at time , (dashed line). (d) Mechanical-hysteresis schematic plot; strain response at time , (-axis), vs. stress or pressure input at time , (-axis).

At intermediate frequencies, a considerable amount of phase lag occurs, as shown schematically in Figure 3(c), showing a more viscous behavior, and the stress-strain profiles can be plotted to exhibit a hysteresis loop as shown in Figure 3(d) where the area inside the loop represents the energy loss. At the frequency where is maximum, the system exhibits the highest damping capacity, and in real materials, this condition is important for damping applications.

3. Method of Simulation

In this section, we describe the setup and method of simulations. In all simulations, we construct a BCC unit cell under hydrostatic pressure (Figure 4) with a periodic grid where and the lattice parameter ; this results in the reference grid spacing, in all dimensions, and the reference volume . The initial density profile is constructed with the one-mode approximation in Equation (6), where the reference density, , is specified. The profile is then relaxed using the evolution equation (Equation (7)), while fixing and , until the equilibrium density profile is reached. At this state, the reference internal pressure, , can be calculated using Equation (8). Also, the numerical method employed to solve the evolution equation, and also the PCD equation (Equation (10)), is the Fourier pseudospectral method where the Fourier transform is used to calculate spatial derivatives and the forward finite difference scheme is used to discretize the time derivative (see Appendix A). For the viscoelastic simulations, the PCD equation is employed where different profiles of are specified and the change in the volume is used to update the grid spacing through the relation

We note that the grid spacing is identical in all directions due to the assumption of a 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 equation. The system response is characterized by the strain defined by which is identical in all directions. The profiles of and are then used to quantify the viscoelastic behavior predicted form the PFC-PCD method.

Two types of are used in this study. For the viscoelastic creep simulation, the applied pressure is set to where is a constant. The values of are chosen such that the pressure deviation from the reference pressure is positive or negative; the former leads to tensile deformation while the latter results in compressive deformation. The quantities and are numerically extracted from the strain response, , where is the values of at large and is the value of at (see Equation (14) and Figure 2).

For the 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-state strain response is reached. The strain response profile at steady state will exhibit the sinusoidal behavior similar to Figure 3(a), and we can numerically measure the strain amplitude in order to calculate the amplitude ratio . When plotting the hysteresis loop as in Figure 3(d), the complex moduli and loss tangent can be calculated by numerically identifying the values of (pressure value at zero strain), (pressure value at maximum strain), and (maximum strain value). Then, we evaluate the expressions , , and . We note that in all hysteresis simulations, is set to .

Finally, the values of and in all simulations are restricted to the range of (−0.248, −0.260) and (0.175, 0.200), respectively. These values of and ensure the stability of the bcc crystal as shown in the phase diagram in Ref. [43].

4. Results and Discussions

4.1. Viscoelastic-Creep Behavior

In this part, we show the viscoelastic-creep behavior exhibited by the PFC-PCD model and demonstrate the dependence of this behavior on and . In Figure 5, we aim to show the strain response from different values of . The graphs show the strain, initially at zero, increase or decrease gradually until the steady-state values (denoted as ) are reached. These delayed responses obtained from the simulations are indicative of the viscoelastic-creep behavior; this means that the PFC model combined with the PCD equation is capable of modeling viscoelastic behavior. Also, depending on the values of , either tensile or compressive deformations are realized. When , or , tensile deformation is obtained while , or , results in compressive deformation. The magnitudes of (steady-state deformation) also scale with the magnitudes of , as expected from the elastic effect in the PFC model (Equation (1)).

Next, Figure 6 shows how viscoelastic behavior is affected by . From the figure, the increase in results in approaching at a faster rate, or smaller ; this means that the increase in leads to the faster strain response or higher degree of elasticity. This result demonstrates that can be used to adjust the degree of viscoelasticity ranging from highly viscous to highly elastic. The capability to control the degree of viscoelastic creep can be used to tune the PFC-PCD model to different types of materials.

The degree of viscoelasticity is also affected by ϵ (inverse of temperature), as shown in Figure 7. From the figure, the values of decreases as the values of ϵ increases for the range of where the solid is stable. Since ϵ is inversely proportional to temperature, this indicates that the system responds faster, or becomes more elastic, as temperature decreases; this prediction agrees with the experimental observation where relaxation time in creep phenomena decreases with temperature.

We note that the results from Figures 6 and 7 assume that and are independent from one another. Nevertheless, since both and affect the viscoelasticity behavior, it is likely that and (or temperature) are related. From the derivation in [37], originates from linear phenomenological law in the classical irreversible thermodynamic framework. Therefore, the dependence of on temperature could be introduced by considering a more rigorous treatment of this phenomenological law.

Lastly, we report how the stiffness changes with and . This is shown in Figures 8(a) and 8(b), where the values of are plot as functions of ϵ and , respectively. Since is the same in all simulations, lower values indicate higher stiffness and vice versa. From Figure 8(a), decreases as ϵ increases while, from Figure 8(b), increases as increases. This indicates that the stiffness is higher at lower temperature and lower atomic densities. These results are qualitatively consistent with the theoretical framework [13, 48] and experimental observation [49–51]. We note that stiffness is not a viscoelastic property and the consistency of the stiffness variation with and originates from the PFC free energy (Equation (1)), not the PCD equation. Nevertheless, these stiffness results demonstrate the application of the PCD equation since it allows convenient calculation of from specified values of .

(a)

(b)

4.2. Hysteresis Behavior

In this part, we report the hysteresis behavior exhibited by the PFC-PCD model and demonstrate the dependency of this behavior on (pressure-oscillation frequency), (PCD parameter), and (inverse of temperature). The steady-state strain response from the PFC-PCD model is shown in Figures 9(a)–9(c) while the strain response from the SLS model is shown in Figures 9(d)–9(f) for comparison. For the and , the PFC-PCD model shows the phase lag between the strain and pressure profiles and also the increase in the phase lag with increasing frequencies; these results qualitatively agree with the result from the SLS model. However, the phase lag from the PFC-PCD model does not decrease at high frequency as in the SLS model as shown for ; the strain response from the PFC-PCD model still shows significant phase lag, while for the SLS model, both the stress and strain profiles are in phase. The results indicate that the PFC-PCD model does not exhibit elastic behavior at very frequency. This suggests that even though the PFC-PCD model is capable of exhibiting viscoelastic response (or phase lag) from the oscillating pressure, the applicability of these phenomena is still limited to the low-frequency excitation.

(a)

(b)

(c)

(d)

(e)

(f)

Another comparison between the PFC-PCD and the SLS model can be shown in Figure 10 where the amplitude ratios are shown. Both the PFC-PCD and the SLS model show reducing amplitudes with increasing frequencies and both amplitude ratios exhibit qualitatively similar functional forms. However, the amplitude ratio from the PFC-PCD model is vanishingly small at high frequencies while the amplitude ratio from the SLS model remains at finite values. The values of the amplitude ratio can also be seen from the hysteresis loop in Figure 11 where the results in Figure 9 are plotted on stress/pressure-strain axes. As the frequency increases, the hysteresis loop rotates counter-clockwise, indicating smaller strain amplitude (or amplitude ratio). At , the hysteresis loop from the PFC-PCD model almost forms a vertical line due to an essentially zero amplitude ratio while the hysteresis loop from the SLS model remains tilted away from the -axis. This extremely small amplitude ratio indicates that the system produces almost zero strain from the input pressure, which leads to unnaturally large values of the moduli ( and ). Therefore, the application of the PFC-PCD method to model dynamic viscoelastic behavior should be limited to low-frequency excitation. Nevertheless, at the low frequency, the prediction of , , and is in good agreement with the results from the SLS model as shown in Figure 12.

(a)

(b)

(a)

(b)

(a)

(b)

To further investigate into the difference between the results from the PFC-PCD and those from the SLS model, we expand the PCD equation in terms of strain and perform Taylor’s expansion to the first order with respect to the strain; the details of the derivation are shown in Appendix C. Compared with the SLS equation, the transformed PCD equation lacks the term with the time derivative of stress. This analysis indicates that the PCD equation is more similar to the Kelvin-Voigt model than the SLS model. The Kelvin-Voigt model consists of a spring and a dashpot in parallel and is capable of reproducing the creep phenomena, but not the stress relaxation. This similarity between the PCD and the Kelvin-Voigt models explains why the PFC-PCD model can reproduce creep phenomena, but not the dynamic behavior similar to the SLS model. We note that the lack of similarity to the SLS model not only applies to our recently proposed PCD equation [36] but also applies to the PCD equation proposed by Kocher and Provatas [38] since the functional forms of both equations are the same. Nevertheless, this analysis suggests the possibility of improving the PFC-PCD model by modifying the PCD equation in such a way that the effect of the term with the time derivative of stress is present.

Similar to the result from the viscoelastic-creep study, the PCD parameter can be used to control the degree of viscoelasticity. Figures 13(a) and 13(b) show the strain response with different values of at a fixed pressure-oscillation frequency. The hysteresis loops from these results are also shown in Figure 14. The results show that decreasing the value of leads to higher phase lag (Figure 13) and higher damping capacity (Figure 14). This indicates that lowering the value leads to the system having a more viscous response. The opposite is also true where increasing leads to a more elastic response from the system.

(a)

(b)

The influence of on , , and is shown in Figure 15. Figure 15(a) shows that an increase in (decrease in temperature) leads to an increase in and a decrease in . As , the values of decrease with increasing , as shown in Figure 15(b). The results show that as temperature decreases, the system behaves in a less viscous manner and with reducing damping capacity. These trends qualitatively agree with polymeric materials below glass transition [52] and metals such as aluminum alloy [53].

(a)

(b)

It should be noted that the particular form of the PFC free energy (Equation (1)) is applicable to crystalline materials such as metals [45]; therefore, it might seem that the results from this work only pertain to these types of materials. Nevertheless, since the viscoelastic behavior originates from the PCD equation, not the PFC free energy, an alternative form of PFC free energy (Equation (1)) can be used and the modified PFC-PCD model should still exhibit similar viscoelastic behavior as well as similar dependence of viscoelastic properties on the PFC-PCD model parameters. The capability to generalize the PFC free energy allows modeling of viscoelastic behavior in different types of materials such as polymeric materials where viscoelasticity is much more pronounced than that in metals.

5. Conclusion

In this study, we aim to show the capability and the limitation of the PFC-PCD model which was first employed to investigate the viscoelastic-creep and mechanical-hysteresis behaviors in a BCC (3D case) unit cell under hydrostatic pressure. To achieve our goal, we implement two types of pressure profiles in our analysis: constant pressure input for viscoelastic-creep study and sinusoidal pressure oscillation for mechanical-hysteresis study where all PFC parameters, temperature and atomic density parameters, are established in a solid stable region. The conclusions for each analysis are summarized as below: 1.In the viscoelastic-creep study, we first investigate the influence of two PCD parameters: magnitude of pressure input and mobility coefficient on viscoelastic-creep behavior. Regarding the results, we find the magnitude of pressure input scales with that of deformation at steady-state time both tensile and compressive loads. The mobility coefficient can be used to control the degree of viscosity of deformation in that higher mobility results in higher elasticity as well as lower mobility leads to higher viscosity. This evidence suggests that we can control both the magnitude of deformation and the degree of viscoelasticity through PCD parameters. Second, we study the impact of PFC parameters: temperature on relaxation time and the magnitude of deformation at steady-state condition as well as atomic density parameter on the magnitude of deformation at steady-state condition. Regarding the results, we find that the system response is faster and has less relaxation time, less deformation, and more stiffness, at lower temperature which is generally consistent with theoretical framework and experimental observation. Also, the system exhibits larger deformation and less stiffness, at higher atomic density that shows an agreement with previous work2.In dynamic behavior study, we first focus on studying the impact of excitation frequency on strain response, amplitude ratio, mechanical hysteresis, and complex moduli. The results predicted by the PFC-PCD model indicate that the system displays more phase lag between sinusoidal pressure oscillation and strain response with reducing amplitude ratio at increased frequency which qualitatively agrees with that of the SLS model. However, at high frequency range, the results exhibited by both models are inconsistent. This finding can also be observed in the mechanical-hysteresis loop. The hysteresis loop exhibited by PFC-PCD becomes a vertical line which indicates the system produce almost zero strain. This finding produces unusually large values of complex moduli unlike the SLS model prediction. Furthermore, we show that the mobility coefficient can also be used to control the degree of viscoelasticity in dynamic behavior under sinusoidal pressure oscillation. The impact of temperature on complex moduli predicted by the PFC-PCD model is also in good agreement with experimental observation

Finally, we show that the PCD model is more similar to the Kelvin-Voigt model than the SLS model regarding Taylor’s expansion in Appendix C. This deviation does not seem to cause any difference in viscoelastic-creep behavior, but the mechanical-hysteresis behavior tends to differ from the SLS model at high frequency range. However, the prediction of the model agrees well with the SLS model at the low-frequency range. Therefore, the application of the PFC-PCD method to predict dynamic viscoelastic behavior should be limited in some certain frequency region. This analysis provides the possibility of enhancing the PFC-PCD model by modifying the PCD equation in a more proper way which might be addressed in future improvement.

Appendix

A. PFC-PCD Equations and Their Numerical Scheme

The PFC-PCD equation consists of two governing equations: Cahn-Hilliard- (CH-) type equation and the pressure-controlled dynamic (PCD) equation. The first equation is involved in dissipative dynamics and a mass conservation equation [2] formulated in partial differential equation (PDE) with four independent variables ; the time variable with three spatial coordinate variables where is a density field variable in three-dimensional space. is the mobility coefficient. is a three-dimensional Laplacian operator.

is the variational derivative of by which defines a chemical potential . The free energy functional can be expressed or in alternative form where .

The free energy density variable, , can be expressed in expansion form as

Denote that where each independent variable of is defined by

To solve Equation (A.1), the variable must be reformulated as a function of and its derivative. Therefore, we need to find by using the functional derivative definition

Evaluate each derivative term on the right-hand side of Equation (A.6):

Substitute each term of Equation (A.7) into Equation (A.6). Then, the chemical potential is literally obtained

Finally, Equation (A.1) becomes where . It is obvious that Equation (A.9) is the sixth-order nonlinear partial differential equation (6OD-NPDE) which can be expressed in a full description as

In order to solve Equation (A.10) numerically, there are two options which can be employed: finite difference method (FD) and Fourier spectral method (FSM). In this study, not only is the periodic boundary condition governed but also the equation itself has a very high derivative order in spatial coordinates which demands a lot of numerical expense if the FD scheme is used. Therefore, the FSM method is more appropriate than FD in terms of numerical accuracy and numerical expense to solve Equation (A.10). To employ the FSM scheme, we begin with Equation (A.9). If we apply the Fourier transform both sides, we obtain

Then, the 6OD-NPDE is literally reduced to an empirical ordinary differential equation which depends on time constraint only, where is the Fourier transform of . is the Fourier transform of which can be described as follows:

is a vector in Fourier space, and each Fourier wave vector component is defined by where is the lattice parameter along the direction and our study uses for . is the number of grid points along the direction in which we employ for . It should be noted that

We employ a semi-implicit scheme [54] to solve Equation (A.11). The forward difference approximation for the time derivative term is given as where and denote at time and , respectively. Substitute Equation (A.15) in Equation (A.11) where and are at time and at time , respectively. Rearranging the above equation, we obtain and it yields

Finally, the variable which is in reciprocal space can be converted into the real space variable by using inverse Fourier transform where is dimensionality of the problem. This variable will be used to calculate in the pressure-controlled dynamic equation further. Next, the second important equation is the PCD equation which was developed regarding the thermodynamics of the hydrostatically stressed crystal solid [39] and classical irreversible thermodynamic frameworks [40–42]. The advantage of this equation is that it allows creating the deformation in the PFC model by specifying the external pressure as an input variable unlike the conventional PFC simulations where the volume, or grid spacing, is an input variable [43]. This equation is given as follows:

To solve the PCD equation, we use the forward difference scheme for the time derivative term.

Then, Equation (A.20) becomes where and are the volume at time and , respectively. Also, and are and at time , respectively. Then, Equation (A.22) yields

Finally, we literally obtain a new volume value at the next time step as we required. It should be noted that is where

All superscripts, , regarding Equation (A.25) which appear in all represent those physical quantities, , at time .

B. Standard Linear Solid Model

In this appendix, we aim to derive the expressions for the complex modulus and strain response from the SLS model under harmonic excitation. First, we derive complex modulus under dynamic excitation. Consider a more general form of the stress (scaled quantity): where is the amplitude, is the angular frequency and is the complex amplitude defined by . At steady state, the strain response will also exhibit a harmonic function with the same frequency; therefore, one can postulate the functional form of the strain as an input function [20]: where is the strain amplitude. Substituting Equations (B.1) and (B.2) into Equation (12), we obtain

This equation defines where each of and is and , respectively. Therefore, we literally obtain

To derive the solution expression, we consider stress as an input function and propose the strain function be written as

Substituting Equations (B.5) and (B.6) into Equation (12) to obtain

For a particular case where , therefore, the strain response corresponding to this case can be obtained by that yields where , which can be reduced to where

It should be noted that .

C. Dynamic Similarity between PCD Equation and SLS Model

In this appendix, we demonstrate how dynamic similarity between the PCD Equation (10) and the SLS model Equation (12) can be established using Taylor’s expansion [55] at infinitesimal deformation around the reference state. First, we begin with PCD equation where and represents the average internal pressure and external applied pressure with where is the crystal length parameter. If we denote that and represents a current volume. It should be noted that where all are spatial coordinates in the Eulerian description. The integration boundary for each coordinate is to where and . Equation (C.2) becomes

If we assume infinitesimal deformation, small strain can be written as where and . Therefore, Equation (C.3) can be written in expansion form [55]: where represents at and represents a higher order derivative term as a function of . In case of infinitesimal deformation, the higher order term, where , can be truncated, then Equation (C.4) is reduced to where . Also, can be written as where represents average reference pressure calculated under an undeformed state and represents the deviation part of from . Since is a composite function, therefore, the time derivative regarding Equation (C.1) can be expressed by the chain rule where and . Substitute Equation (C.5)–(C.8) into Equation (C.1), we obtain that yields

Note that . This equation can be reduced to

Regarding binomial expansion, we know that for the first term approximation if . Equation (C.11) becomes

If we neglect and and use the fact that , Equation (C.12) turns into Then, Equation (C.13) finally yields the simplified model of PCD equation as which obviously has structure like a Kelvin-Voigt model: where and . Equation (C.15) is slightly different from the SLS model Equation (12) in that the term is introduced on the right-hand side in the SLS model. This deviation does not cause the difference in viscoelastic-creep behavior produced by both the PCD model and the SLS model solution under constant pressure. However, in case of the mechanical hysteresis, although this behavior produced by the PCD equation is deviated with that of the SLS model and behaves like the Kelvin-Voigt model at an elevated frequency range, it consistently predicts similar results with the SLS model at low frequency. This can be described by the solution obtained from both models under the sinusoidal load, . According to Equation (B.9), the strain solution produced by the SLS model is . However, if , the solution will converge to which is almost identical to that of a simplified model of the PCD equation.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the National Science and Technology Development Agency (NSTDA)-University-Industry Research Collaboration (NUI-RC) program under Grant No. NUI-RC-M33-22-59-001D; the Thailand Research Fund (TRF) under Grant No. TRG5880008; the New Researcher Scholarship of Coordinating Center for Thai Government Science and Technology Scholarship Students, Ministry of Science and Technology Grant No. FDA-CO-2561-840g-TH; King Mongkut’s Institute of Technology Ladkrabang Research Fund under Grant No. KREF046103; and the High Performance Computing Services from (Thailand) National Electronics and Computer Technology Center (NECTEC). The authors are also indebted to Dr. Sasawat Mahabunphachai for his helpful discussion and encouragement on the research.

References

K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, “Modeling elasticity in crystal growth,” Physical Review Letters, vol. 88, no. 24, article 245701, 2002.
View at: Publisher Site | Google Scholar
K. R. Elder and M. Grant, “Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals,” Physical Review E, vol. 70, no. 5, article 051605, 2004.
View at: Publisher Site | Google Scholar
J. Berry, M. Grant, and K. R. Elder, “Diffusive atomistic dynamics of edge dislocations in two dimensions,” Physical Review E, vol. 73, no. 3, article 031609, 2006.
View at: Publisher Site | Google Scholar
K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, “Phase-field crystal modeling and classical density functional theory of freezing,” Physical Review B, vol. 75, no. 6, article 064107, 2007.
View at: Publisher Site | Google Scholar
T. Pusztai, G. Tegze, G. I. Tóth et al., “Phase-field approach to polycrystalline solidification including heterogeneous and homogeneous nucleation,” Journal of Physics: Condensed Matter, vol. 20, no. 40, article 404205, 2008.
View at: Publisher Site | Google Scholar
P. Y. Chan, N. Goldenfeld, and J. Dantzig, “Molecular dynamics on diffusive time scales from the phase-field-crystal equation,” Physical Review E, vol. 79, no. 3, article 035701, 2009.
View at: Publisher Site | Google Scholar
T. Hirouchi, T. Takaki, and Y. Tomita, “Effects of temperature and grain size on phase-field-crystal deformation simulation,” International Journal of Mechanical Sciences, vol. 52, no. 2, pp. 309–319, 2010.
View at: Publisher Site | Google Scholar
J. Berry, K. R. Elder, and M. Grant, “Melting at dislocations and grain boundaries: a phase field crystal study,” Physical Review B, vol. 77, no. 22, article 224114, 2008.
View at: Publisher Site | Google Scholar
J. Mellenthin, A. Karma, and M. Plapp, “Phase-field crystal study of grain-boundary premelting,” Physical Review B, vol. 78, no. 18, article 184110, 2008.
View at: Publisher Site | Google Scholar
P. Stefanovic, M. Haataja, and N. Provatas, “Phase-field crystals with elastic interactions,” Physical Review Letters, vol. 96, no. 22, article 225504, 2006.
View at: Publisher Site | Google Scholar
P. Y. Chan, G. Tsekenis, J. Dantzig, K. A. Dahmen, and N. Goldenfeld, “Plasticity and dislocation dynamics in a phase field crystal model,” Physical Review Letters, vol. 105, no. 1, article 015502, 2010.
View at: Publisher Site | Google Scholar
S. Hu and S. Wang, “The influences of crystal orientation and crack interaction on the initiation of growth and propagation mode of microcrack: a phase-field-crystal study,” Physica B: Condensed Matter, vol. 552, pp. 104–109, 2019.
View at: Publisher Site | Google Scholar
N. Pisutha-Arnond, V. Chan, K. Elder, and K. Thornton, “Calculations of isothermal elastic constants in the phase-field crystal model,” Physical Review B, vol. 87, no. 1, article 014103, 2013.
View at: Publisher Site | Google Scholar
Z.-L. Wang, Z.-F. Huang, and Z. Liu, “Elastic constants of stressed and unstressed materials in the phase-field crystal model,” Physical Review B, vol. 97, no. 14, article 144112, 2018.
View at: Publisher Site | Google Scholar
W. Zhou, J. Wang, Z. Wang et al., “Elastic strain response in the modified phase-field-crystal model,” Chinese Physics B, vol. 26, no. 9, article 090702, 2017.
View at: Publisher Site | Google Scholar
P. Stefanovic, M. Haataja, and N. Provatas, “Phase field crystal study of deformation and plasticity in nanocrystalline materials,” Physical Review E, vol. 80, no. 4, article 046107, 2009.
View at: Publisher Site | Google Scholar
M. Greenwood, N. Ofori-Opoku, J. Rottler, and N. Provatas, “Modeling structural transformations in binary alloys with phase field crystals,” Physical Review B, vol. 84, no. 6, article 064104, 2011.
View at: Publisher Site | Google Scholar
Y.-L. Lu, T.-T. Hu, G.-M. Lu, and Z. Chen, “Phase-field crystal study of segregation induced grain-boundary premelting in binary alloys,” Physica B: Condensed Matter, vol. 451, pp. 128–133, 2014.
View at: Publisher Site | Google Scholar
L. S. Roderic, Viscoelastic Solids, CRC press, Boca Raton, FL, USA, 1998.
W. N. Findley and F. A. Davis, Creep and Relaxation of Nonlinear Viscoelastic Materials, Courier Corporation, North Chelmsford, MA, USA, 2013.
J. D. Ferry, Viscoelastic Properties of Polymers, John Wiley & Sons, Hoboken, NJ, USA, 1980.
A. D. Drozdov and J. C. de Christiansen, “Viscoelasticity and viscoplasticity of semicrystalline polymers: structure-property relations for high-density polyethylene,” Computational materials science, vol. 39, no. 4, pp. 729–751, 2007.
View at: Publisher Site | Google Scholar
E. Frey, K. Kroy, and J. Wilhelm, Viscoelasticity of Biopolymer Networks and Statistical Mechanics of Semiflexible Polymers, Elsevier, JAI, New York City, NY, USA, 1999.
H. Di Benedetto, F. Olard, C. Sauzéat, and B. Delaporte, “Linear viscoelastic behaviour of bituminous materials: from binders to mixes,” Road Materials and Pavement Design, vol. 5, supplement 1, pp. 163–202, 2004.
View at: Publisher Site | Google Scholar
I. Ritchie, K. Sprungmann, and M. Sahoo, “Internal friction in sonoston - a high damping Mn/Cu-based alloy for marine propeller applications,” Le Journal de Physique Colloques, vol. 46, no. C10, pp. C10–409, 1985.
View at: Publisher Site | Google Scholar
J. Kim and J.-Y. Choi, “Performance test for transmitted noise reduction of smart panels using piezoelectric shunt damping,” Smart Materials and Structures, vol. 14, no. 4, pp. 587–593, 2005.
View at: Publisher Site | Google Scholar
R. Fan, G. Meng, J. Yang, and C. He, “Experimental study of the effect of viscoelastic damping materials on noise and vibration reduction within railway vehicles,” Journal of Sound and Vibration, vol. 319, no. 1-2, pp. 58–76, 2009.
View at: Publisher Site | Google Scholar
S. Elias and V. Matsagar, “Research developments in vibration control of structures using passive tuned mass dampers,” Annual Reviews in Control, vol. 44, pp. 129–156, 2017.
View at: Publisher Site | Google Scholar
K. Otsuka and C. M. Wayman, Shape Memory Materials, Cambridge university press, Cambridge, UK, 1999.
C. Rosen, B. V. Hiremath, and R. Newnham, Piezoelectricity, Springer Science & Business Media, Berlin, Germany, 1992.
S. R. White, N. Sottos, P. Geubelle et al., “Autonomic healing of polymer composites,” Nature, vol. 409, no. 6822, pp. 794–797, 2001.
View at: Publisher Site | Google Scholar
R. O. St, C. Weiss, J. L. Denlinger, and E. A. Balazs, “A preliminary assessment of Na-hyaluronate injection into “no man’s land” for primary flexor tendon repair,” Clinical Orthopaedics and Related Research, vol. &NA;, no. 146, p. 269???275, 1980.
View at: Publisher Site | Google Scholar
V. B. Braginsky, V. Mitrofanov, and V. I. Panov, Systems with Small Dissipation, University of Chicago Press, Chicago, IL, USA, 1985.
H. G. Craighead, “Nanoelectromechanical Systems,” Science, vol. 290, no. 5496, pp. 1532–1535, 2000.
View at: Publisher Site | Google Scholar
R. Lifshitz and M. L. Roukes, “Thermoelastic damping in micro- and nanomechanical systems,” Physical Review B, vol. 61, no. 8, pp. 5600–5609, 2000.
View at: Publisher Site | Google Scholar
E. Jatupon and P. Nirand, “Investigation on viscoelastic-creep behavior of the phase-field crystal method,” IOP conference series: materials science and engineering, vol. 361, article 012009, 2018.
View at: Publisher Site | Google Scholar
J. Em-Udom and N. Pisutha-Arnond, “Pressure-controlled dynamic equation for the phase-field crystal method,” IOP conference series: materials science and engineering, vol. 361, article 012006, 2018.
View at: Publisher Site | Google Scholar
G. Kocher and N. Provatas, “New density functional approach for solid-liquid-vapor transitions in pure materials,” Physical Review Letters, vol. 114, no. 15, article 155501, 2015.
View at: Publisher Site | Google Scholar
F. Larché and J. W. Cahn, “Eine Lineare theorie des thermodynamischen Gleichgewichts von festkorpern unter spannung,” Acta metallurgica, vol. 21, no. 8, pp. 1051–1063, 1973.
View at: Publisher Site | Google Scholar
S. R. De Groot and P. Mazur, Non-Equilibrium Thermodynamics, Courier Corporation, North Chelmsford, MA, USA, 2013.
Y. Mishin, J. Warren, R. Sekerka, and W. Boettinger, “Irreversible thermodynamics of creep in crystalline solids,” Physical Review B, vol. 88, no. 18, article 184303, 2013.
View at: Publisher Site | Google Scholar
G. Lebon, D. Jou, and J. Casas-Vázquez, Understanding Non-Equilibrium Thermodynamics, Springer, Berlin, Germany, 2008.
T. Hirouchi, T. Takaki, and Y. Tomita, “Development of numerical scheme for phase field crystal deformation simulation,” Computational Materials Science, vol. 44, no. 4, pp. 1192–1197, 2009.
View at: Publisher Site | Google Scholar
N. Provatas and K. Elder, Phase-Field Methods in Materials Science and Engineering, John Wiley & Sons, Hoboken, NJ, USA, 2011.
K.-A. Wu and A. Karma, “Phase-field crystal modeling of equilibrium bcc-liquid interfaces,” Physical Review B, vol. 76, no. 18, article 184107, 2007.
View at: Publisher Site | Google Scholar
P. Voorhees and W. C. Johnson, Solid State Physics-Advances in Research and Applications, vol. 59, Academic Press Inc, Cambridge, MA, USA, 2004.
R. S. Lakes, Viscoelastic Materials, Cambridge University Press, Cambridge, UK, 2009.
J. Wachtman Jr., W. Tefft, D. Lam Jr., and C. Apstein, “Exponential temperature dependence of Young's modulus for several oxides,” Physical review, vol. 122, no. 6, pp. 1754–1759, 1961.
View at: Publisher Site | Google Scholar
P. M. Sutton, “The variation of the elastic constants of crystalline aluminum with temperature between 63°K and 773°K,” Physical Review, vol. 91, no. 4, pp. 816–821, 1953.
View at: Publisher Site | Google Scholar
W. Chancellor, A. Wolfenden, and G. Ludtka, “Temperature dependence of Young's modulus and shear modulus in uranium-2.4 wt% niobium alloy,” Journal of Nuclear Materials, vol. 171, no. 2-3, pp. 389–394, 1990.
View at: Publisher Site | Google Scholar
R. Bruls, H. Hintzen, G. De With, and R. Metselaar, “The temperature dependence of the Young's modulus of MgSiN₂, AlN and Si₃N₄,” Journal of the European Ceramic Society, vol. 21, no. 3, pp. 263–268, 2001.
View at: Publisher Site | Google Scholar
D. I. Jones, Handbook of Viscoelastic Vibration Damping, John Wiley & Sons, Hoboken, NJ, USA, 2001.
L. Licitra, D. D. Luong, O. M. Strbik III, and N. Gupta, “Dynamic properties of alumina hollow particle filled aluminum alloy A356 matrix syntactic foams,” Materials & Design, vol. 66, pp. 504–515, 2015.
View at: Publisher Site | Google Scholar
L. Q. Chen and J. Shen, “Applications of semi-implicit Fourier-spectral method to phase field equations,” Computer Physics Communications, vol. 108, no. 2-3, pp. 147–158, 1998.
View at: Publisher Site | Google Scholar
G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists, AAPT, Baltimore, MA, USA, 1999.

Copyright

Copyright © 2020 J. Em-Udom and N. Pisutha-Arnond. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

716

Downloads

862

Citations