**2. Materials and Methods**

#### *2.1. Principle of Compressional Viscoelastography* **2. Materials and Methods**

*Materials* **2021**, *14*, x FOR PEER REVIEW 3 of 20

The design of the compressional viscoelastography system is shown in Figure 2 (refer to references [33,34] for a more detailed introduction of the design of this system). The system uses a circular pressure plate to apply a uniform compressional pressure on the surface of the material. The pressure and back plates are connected to each other and attached to the ultrasound transducer. The magnitude of the compressional pressure can be measured by load cells attached between the two plates. The system has a feedback control unit that allows for the simultaneous monitoring and control of the compressional pressure. The movement of the transducer is controlled by a stepper motor. The strain of an element within the material is measured using ultrasound imaging. *2.1. Principle of Compressional Viscoelastography*  The design of the compressional viscoelastography system is shown in Figure 2 (refer to references [33,34] for a more detailed introduction of the design of this system). The system uses a circular pressure plate to apply a uniform compressional pressure on the surface of the material. The pressure and back plates are connected to each other and attached to the ultrasound transducer. The magnitude of the compressional pressure can be measured by load cells attached between the two plates. The system has a feedback control unit that allows for the simultaneous monitoring and control of the compressional

In the measurement, the transducer is moved downward to compress the material at a constant loading rate until the preset pressure level is reached. The loading rate must be high to mimic a step load. Once the preset pressure level is reached, the compressional pressure is maintained as constant through a feedback control system. During the period when the compressional pressure is maintained as constant, the stress of an element within the material is constant as well. Therefore, each element within the material exhibits creep behavior. The measurement processes are illustrated in Figure 3. pressure. The movement of the transducer is controlled by a stepper motor. The strain of an element within the material is measured using ultrasound imaging. In the measurement, the transducer is moved downward to compress the material at a constant loading rate until the preset pressure level is reached. The loading rate must be high to mimic a step load. Once the preset pressure level is reached, the compressional pressure is maintained as constant through a feedback control system. During the period when the compressional pressure is maintained as constant, the stress of an element

In the present study, the Maxwell representation of the standard linear solid model (hereafter referred to as the standard linear solid model for brevity), as shown in Figure 4, is used to describe the viscoelastic behaviors of materials. If a step stress excitation is used to excite the material, the creep behavior (increasing strain over time) described by the standard linear solid model is [35]: within the material is constant as well. Therefore, each element within the material exhibits creep behavior. The measurement processes are illustrated in Figure 3. In the present study, the Maxwell representation of the standard linear solid model (hereafter referred to as the standard linear solid model for brevity), as shown in Figure 4, is used to describe the viscoelastic behaviors of materials. If a step stress excitation is used

$$
\begin{array}{c}
\epsilon(t) = \frac{\sigma\_0}{\sigma\_0 E\_1} \left(1 - \frac{E\_2}{E\_2 + E\_2} \cdot e^{-\frac{1}{\tau\_c}t} \right) \\
\epsilon(t) = \frac{\sigma\_0 E\_1}{E\_1} \left(1 - \frac{E\_2 + E\_2}{E\_1 + E\_2} \cdot e^{-\frac{1}{\tau\_c}t} \right)
\end{array}
\tag{1}
$$

