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

Authors: J. Em-Udom; N. Pisutha-Arnond
Journal: Advances in Mathematical Physics
 (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.


---

## Body

## 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 AppendixC. 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]:
(1)ℱ=∫wρdr,wρ≡ρ2at+λq02+∇22ρ+gtρ44,wρ is the free energy density; at, gt, λ, and q0 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
(2)ϕr,ρ¯=ρ¯+∑jAjeiGj⋅r+c.c.,where Gj is a reciprocal lattice vector (RLV), Aj is an amplitude of the density wave corresponding to Gj, ρ¯ is an average atomic number density, and c.c. denotes a complex conjugate. The minimization of F essentially results in the truncation of the summation in Equation (2) to terms with small values of Gj as the terms with higher values of Gj have vanishingly small Aj. 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 110 RLVs or Gj=2π√2/La, where La is the lattice parameter. By assuming that Aj of the terms with the same Gj are equivalent, one obtains
(3)ρone=ρ¯+4Ascosqar1cosqar2+cosqar1cosqar3+cosqar2cosqar3,where qa=2π/La and As is the density-wave amplitude.To simplify the expressions in Equation (1), one can nondimensionalize the variables by using the following substitutions [2]:
(4)r~≡q0r,∈~≡−atλq04,ρ~≡gtλq04ρ,ℱ~≡gtλ2q08−dℱ,w~≡gtλ2q08w,where the quantities with tildes are nondimensional and d=3 is the dimensionality of the problem. The expressions in Equation (1) simplify to
(5)ℱ~=∫w~ρ~dr~,w~ρ~=ρ~2−∈~+1+∇~22ρ~+ρ~44,where ϵ can be interpreted as the degree of undercooling which is inversely proportional to temperature. The one-mode approximation of the BCC crystal becomes
(6)ρ~one=ρ¯~+4A~scosq~ar~1cosq~ar~2+cosq~ar~1cosq~ar~3+cosq~ar~2cosq~ar~3,where q~a=qa/q0=2π/Laq0=2π/L~a. The value of q~a can be determined from minimization of the free energy density [44] leading to q~a=1/2 and, consequently, L~a=2π2.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].
(7)∂ρ~∂t~=L~μ∇~2μ~,μ~=1−∈~ρ~+ρ~3+2∇~2ρ~+∇~4ρ~,where L~μ 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(8)P¯~int=1V~∫−f~+μ~ρ~dr~,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
(9)dV~dt~=L~def′P¯~int−P¯~ext,where P¯~ext is the externally applied pressure and L~def′ is the mobility coefficient. This equation can be used in conjunction with Equation (7) to simulate the system under specified P¯~ext. If the external pressure is higher (lower) than the internal pressure, or P¯~ext>P¯~intP¯~ext<P¯~int, the value of dV~/dt~ 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
(10)dV~dt~=L~defV~P¯~int−P¯~ext,where L~def 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 (V~) on the right-hand side of Equation (10). This extra volume term results in a higher deformation rate in tension (due to increasing V~) than the deformation rate in compression (due to decreasing V~). 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 Figure1. 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]
(11)E1+E2dedt+E1E2ηe=E1ησ+dσdt,where e is the strain response from the applied stress σ. Equation (11) can be written in simple form as follows:
(12)dedt+Ke=σ′+τdσ′dt,where K=E1E2/ηE1+E2 and τ=η/E1. The variable σ′ is the scaled quantity defined as σ′=σ/τE1+E2.Figure 1
Schematic illustration of standard linear solid (SLS) model.For the creep behavior, we consider the case where the system (initially at zero stress and strain) is subject to constant stress:(13)σ′t=σ¯′.The corresponding strain is then(14)et=e∞1−exp−tτc,where e∞=σ¯′/K is the deformation at t=∞ and τc=1/K is the relaxation time; the subscript “c” indicates the creep phenomena. The quantity e∞ is the deformation at steady state while τc 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
(15)σ′t=σA′sinωt,where σA′ is the stress amplitude, ω=2πf is the angular frequency, and f is the frequency. At steady state, the strain response yields the following solution (see Appendix B):
(16)et=eAsinωt−δ,Figure 2
Viscoelastic creep under step load at the samee∞=2·10−3 condition and τc=9.0 where t (x-axis) is the time variable and et (y-axis) is the strain response at time t.where(17)eA=σA′1+τω2K2+ω2,is the strain amplitude. As a result, the amplitude ratio, Ar, is
(18)Ar=σA′σA1+τω2K2+ω2.Also, the loss tangent factor,(19)δ=arctan1−τKK+τω2ω,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)
(20)Gcomplex′≡σ′eA=K+jω1+jτω=Gstorage′+jGloss′.The quantityGcomplex′ is a property of viscoelastic materials and can be written further as Gcomplex′=Gstorage′+jGloss′, where Gstorage′ is the storage modulus and Gloss′ is the loss modulus. The quantity Gstorage′ measures the stored energy and is an elastic response of material while Gloss′ measures energy dissipation and is a viscous response of a material. The expressions for Gstorage′ and Gloss′ can be obtained from rearranging the last equality in Equation (20), yielding
(21)Gstorage′=K+τω21+τω2,Gloss′=1−τKω1+τω2⋅The loss tangent,tanδ, is the ratio of storage and loss modulus or tanδ=Gloss′/Gstorage′; this quantity is a measure of damping in the material.The behavior of the SLS model, a function of excitation frequency, is shown in Figure3(a). At low frequencies, tanδ 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 Gstorage′ 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 tanδ; 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 Gstorage′ is relatively large.Figure 3
Hysteresis behavior regarding the SLS model. (a) Complex moduli, storage modulusGstorage′ (solid line on y-axis), loss modulus Gloss′ (dashed line on y-axis) and loss tangent factor tanδ (dense dashed line on y-axis) established by the SLS model at conditions σA′=0.00075, K=0.3750, and τ=0.1667 and frequency, f (x-axis), ranges between 0.01 and 100: box area depicts Figure 12(b). (b) Amplitude ratio, Ar (y-axis), at frequency, f (x-axis), ranges between 0.004 and 10. (c) System response schematic demonstrates the phase lag between the stress or pressure input at time t, σ′t (solid line), and the strain response at time t, et (dashed line). (d) Mechanical-hysteresis schematic plot; strain response at time t, et (x-axis), vs. stress or pressure input at time t, σ′t (y-axis).
(a)
(b)
(c)
(d)At intermediate frequencies, a considerable amount of phase lag occurs, as shown schematically in Figure3(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 tanδ is maximum, the system exhibits the highest damping capacity, and in real materials, this condition is important for damping applications.


## 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]:
(1)ℱ=∫wρdr,wρ≡ρ2at+λq02+∇22ρ+gtρ44,wρ is the free energy density; at, gt, λ, and q0 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
(2)ϕr,ρ¯=ρ¯+∑jAjeiGj⋅r+c.c.,where Gj is a reciprocal lattice vector (RLV), Aj is an amplitude of the density wave corresponding to Gj, ρ¯ is an average atomic number density, and c.c. denotes a complex conjugate. The minimization of F essentially results in the truncation of the summation in Equation (2) to terms with small values of Gj as the terms with higher values of Gj have vanishingly small Aj. 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 110 RLVs or Gj=2π√2/La, where La is the lattice parameter. By assuming that Aj of the terms with the same Gj are equivalent, one obtains
(3)ρone=ρ¯+4Ascosqar1cosqar2+cosqar1cosqar3+cosqar2cosqar3,where qa=2π/La and As is the density-wave amplitude.To simplify the expressions in Equation (1), one can nondimensionalize the variables by using the following substitutions [2]:
(4)r~≡q0r,∈~≡−atλq04,ρ~≡gtλq04ρ,ℱ~≡gtλ2q08−dℱ,w~≡gtλ2q08w,where the quantities with tildes are nondimensional and d=3 is the dimensionality of the problem. The expressions in Equation (1) simplify to
(5)ℱ~=∫w~ρ~dr~,w~ρ~=ρ~2−∈~+1+∇~22ρ~+ρ~44,where ϵ can be interpreted as the degree of undercooling which is inversely proportional to temperature. The one-mode approximation of the BCC crystal becomes
(6)ρ~one=ρ¯~+4A~scosq~ar~1cosq~ar~2+cosq~ar~1cosq~ar~3+cosq~ar~2cosq~ar~3,where q~a=qa/q0=2π/Laq0=2π/L~a. The value of q~a can be determined from minimization of the free energy density [44] leading to q~a=1/2 and, consequently, L~a=2π2.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].
(7)∂ρ~∂t~=L~μ∇~2μ~,μ~=1−∈~ρ~+ρ~3+2∇~2ρ~+∇~4ρ~,where L~μ 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(8)P¯~int=1V~∫−f~+μ~ρ~dr~,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
(9)dV~dt~=L~def′P¯~int−P¯~ext,where P¯~ext is the externally applied pressure and L~def′ is the mobility coefficient. This equation can be used in conjunction with Equation (7) to simulate the system under specified P¯~ext. If the external pressure is higher (lower) than the internal pressure, or P¯~ext>P¯~intP¯~ext<P¯~int, the value of dV~/dt~ 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
(10)dV~dt~=L~defV~P¯~int−P¯~ext,where L~def 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 (V~) on the right-hand side of Equation (10). This extra volume term results in a higher deformation rate in tension (due to increasing V~) than the deformation rate in compression (due to decreasing V~). 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 Figure1. 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]
(11)E1+E2dedt+E1E2ηe=E1ησ+dσdt,where e is the strain response from the applied stress σ. Equation (11) can be written in simple form as follows:
(12)dedt+Ke=σ′+τdσ′dt,where K=E1E2/ηE1+E2 and τ=η/E1. The variable σ′ is the scaled quantity defined as σ′=σ/τE1+E2.Figure 1
Schematic illustration of standard linear solid (SLS) model.For the creep behavior, we consider the case where the system (initially at zero stress and strain) is subject to constant stress:(13)σ′t=σ¯′.The corresponding strain is then(14)et=e∞1−exp−tτc,where e∞=σ¯′/K is the deformation at t=∞ and τc=1/K is the relaxation time; the subscript “c” indicates the creep phenomena. The quantity e∞ is the deformation at steady state while τc 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
(15)σ′t=σA′sinωt,where σA′ is the stress amplitude, ω=2πf is the angular frequency, and f is the frequency. At steady state, the strain response yields the following solution (see Appendix B):
(16)et=eAsinωt−δ,Figure 2
Viscoelastic creep under step load at the samee∞=2·10−3 condition and τc=9.0 where t (x-axis) is the time variable and et (y-axis) is the strain response at time t.where(17)eA=σA′1+τω2K2+ω2,is the strain amplitude. As a result, the amplitude ratio, Ar, is
(18)Ar=σA′σA1+τω2K2+ω2.Also, the loss tangent factor,(19)δ=arctan1−τKK+τω2ω,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)
(20)Gcomplex′≡σ′eA=K+jω1+jτω=Gstorage′+jGloss′.The quantityGcomplex′ is a property of viscoelastic materials and can be written further as Gcomplex′=Gstorage′+jGloss′, where Gstorage′ is the storage modulus and Gloss′ is the loss modulus. The quantity Gstorage′ measures the stored energy and is an elastic response of material while Gloss′ measures energy dissipation and is a viscous response of a material. The expressions for Gstorage′ and Gloss′ can be obtained from rearranging the last equality in Equation (20), yielding
(21)Gstorage′=K+τω21+τω2,Gloss′=1−τKω1+τω2⋅The loss tangent,tanδ, is the ratio of storage and loss modulus or tanδ=Gloss′/Gstorage′; this quantity is a measure of damping in the material.The behavior of the SLS model, a function of excitation frequency, is shown in Figure3(a). At low frequencies, tanδ 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 Gstorage′ 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 tanδ; 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 Gstorage′ is relatively large.Figure 3
Hysteresis behavior regarding the SLS model. (a) Complex moduli, storage modulusGstorage′ (solid line on y-axis), loss modulus Gloss′ (dashed line on y-axis) and loss tangent factor tanδ (dense dashed line on y-axis) established by the SLS model at conditions σA′=0.00075, K=0.3750, and τ=0.1667 and frequency, f (x-axis), ranges between 0.01 and 100: box area depicts Figure 12(b). (b) Amplitude ratio, Ar (y-axis), at frequency, f (x-axis), ranges between 0.004 and 10. (c) System response schematic demonstrates the phase lag between the stress or pressure input at time t, σ′t (solid line), and the strain response at time t, et (dashed line). (d) Mechanical-hysteresis schematic plot; strain response at time t, et (x-axis), vs. stress or pressure input at time t, σ′t (y-axis).
(a)
(b)
(c)
(d)At intermediate frequencies, a considerable amount of phase lag occurs, as shown schematically in Figure3(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 tanδ 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 (Figure4) with a periodic N×N×N grid where N=16 and the lattice parameter La=2π2; this results in the reference grid spacing, Δr0=2π/8 in all dimensions, and the reference volume V0=2π23. The initial density profile is constructed with the one-mode approximation in Equation (6), where the reference density, ρ~0, is specified. The profile is then relaxed using the evolution equation (Equation (7)), while fixing Δr=Δr0 and ρ¯=ρ¯0, until the equilibrium density profile is reached. At this state, the reference internal pressure, P¯0, 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 P¯extt are specified and the change in the volume is used to update the grid spacing through the relation
(22)Δrt=VtN31/3.Figure 4
BCC unit cell under hydrostatic pressure schematic using isosurface construction;σii=P¯ext depicts stress in the ii direction where all σij components, shear stress components i≠j, are equal to zero.We note that the grid spacingΔrt 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
(23)ρ¯t=ρ¯0V0Vt.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(24)et=Δrt−Δr0Δr0,which is identical in all directions. The profiles of et and P¯extt are then used to quantify the viscoelastic behavior predicted form the PFC-PCD method.Two types ofP¯extt are used in this study. For the viscoelastic creep simulation, the applied pressure is set to
(25)P¯extt=P¯ext,where P¯ext is a constant. The values of P¯ext are chosen such that the pressure deviation from the reference pressure
(26)ΔP¯0≡P¯0−P¯extP¯0,is positive or negative; the former leads to tensile deformation while the latter results in compressive deformation. The quantities e∞ and τc are numerically extracted from the strain response, et, where e∞ is the values of et at large t and τc is the value of t at eτc=e∞1−exp−1≈0.6321e∞ (see Equation (14) and Figure 2).For the hysteresis simulation, the applied pressure is set to(27)P¯extt=P¯0−P¯Asinωt,where P¯A is the pressure amplitude. The pressure deviation from the reference pressure will then be
(28)ΔP¯0t≡P¯0−P¯exttP¯0=P¯AP¯0sinωt.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 Figure3(a), and we can numerically measure the strain amplitude eA in order to calculate the amplitude ratio Ar=eA/P¯A. 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 A (pressure value at zero strain), B (pressure value at maximum strain), and C (maximum strain value). Then, we evaluate the expressions G′storage=B/C, G′loss=A/C, and tanδ=G′loss/G′storage. We note that in all hysteresis simulations, P¯A is set to 0.02P¯0.Finally, the values ofP¯0 and ϵ in all simulations are restricted to the range of (−0.248, −0.260) and (0.175, 0.200), respectively. These values of P¯0 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ΔP¯0 and Ldef. In Figure 5, we aim to show the strain response from different values of ΔP¯0. The graphs show the strain, initially at zero, increase or decrease gradually until the steady-state values (denoted as e∞) 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 ΔP¯0, either tensile or compressive deformations are realized. When ΔP¯0>0, or P¯ext<P¯0, tensile deformation is obtained while ΔP¯0<0, or P¯ext>P¯0, results in compressive deformation. The magnitudes of e∞ (steady-state deformation) also scale with the magnitudes of ΔP¯0, as expected from the elastic effect in the PFC model (Equation (1)).Figure 5
Viscoelastic-creep behaviors,et, under different values of constant pressure, ΔP¯0, at conditions ρ¯0=0.2600, ϵ=0.2000, and Ldef=1 where t (x-axis) is the time variable and et (y-axis) is the strain response at time t. The results show the capability of the PCD equation to control e∞ through ΔP¯0.Next, Figure6 shows how viscoelastic behavior is affected by Ldef. From the figure, the increase in Ldef results in et approaching e∞ at a faster rate, or smaller τc; this means that the increase in Ldef leads to the faster strain response or higher degree of elasticity. This result demonstrates that Ldef 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.Figure 6
Viscoelastic-creep behaviors,et, under different values of Ldef at conditions ρ¯0=0.2600, ϵ=0.2000, and ΔP¯0=0.03 where t (x-axis) is the time variable and et (y-axis) is the strain response at time t. The results demonstrate the capability of the PCD equation to control the relaxation time, τc, through Ldef.The degree of viscoelasticity is also affected byϵ (inverse of temperature), as shown in Figure 7. From the figure, the values of τc decreases as the values of ϵ increases for the range of P¯0 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.Figure 7
Viscoelastic-creep parametric study under step pressure; influence of temperature,ϵ (x-axis); parameter on relaxation time, τc (y-axis); parameter at ϵ parameter ranges between 0.175 and 0.200 with three different atomic density, ρ¯0, values. All numerical experiments are based on step pressure input (ΔP¯0=0.02). It should be noted that the ϵ parameter is inversely related to temperature.We note that the results from Figures6 and 7 assume that Ldef and ϵ are independent from one another. Nevertheless, since both Ldef and ϵ affect the viscoelasticity behavior, it is likely that Ldef and ϵ (or temperature) are related. From the derivation in [37], Ldef originates from linear phenomenological law in the classical irreversible thermodynamic framework. Therefore, the dependence of Ldef on temperature could be introduced by considering a more rigorous treatment of this phenomenological law.Lastly, we report how the stiffness changes withϵ and ρ¯0. This is shown in Figures 8(a) and 8(b), where the values of e∞ are plot as functions of ϵ and ρ¯0, respectively. Since ΔP¯0 is the same in all simulations, lower e∞ values indicate higher stiffness and vice versa. From Figure 8(a), e∞ decreases as ϵ increases while, from Figure 8(b), e∞ increases as ρ¯0 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 ρ¯0 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 e∞ from specified values of ΔP¯0.Figure 8
Viscoelastic-creep parametric study under step pressure. (a) Influence of the temperature,ϵ (x-axis), parameter on e∞ (y-axis) at the ϵ parameter ranges between 0.175 and 0.200 with three different ρ¯0 values. (b) Influence of the atomic density, ρ¯0 (x-axis), parameter on e∞ (y-axis) at the ρ¯0 parameter ranges between 0.248 and 0.260 with three different ϵ values. All numerical experiments are based on step pressure input (ΔP¯0=0.02). It should be noted that the ϵ parameter is inversely related to temperature.
(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 onf (pressure-oscillation frequency), Ldef (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 f=0.004 and f=0.0232, 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 f=100; 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.Figure 9
Influence of frequency,f, on system response; sinusoidal pressure input ΔP¯0t (solid line on y-axis) and strain response et (dashed line on y-axis) established by the PCD equation: (a) f=0.004, (b) f=0.0232, and (c) f=100 at conditions ρ¯0=0.2600, ϵ=0.1750, and Ldef=1; sinusoidal input σ′t (solid line on y-axis) and strain response et (dashed line on y-axis) established by the SLS model regarding equation (16): (d) f=0.004, (e) f=0.0232, and (f) f=100 at conditions σA′=0.00075, K=0.3750, and τ=1.6670.
(a)
(b)
(c)
(d)
(e)
(f)Another comparison between the PFC-PCD and the SLS model can be shown in Figure10 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 f=100, 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 y-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 (G′storage and G′loss). 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 G′storage, G′loss, and tanδ is in good agreement with the results from the SLS model as shown in Figure 12.Figure 10
Influence of frequency,f (x-axis), on the amplitude ratio, Ar (y-axis), parameter under sinusoidal input (a) PCD equation and (b) SLS model, at conditions ρ¯0=0.2600 and ϵ=0.1750.
(a)
(b)Figure 11
Influence of frequency,f on mechanical-hysteresis behavior. (a) Mechanical hysteresis established by the PCD equation at conditions ρ¯0=0.2600 and ϵ=0.1750 and (b) mechanical hysteresis established by the SLS model at conditions σA′=0.00075, K=0.3750, and τ=1.6670 where et (x-axis) is the strain response at time t whereas ΔP¯0t (y-axis) and σ′t (y-axis) are the pressure input at time t established by the PCD equation and SLS model, respectively.
(a)
(b)Figure 12
Influence of frequency,f (x-axis), on complex moduli (y-axis); G′storage (storage modulus), G′loss (loss modulus), and tanδ (loss tangent factor); (a) complex moduli (y-axis) established by the PCD equation at conditions ρ¯0=0.2600 and ϵ=0.1750; (b) complex moduli (y-axis) established by the SLS model at conditions σA′=0.00075, K=0.3750, and τ=0.1667 and f ranges between 0.01 and 0.18.
(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 AppendixC. 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 parameterLdef can be used to control the degree of viscoelasticity. Figures 13(a) and 13(b) show the strain response with different values of Ldef 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 Ldef leads to higher phase lag (Figure 13) and higher damping capacity (Figure 14). This indicates that lowering the Ldef value leads to the system having a more viscous response. The opposite is also true where increasing Ldef leads to a more elastic response from the system.Figure 13
Influence of theLdef parameter on system response; sinusoidal pressure input ΔP¯0t (solid line on y-axis) and strain response et (dashed line on y-axis) at time t (x-axis); (a) Ldef=10 and (b) Ldef=1, at ρ¯0=0.2600, ϵ=0.1750, and f=0.02.
(a)
(b)Figure 14
Influence ofLdef on mechanical-hysteresis behavior at conditions ρ¯0=0.2600, ϵ=0.1750, and f=0.02 where et (x-axis) is the strain response at time t and ΔP¯0t (y-axis) is the pressure input at time t established by the PCD equation.The influence ofϵ on G′storage, G′loss, and tanδ is shown in Figure 15. Figure 15(a) shows that an increase in ϵ (decrease in temperature) leads to an increase in G′storage and a decrease in G′loss. As tanδ=G′loss/G′storage, the values of tanδ 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].Figure 15
Influence of temperature,ϵ (x-axis), parameter on (a) storage modulus, G′storage (solid line on y-axis), and loss modulus, G′loss (dashed line on y-axis), and (b) loss tangent factor, tanδ (y-axis), at ρ¯0=0.2600 and f=0.02. It should be noted that the ϵ parameter is inversely related to temperature.
(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.


## 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ΔP¯0 and Ldef. In Figure 5, we aim to show the strain response from different values of ΔP¯0. The graphs show the strain, initially at zero, increase or decrease gradually until the steady-state values (denoted as e∞) 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 ΔP¯0, either tensile or compressive deformations are realized. When ΔP¯0>0, or P¯ext<P¯0, tensile deformation is obtained while ΔP¯0<0, or P¯ext>P¯0, results in compressive deformation. The magnitudes of e∞ (steady-state deformation) also scale with the magnitudes of ΔP¯0, as expected from the elastic effect in the PFC model (Equation (1)).Figure 5
Viscoelastic-creep behaviors,et, under different values of constant pressure, ΔP¯0, at conditions ρ¯0=0.2600, ϵ=0.2000, and Ldef=1 where t (x-axis) is the time variable and et (y-axis) is the strain response at time t. The results show the capability of the PCD equation to control e∞ through ΔP¯0.Next, Figure6 shows how viscoelastic behavior is affected by Ldef. From the figure, the increase in Ldef results in et approaching e∞ at a faster rate, or smaller τc; this means that the increase in Ldef leads to the faster strain response or higher degree of elasticity. This result demonstrates that Ldef 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.Figure 6
Viscoelastic-creep behaviors,et, under different values of Ldef at conditions ρ¯0=0.2600, ϵ=0.2000, and ΔP¯0=0.03 where t (x-axis) is the time variable and et (y-axis) is the strain response at time t. The results demonstrate the capability of the PCD equation to control the relaxation time, τc, through Ldef.The degree of viscoelasticity is also affected byϵ (inverse of temperature), as shown in Figure 7. From the figure, the values of τc decreases as the values of ϵ increases for the range of P¯0 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.Figure 7
Viscoelastic-creep parametric study under step pressure; influence of temperature,ϵ (x-axis); parameter on relaxation time, τc (y-axis); parameter at ϵ parameter ranges between 0.175 and 0.200 with three different atomic density, ρ¯0, values. All numerical experiments are based on step pressure input (ΔP¯0=0.02). It should be noted that the ϵ parameter is inversely related to temperature.We note that the results from Figures6 and 7 assume that Ldef and ϵ are independent from one another. Nevertheless, since both Ldef and ϵ affect the viscoelasticity behavior, it is likely that Ldef and ϵ (or temperature) are related. From the derivation in [37], Ldef originates from linear phenomenological law in the classical irreversible thermodynamic framework. Therefore, the dependence of Ldef on temperature could be introduced by considering a more rigorous treatment of this phenomenological law.Lastly, we report how the stiffness changes withϵ and ρ¯0. This is shown in Figures 8(a) and 8(b), where the values of e∞ are plot as functions of ϵ and ρ¯0, respectively. Since ΔP¯0 is the same in all simulations, lower e∞ values indicate higher stiffness and vice versa. From Figure 8(a), e∞ decreases as ϵ increases while, from Figure 8(b), e∞ increases as ρ¯0 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 ρ¯0 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 e∞ from specified values of ΔP¯0.Figure 8
Viscoelastic-creep parametric study under step pressure. (a) Influence of the temperature,ϵ (x-axis), parameter on e∞ (y-axis) at the ϵ parameter ranges between 0.175 and 0.200 with three different ρ¯0 values. (b) Influence of the atomic density, ρ¯0 (x-axis), parameter on e∞ (y-axis) at the ρ¯0 parameter ranges between 0.248 and 0.260 with three different ϵ values. All numerical experiments are based on step pressure input (ΔP¯0=0.02). It should be noted that the ϵ parameter is inversely related to temperature.
(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 onf (pressure-oscillation frequency), Ldef (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 f=0.004 and f=0.0232, 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 f=100; 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.Figure 9
Influence of frequency,f, on system response; sinusoidal pressure input ΔP¯0t (solid line on y-axis) and strain response et (dashed line on y-axis) established by the PCD equation: (a) f=0.004, (b) f=0.0232, and (c) f=100 at conditions ρ¯0=0.2600, ϵ=0.1750, and Ldef=1; sinusoidal input σ′t (solid line on y-axis) and strain response et (dashed line on y-axis) established by the SLS model regarding equation (16): (d) f=0.004, (e) f=0.0232, and (f) f=100 at conditions σA′=0.00075, K=0.3750, and τ=1.6670.
(a)
(b)
(c)
(d)
(e)
(f)Another comparison between the PFC-PCD and the SLS model can be shown in Figure10 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 f=100, 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 y-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 (G′storage and G′loss). 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 G′storage, G′loss, and tanδ is in good agreement with the results from the SLS model as shown in Figure 12.Figure 10
Influence of frequency,f (x-axis), on the amplitude ratio, Ar (y-axis), parameter under sinusoidal input (a) PCD equation and (b) SLS model, at conditions ρ¯0=0.2600 and ϵ=0.1750.
(a)
(b)Figure 11
Influence of frequency,f on mechanical-hysteresis behavior. (a) Mechanical hysteresis established by the PCD equation at conditions ρ¯0=0.2600 and ϵ=0.1750 and (b) mechanical hysteresis established by the SLS model at conditions σA′=0.00075, K=0.3750, and τ=1.6670 where et (x-axis) is the strain response at time t whereas ΔP¯0t (y-axis) and σ′t (y-axis) are the pressure input at time t established by the PCD equation and SLS model, respectively.
(a)
(b)Figure 12
Influence of frequency,f (x-axis), on complex moduli (y-axis); G′storage (storage modulus), G′loss (loss modulus), and tanδ (loss tangent factor); (a) complex moduli (y-axis) established by the PCD equation at conditions ρ¯0=0.2600 and ϵ=0.1750; (b) complex moduli (y-axis) established by the SLS model at conditions σA′=0.00075, K=0.3750, and τ=0.1667 and f ranges between 0.01 and 0.18.
(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 AppendixC. 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 parameterLdef can be used to control the degree of viscoelasticity. Figures 13(a) and 13(b) show the strain response with different values of Ldef 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 Ldef leads to higher phase lag (Figure 13) and higher damping capacity (Figure 14). This indicates that lowering the Ldef value leads to the system having a more viscous response. The opposite is also true where increasing Ldef leads to a more elastic response from the system.Figure 13
Influence of theLdef parameter on system response; sinusoidal pressure input ΔP¯0t (solid line on y-axis) and strain response et (dashed line on y-axis) at time t (x-axis); (a) Ldef=10 and (b) Ldef=1, at ρ¯0=0.2600, ϵ=0.1750, and f=0.02.
(a)
(b)Figure 14
Influence ofLdef on mechanical-hysteresis behavior at conditions ρ¯0=0.2600, ϵ=0.1750, and f=0.02 where et (x-axis) is the strain response at time t and ΔP¯0t (y-axis) is the pressure input at time t established by the PCD equation.The influence ofϵ on G′storage, G′loss, and tanδ is shown in Figure 15. Figure 15(a) shows that an increase in ϵ (decrease in temperature) leads to an increase in G′storage and a decrease in G′loss. As tanδ=G′loss/G′storage, the values of tanδ 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].Figure 15
Influence of temperature,ϵ (x-axis), parameter on (a) storage modulus, G′storage (solid line on y-axis), and loss modulus, G′loss (dashed line on y-axis), and (b) loss tangent factor, tanδ (y-axis), at ρ¯0=0.2600 and f=0.02. It should be noted that the ϵ parameter is inversely related to temperature.
(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 work
2.
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 observationFinally, we show that the PCD model is more similar to the Kelvin-Voigt model than the SLS model regarding Taylor’s expansion in AppendixC. 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.



---

*Source: 2821402-2020-02-01.xml*