Solid-state refrigeration of ferroelectric materials has shown its low cooling energy density and high efficiency for miniaturization of electronic components. However, the low electrocaloric effect (ECE) and the limitation of the choice of ferroelectric materials hinder the application of the large ECE. With the giant ECE (12K) in the Pb(Zr,Ti)O₃ and polyvinylidene fluoride (PVDF) thin films (TFs) being obtained (Mischenko et al., 2006; Neese et al., 2008), the ferroelectric refrigeration has received more attention. After that, tremendous ECEs have been widely available in typical ferroelectrics such as PbTiO3 (PTO) (Qiu and Jiang, 2008), BaTiO3 (BTO) (Huang et al., 2018a), and lead-free K_1-x Na _x NbO₃ (KNN) (Zhao et al., 2020) by doping elements. A large number of experiments have found that the phase transition temperature of ferroelectric materials will sharply reduce to room temperature by element doping, which will also increase the operating temperature range of large ECEs .

With the development of experimental research on the ECE, multilayer composites (Qian et al., 2021; Liu et al., 2022a) and ceramic–polymer ferroelectric nanocomposites (Lich et al., 2020; Liu Y. et al., 2022b) can integrate the performance advantages of ferroelectric materials to strengthen the ECE and EC. Recently, molecular ferroelectric ImClO₄ and two-dimensional ferroelectric CuInP₂S₆ have shown excellent performance in electrocaloric refrigeration. Superior EC strength ( ∆ T / ∆ E = 0.84   K m M V − 1 ) has been attained in molecular ferroelectric ImClO₄ with the unique polarization mechanism generated by the order–disorder behavior (Li et al., 2020). With the van der Waals layered structure in two-dimensional ferroelectric CuInP₂S₆, a large EC strength ( ∆ T / ∆ E = 0.295   K m M V − 1 ) has been achieved (Si et al., 2019). These new electrocaloric materials are expected to become the development trend of high-performance ferroelectric solid-state refrigeration materials.

Different from ECEs actuated by ferroelectric–paraelectric phase transitions from the intrinsic factors of ion doping and material lamination, external stress–strain regulation can also significantly affect the change of the domain structure, eventually leading to large ECEs (Wang et al., 2014; Liu and Wang, 2015). The large ECEs generated by multidomain-to-monodomain transitions of PTO ferroelectrics under mismatch strains have been obtained by phase-field simulation (Hou et al., 2018). Most of the results show that in-plane polarization is produced by tensile mismatch strain, while out-of-plane polarization is affected by compressive strain. Different from the mismatch strain, the appropriate stress can also enhance the EC strength and increase the working temperature range of high EC strength (Wu et al., 2015).

Furthermore, the designed shape and structure of ferroelectrics can dramatically affect the ECEs. When the ferroelectric material contains cracks under the action of the electric field, a high electrostatic energy is generated in the cracks and the stress concentration near the tip of the cracks causes domain structure transition, thereby generating a giant ECE (Huang et al., 2021). So far, there are few research works on ferroelectric ECE with defects. In this paper, large ECEs induced by multidomain-to-monodomain transition in ferroelectrics with electrical inclusions are investigated based on the phase-field method. The main contents of the paper are arranged as follows: through the concise introduction, the phase-field model of ferroelectrics containing elliptical electric inclusion energy is given, and the unique electrical boundary conditions of the ferroelectric model containing electric inclusions are set in Section 2. The different giant ECEs formed by the ferroelectric model containing air, silicone oil, and water under the applied electric field are simulated, and the EC strength of the three models is analyzed in Section 3. Finally, the paper is summarized in Section 4.

Based on the Ginzburg–Landau phase transition theory, the total energy of the system is expressed as a function of time-dependent order parameters (Huang et al., 2018b). In particular, in order to simulate the influence of electrical inclusions on the evolution of the domain structure, the electrical inclusions are simulated as dielectrics filled with air, silicone oil, and water, which have negligible elastic energy. Then, the total electrical enthalpy h of the phase-field model includes Landau energy, gradient energy, stress–strain energy, electrostriction energy, electrostatic energy of ferroelectric solids, and electric inclusion, which can be expressed as follows: h P i , P i , j , ε i j , E i = a i P i 2 + a i j e P i 2 P j 2 + a i j k P i 2 P j 2 P k 2 + 1 2 G i j k l P i , j P k , l + 1 2 c i j k l ε i j ε k l − q i j k l ε i j P k P l − 1 2 κ 0 E i E i − E i P i − 1 2 κ c E i E i , (1) where a i , a i j e , and a i j k are the dielectric stiffness related to Landau energy and P i is the polarization regarded as the first-order parameter. G i j k l and P i , j denote the gradient coefficients and polarization gradient related to gradient energy, respectively. As the field variable of stress–strain energy, c i j k l and ε i j are expressed as the elastic stiffness constants and the strain, respectively. Similarly, q i j k l and E i are the field variables of electrostatic energy and denote the electrostrictive coefficients and electric field, respectively. It should be noted that κ 0 and κ c are represented as the dielectric constant of matrix and electric inclusion, respectively.