where *ǫ*(*t*) is the strain over time of an element within the material; *τ<sup>C</sup>* = *η*(*E*<sup>1</sup> + *E*2)/*E*1*E*<sup>2</sup> is the creep time constant (also called the retardation time constant); *E*1, *E*2, and *η* are parameters relevant to viscoelastic properties; *σ*<sup>0</sup> is the stress value of an element within the material following the step stress excitation at *t* = 0 (i.e., the beginning of creep); and *σ*<sup>0</sup> is the constant stress value during creep. In the present study, *σ*<sup>0</sup> in Equation (1) is assumed to be the magnitude of uniform compressional pressure applied on the surface of the material, although this assumption may not be perfectly accurate for each element within the material since the appropriateness of this assumption depends on the loading and boundary conditions. ଵܧ <sup>ଶ</sup>ܧ + <sup>ଵ</sup>ܧ where ߳(ݐ (is the strain over time of an element within the material; ߬ = ߟ)ܧ<sup>ଵ</sup> + ܧ<sup>ଶ</sup> )/ <sup>ଶ</sup>ܧଵܧ is the creep time constant (also called the retardation time constant); ܧ<sup>ଵ</sup> <sup>ଶ</sup>ܧ , ߟ and , are parameters relevant to viscoelastic properties; ߪ is the stress value of an element within the material following the step stress excitation at ݐ = 0) i.e., the beginning of creep); and ߪ is the constant stress value during creep. In the present study, ߪ in Equation (1) is assumed to be the magnitude of uniform compressional pressure applied on the surface of the material, although this assumption may not be perfectly accurate for each element within the material since the appropriateness of this assumption depends on the

If Equation (1) is used to curve-fit the creep curve of each element within the material, the *E*1, *E*2, and *τ<sup>C</sup>* of each element can be quantitatively evaluated. Consequently, the spatial distribution of the viscoelastic properties of the material can be obtained, as described in the following section. loading and boundary conditions. If Equation (1) is used to curve-fit the creep curve of each element within the material, the ܧ<sup>ଵ</sup> <sup>ଶ</sup>ܧ , , and ߬ of each element can be quantitatively evaluated. Consequently, the spatial distribution of the viscoelastic properties of the material can be obtained, as de-

scribed in the following section.

**Figure 2.** *Cont.*

**Figure 2.** The design of a compressional viscoelastography system with different experimental settings. (**a**) The surface area of the pressure plate is half of that of the top surface of the sample, and the sample is placed directly on the testing platform without any fixation. This setup is associated with the condition in Figure 5d. (**b**) The surface area of the pressure plate is half of that of the top surface of the sample, and the sample is secured by a sample container. This setup is associated with the condition in Figure 5f. (**c**) The surface area of the pressure plate is equal to that of the top surface of the sample, and the sample is placed directly on the testing platform without any fixation. This setup is associated with the condition in Figure 5a. (**d**) The surface area of the pressure plate is equal to that of the top surface of the sample, and the sample is secured by a sample container. This setup is associated with the condition in Figure 5c. **Figure 2.** The design of a compressional viscoelastography system with different experimental settings. (**a**) The surface area of the pressure plate is half of that of the top surface of the sample, and the sample is placed directly on the testing platform without any fixation. This setup is associated with the condition in Figure 5d. (**b**) The surface area of the pressure plate is half of that of the top surface of the sample, and the sample is secured by a sample container. This setup is associated with the condition in Figure 5f. (**c**) The surface area of the pressure plate is equal to that of the top surface of the sample, and the sample is placed directly on the testing platform without any fixation. This setup is associated with the condition in Figure 5a. (**d**) The surface area of the pressure plate is equal to that of the top surface of the sample, and the sample is secured by a sample container. This setup is associated with the condition in Figure 5c. *Materials* **2021**, *14*, x FOR PEER REVIEW 5 of 20

**Figure 3.** Illustration of the measurement processes of compressional viscoelastography. **Figure 3.** Illustration of the measurement processes of compressional viscoelastography.

**Figure 4.** Illustration of the Maxwell representation of the standard linear solid model. ܧଵ, ܧଶ, and ߟ are three viscoelastic properties. **Figure 4.** Illustration of the Maxwell representation of the standard linear solid model. *E*<sup>1</sup> , *E*2, and *η* are three viscoelastic properties.

### *2.2. Computational Simulations and Data Analysis 2.2. Computational Simulations and Data Analysis*

Finite element analysis using ABAQUS 2019 (Dassault Systems, Simulia Corporation, Johnson, RI, USA) is used to investigate the measurement accuracy of compressional viscoelastography when measuring the viscoelastic properties of the material under different loading and boundary conditions. Finite element analysis using ABAQUS 2019 (Dassault Systems, Simulia Corporation, Johnson, RI, USA) is used to investigate the measurement accuracy of compressional viscoelastography when measuring the viscoelastic properties of the material under different loading and boundary conditions.

Six simulation tests were run in this study. The specific settings regarding the loading and boundary conditions in each simulation test are described below and illustrated in Figure 5: Six simulation tests were run in this study. The specific settings regarding the loading and boundary conditions in each simulation test are described below and illustrated in Figure 5:


The six simulation tests are different in their loading and boundary conditions as described above, but they have some common settings. In each simulation test, an axisymmetric finite element model is used, and the radius and thickness of the axisymmetric model are 50 mm and 50 mm, respectively. Square finite elements (0.5 mm × 0.5 mm) are used to mesh the model. The radius and thickness of the imaging region (the region of the model where data are taken to produce the image, which is chosen as the region of the The six simulation tests are different in their loading and boundary conditions as described above, but they have some common settings. In each simulation test, an axisymmetric finite element model is used, and the radius and thickness of the axisymmetric model are 50 mm and 50 mm, respectively. Square finite elements (0.5 mm × 0.5 mm) are used to mesh the model. The radius and thickness of the imaging region (the region of the model where data are taken to produce the image, which is chosen as the region

of the model just below the transducer, with a length of 50 mm) are 25 mm and 50 mm, respectively, as shown in Figure 2a.

The model is made of linear viscoelastic material. The material is also assumed to be incompressible, isotropic, and homogeneous. The mechanical properties of the material are defined by four properties, including the modulus of elasticity (*E*), Poisson's ratio (set as a constant of 0.495, the maximum Poisson's ratio that can be set in ABAQUS), and two parameters in the one-branch dimensionless relaxation modulus:

$$g\_R(t) = 1 - g\left(1 - e^{\frac{-t}{T\_R}}\right) \tag{2}$$

where *g* is a material constant (0 < *g* < 1), and *τ<sup>R</sup>* is the relaxation time constant.

The loading rate for applying the compressional pressure is designed to be high (the time duration from zero to maximum pressure is 1/6 s) to mimic a step load. Once the maximum pressure (set as a constant of 1000 Pa in each simulation test) is reached, the maximum pressure is then maintained as a constant for a period of time. During this period, each element in the model exhibits creep behavior; that is, the strain response of each element increases with time until the steady state is reached. The creep curve of each element is sent to MATLAB (R2019a, Mathworks, Natick, MA, USA) for analysis. The *E*1, *E*2, and *τ<sup>C</sup>* of each element can be obtained by using Equation (1) to curve-fit the creep curve of each element. *η* can then be obtained by using the definition of the creep time constant *τ<sup>C</sup>* = *η*(*E*<sup>1</sup> + *E*2)/*E*1*E*2. It has been reported that *E* is equal to *E*1, while *g* is equal to *E*2/(*E*<sup>1</sup> + *E*2), and *τ<sup>R</sup>* is equal to *η*/*E*<sup>2</sup> [35]. Consequently, the three mechanical properties of each element (*E*, *g*, and *τR*) set in ABAQUS can be obtained.