Based on the basic principle of thermodynamics and Ginzburg–Landau equation, the mechanical equilibrium equation, electrical equilibrium equation, and phase-field dynamic equation can be expressed as follows (Huang et al., 2018a): ∂ ∂ x j ∂ h ∂ ε i j = 0 , (2) ∂ ∂ x j − ∂ h ∂ E i = 0 , (3) ∂ P ∂ t = − L ∂ h ∂ P i − ∂ h ∂ P i , j , j , (4) where x j and t represent spatial and temporal, respectively. L represents the dynamic coefficient of the phase field. Here, the case with mechanical volume force and body charge is not considered. Thus, by integrating Eqs 2–4, the following weak solution forms are obtained: ∫ V ∂ h ∂ ε i j δ ε i j + ∂ h ∂ E i δ E i + 1 L ∂ P i ∂ t + ∂ h ∂ P i δ P i + ∂ h ∂ P i , j δ P i , j d V = ∫ s T i δ u i − ω δ ϕ + π i δ P i d S , (5) where u i and ϕ , as basic constant variables, represent the displacement and electric potential, respectively. T i and ω , as dependent variables, denote the surface traction and surface charge, respectively. In order to simplify the expression, π i is used instead of ∂ h ∂ P i , j n j to express the surface gradient energy flux. Based on the variational principle, the non-linear multifield coupling finite element method is adopted (Wang and Kamlah, 2009). The three-dimensional 8-node isoparametric element is used for spatial discretization, and the backward Euler iteration method is used for time integration to solve the non-linear equations. Finally, the physical field solutions are obtained by post-processing programming.

Based on the phenomenological theory and thermodynamic principle, through the Maxwell relationship, when the external electric field is applied, the adiabatic temperature change ( Δ T ) is calculated by Δ T = − T ∫ E a E b 1 C E T , E ∂ P i ∂ T E d E i (6) where C E T , E is the thermal capacity of ferroelectric materials. Based on the relevant experimental results (Liu and Wang, 2015), the thermal capacity of PTO can be set as a constant ( 3.9 × 10 6 J m − 3 K ) in this paper by ignoring the effect of the polarization field.

Figure 1 shows the overall schematic diagram of the ferroelectric plate with an electric inclusion under the applied electric field. As the phase-field model of this paper, the dimensionless sizes are set as follows: the length of the ferroelectric plate is 120, the width is 90, and the thickness is 1. The minimal thickness is set to degenerate the 3-D problem into the plane stress problem, which greatly simplifies the calculation process. The electric inclusion at the center is a flat ellipse, the long half axis of the electric inclusion is 20, and the short half axis is 1. It should be pointed out that in order to only reflect the influence of electrical inclusions on the domain structure, zero stress or strain loading is applied. The upper and lower electric potential values are set as V e , * and − V e , * , respectively, where superscripts e and * denote applied loading and dimensionless variable, respectively. Thus, the electric loads can be gained as E e , * = V e , * / 45 . The electrical boundary conditions of other surfaces can be set as the electrical open circuit (i.e., D i n i = 0 ). Since the model does not consider the application of mechanical loads, the traction-free condition is applied on each surface (i.e., σ i j n j = 0 ). The boundary condition for the polarization field is set to polarization-free on each surface ( i . e ., ∂ h ∂ P i , j n j = 0 ) .

Figure 2 shows the mesh division of the model and the mesh densification at the crack tip, which results in the refined mesh size being much smaller than the domain wall size to ensure the correctness of the phase-field model. The material parameters of PTO in this paper are consistent with those in the previous work (Huang et al., 2018b). The dielectric constant of the medium in inclusion as air, oil, or water can be set as κ c , a i r = 8.854 × 10 − 12 C V − 1 m − 1 , κ c , o i l = 2.12 × 10 − 11 C V − 1 m − 1 , and κ c , w a t e r = 7.099 × 10 − 10 C V − 1 m − 1 , respectively. In order to be consistent with the material parameters of dimensionless stiffening, we set the following dimensionless field variables: E e , * = E e / E 0 E 0 = 1306   k V   c m − 1 , P y * = P y / P 0 P 0 = 0.757 C   m − 2 . (7)

The initial polarization field direction is randomly distributed and the initial value is two orders of magnitude lower than that of the simulated polarization field, which is convenient to simulate the growth and enlargement of the ferroelectric domain structure, and finally form a steady-state evolution process. The dimensionless time step of the computational simulation is taken as ∆ t = 0.04 .

First, the ECE of the electric inclusion model containing air is simulated in this paper. Figure 3 shows the curve of the average polarization field along the y direction with the temperature under the electric loads. The results show that the polarization field increases suddenly at A ₁ - A ₂ and B ₁ - B ₂ , resulting in negative giant ECEs, which is consistent with the previous simulation results (Hou et al., 2018). The results also show that when the temperature exceeds 500°C, ferroelectric phase-to-paraelectric phase transition occurs in ferroelectric materials containing air inclusions, and the polarization field decreases sharply, resulting in positive giant ECEs.

In order to reveal the mechanism of large ECE through microstructure evolution, the detailed domain structures of PTO with an air inclusion are given in Figure 4 ( a ₁)—( d ₁), respectively. Under zero electric loads, the domain structure is in a head-to-tail state at different temperatures, which can greatly reduce the depolarization energy, as shown in Figure 4 ( a ₁) and ( c ₁). Figures 4 ( b ₁) and ( d ₁) exhibit the domain transition from the multidomain to monodomain state along the direction of the electric field ( E e , * = 0.7 ) during the ambient temperature rising from 90°C to 95°C, and only a winged 90° domain change is formed near the air inclusion. The giant strain energy and electrostatic energy could break through the energy barrier, and the sharp increase in the polarization field eventually leads to a negative giant ECE.

Figure 5 indicates that when the applied electric field changes, the average adiabatic temperature change has the largest negative ECE near the ambient temperature field ( E e , * = 0.5 ) of 150°C, and the temperature domain change produces a negative giant ECE ( ∆ T = − 9.79 ℃ ). When the temperature is close to the phase transition temperature of 500°C, a positive giant ECE is generated, and the temperature phase transition produces a positive giant ECE ( ∆ T = 14.83 ℃ ). By increasing of the applied electric loads, a high electrostatic energy is generated in the air inclusion, and as a result, the domain transition temperature moves sharply to room temperature. The large negative ECE is generated ( ∆ T = − 6.60 ℃ ) near the temperature field of 90°C. It can be found that although the negative ECE obtained by increasing the electric field is slightly reduced, the temperature is greatly reduced, which is more conducive to engineering practice.

As we know, the air inclusion is apt to dielectric breakdown and produces a discharge phenomenon under high applied electric loads. Therefore, the inclusion is filled with silicone oil or water, which can effectively protect the ferroelectric material and improve the tolerance of the applied electric field strength. Figure 6 shows the temperature dependence of the polarization component P ₂ of PTO with oil and water inclusions under different electric loads. The results show that when the electric field is applied, the polarization field of the model containing silicone oil inclusion increases abruptly at C ₁ - C ₂ and the water-inclusion model changes abruptly at D ₁ - D ₂ , all of which produce negative giant ECEs. Interestingly, the temperature for the giant ECE of the water-inclusion model is 50°C, which effectively controls the temperature around the room temperature, which is convenient for experimental and engineering applications.

The detailed domain structures of PTO with oil and water inclusions at C ₁ - C ₂ and D ₁ - D ₂ to explain the mechanism of the giant ECE are shown in Figures 7, 8, respectively. The evolution of the microstructure reveals that the multidomain-to-monodomain transition is the result of energy minimization of the whole system under the coupling multifield, including the electric inclusion with high electrostatic energy, which explains the macroscopic behavior of the P-T curve in Figure 6 from the microscopic mechanism.

The adiabatic temperature change under different electric loads is shown in Figure 9 for the oil-inclusion and water-inclusion models. As shown in Figure 9, the peak value of negative ATC of the water-inclusion model at 50°C with ∆ T = − 9.30 ℃ under relatively lower applied electric field ∆ E = 39.18 M V m − 1 . Table 1 summarizes the comparison of the giant ECE characterization parameters of the three electrical inclusion models with the literature results in detail. The results indicate that giant ECEs can be realized in the water-inclusion model near room temperature by effectively selecting the appropriate applied electric field increment. The current work may provide a fundamental understanding of domain transition and ECEs impacted by the electric inclusion and applied electric field.

In summary, the ECE of ferroelectric solids with different electric inclusions under electric loads is studied using a phase-field method. Numerical simulations show that the electric inclusion models containing air and silicone oil produce high electrostatic energy under applied electric field, and negative large ECEs are produced, but the domain transition temperature is much higher than room temperature. Because the dielectric constant of water and PTO are on the same order of magnitude, there is no electric field concentration inside the water inclusion, which only affects the distribution of the polarization field of the model as a crack; thus, a ferroelectric solid with a water inclusion under relatively lower applied electric field ∆ E = 39.18 M V m − 1 , an unusual ultrahigh negative ECE (about −9.30 K), and a large EC strength ∆ T / ∆ E = 0.237   K m M V − 1 are obtained near room temperature (about 50°C). The present study provides a feasible way to regulate ECEs near room temperature under low electric loads by controlling the dielectric constant of electric inclusions in ferroelectrics.