Once the value of a mechanical property (*E*, *g*, or *τR*) of each element is assigned with a specific color, the 2D spatial distribution map of that mechanical property of the axisymmetric image region can be constructed (since the image region is quadrilateral and the elements within are square). However, this map of the axisymmetric image region is just half of the map that should be obtained in reality (since the model is axisymmetric). Therefore, this map of the axisymmetric image region (called the original map here) is reflected, and then combined with the original map to produce a full map.

The corresponding error map is also constructed for each mechanical property. The error value at each element in the error map is calculated by:

$$\text{error} = \frac{|\text{simulation value} - \text{theoretical value}|}{\text{theoretical value}} \tag{3}$$

where the simulation value is the value of the mechanical property of an element obtained from the simulation, and the theoretical value is the value of the mechanical property set in ABAQUS. If the error is less than 10% [36], the measurement is considered to be accurate.

In each simulation test, the mechanical properties are fixed and set as *E* = 10 kPa, *τ<sup>R</sup>* = 5 s, and *g* = 0.8 (because, in our preliminary study, it was found that the pattern of the map of a mechanical property is similar regardless of the mechanical properties; in other words, the pattern of the map of a mechanical property is sensitive to the loading and boundary conditions but insensitive to the mechanical properties). In the 2D spatial distribution map of each mechanical property in each simulation test, the percentage of the region in the map consisting of elements having the simulation value within ±10% [36] of the theoretical value set in ABAQUS is calculated. The larger this percentage, the more accurate the measurement, because the wider the region in the 2D spatial distribution map, the more accurate the value of the mechanical property.

ଶܧ

<sup>ଶ</sup>ܧ + <sup>ଵ</sup>ܧ)ߟ = ߬ constant

equal to ܧଶ/(ܧ<sup>ଵ</sup> + ܧ<sup>ଶ</sup>

*Materials* **2021**, *14*, x FOR PEER REVIEW 8 of 20

parameters in the one-branch dimensionless relaxation modulus:

݃ோ

<sup>ଶ</sup>ܧଵܧ/(

), and ߬ோ is equal to ߟ/ܧ<sup>ଶ</sup>

cal properties of each element (ܧ, ݃, and ߬ோ) set in ABAQUS can be obtained.

The model is made of linear viscoelastic material. The material is also assumed to be incompressible, isotropic, and homogeneous. The mechanical properties of the material are defined by four properties, including the modulus of elasticity (ܧ(, Poisson's ratio (set as a constant of 0.495, the maximum Poisson's ratio that can be set in ABAQUS), and two

The loading rate for applying the compressional pressure is designed to be high (the time duration from zero to maximum pressure is 1/6 s) to mimic a step load. Once the maximum pressure (set as a constant of 1000 Pa in each simulation test) is reached, the maximum pressure is then maintained as a constant for a period of time. During this period, each element in the model exhibits creep behavior; that is, the strain response of each element increases with time until the steady state is reached. The creep curve of each element is sent to MATLAB (R2019a, Mathworks, Natick, MA, USA) for analysis. The ܧ<sup>ଵ</sup>

, and ߬ of each element can be obtained by using Equation (1) to curve-fit the creep curve of each element. ߟ can then be obtained by using the definition of the creep time

. It has been reported that ܧ is equal to ܧ<sup>ଵ</sup>

<sup>ఛ</sup>ೃ൰ ( 2 )

[35]. Consequently, the three mechani-

,

, while ݃ is

(ݐ = (1 − ݃ ൬1 − ݁ି௧

where ݃ is a material constant (0 < ݃ < 1), and ߬ோ is the relaxation time constant.

**Figure 5.** Illustration of the loading and boundary conditions in each simulation test. (**a**) The uniform compressional pressure is exerted on the entire top surface. The bottom is fixed along the vertical direction while the side is not fixed. (**b**) The uniform compressional pressure is exerted on the entire top surface. The bottom is fixed along all directions while the side is not fixed. (**c**) The uniform compressional pressure is exerted on the entire top surface. The bottom is fixed along the vertical direction while the side is fixed along the horizontal direction. (**d**) The uniform compressional pressure is exerted on half of the top surface. The bottom is fixed along the vertical direction while the side is not fixed. (**e**) The uniform compressional pressure is exerted on half of the top surface. The bottom is fixed along all directions while the side is not fixed. (**f**) The uniform compressional pressure is exerted on half of the top surface. The bottom is fixed along the vertical direction while the side is fixed along the horizontal direction.
