Modelling of Non-Newtonian Starved Thermal-elastohydrodynamic Lubrication of Heterogeneous Materials in Impact Motion

Abstract

This study presents a numerical model for the thermal-elastohydrodynamic lubrication of heterogeneous materials in impact motion, in which a rigid ball bounces on a starved non-Newtonian oil-covered plane surface of an elastic semi-infinite heterogeneous solid with inhomogeneous inclusions. The impact–rebound process and the microscopic response of the subsurface inhomogeneous inclusions are investigated. The inclusions are homogenized according to Eshelby’s equivalent inclusion method. The Elrod algorithm is adopted to determine the lubrication starvation based on the solutions of pressure and film thickness, while the lubricant velocity and shear rate of the non-Newtonian lubricant are derived by using the separation flow method. The dynamic response of the cases subjected to constant impact mass, momentum, and energy is discussed to reveal the influence of the initial drop height on the impact–rebound process. The results imply that the inclusion disturbs the subsurface stress field and affects the dynamic response of the contact system when the surface pressure is high. The impact energy is the decisive factor for the stress peak, maximum hydrodynamic force, and restitution coefficient, while the dynamic response during the early approaching process is controlled by the drop height.

1 Introduction

Impact motion is involved in many manufacturing processes, such as shot-peening and ball-forming, causing dynamic loads on the machine transmission components. Fluid lubricants are usually added to build an interface layer to improve the machine durability and alleviate the frictional loss. The rapid variation of loads results in the fluctuations of film thickness, pressure, and temperature. The lubrication in impact motion is caused by the normal squeeze motion when an impact load is applied, and the surface deformation and viscosity due to local high pressure should always be taken into account for the analysis of the lubrication film.

The first numerical model for the elastohydrodynamic lubrication (EHL) of materials in impact motion was developed by Christensen [1]. To reduce numerical difficulties, the squeeze action resulting from the change of deformation was neglected. However, this action was shown to play important roles in the dynamic response during the initial stage of the impact motion [2] and the whole impact–rebound process [3,4,5]. In these theoretical studies, the hydrodynamic force and restitution coefficient were applied to evaluate the energy dissipation. Additionally, their works showed that the lubricant was entrapped at the contact interface and formed a central dimple, while the minimum film thickness occurred at the contact periphery. Similar to the pressure–distance curve at the steady-state rolling condition, the pressure–time curve of the lubricant film has both a peak and a spike. Safa et al. [6] verified these results experimentally and further pointed out that the thickness of the central oil film in the contact remains constant in the major part of the impact–rebound process. Later, the EHL of the material in impact motion with a small initial impact gap was investigated numerically [7]. It was found that a periphery dimple occurs instead of the common central dimple, and the location of the minimum film thickness moves from the central region to the outer region. Within the impact–rebound cycle, the central film thickness remains almost constant during the load-sustaining stage. Recently, a unified Reynolds equation for the impact problem and the rolling/sliding problem was proposed based on the impact velocity [8]. The results showed that the impact velocity evolved linearly with time during a significant part of the impact process.

With the development of advanced computer and numerical simulation technologies, thermal effects [9] and lubrication starvation [10] were taken into account to improve the theoretical model. Wang et al. [9] proposed a model for Newtonian EHL of homogeneous materials in fully flooded conditions. The central pressure, film thickness, and temperature were investigated. It was concluded that the thermal effect was either influenced by the impact gap or the initial velocity. Wu et al. [10] rewrote the Reynolds equation by incorporating the Elrod algorithm and found that the rebound was delayed by the oil starvation. The starvation condition refers to a lubrication model, in which the amount of supplied oil is insufficient to fill the inlet gap of the lubricated conjunction. In the early numerical studies of EHL problems, the meniscus inlet boundary position was introduced as an known parameter to indicate the starvation conditions and predict the inlet oil layer thickness [11]. However, the meniscus inlet boundary position depends on the inlet oil thickness in reality. In 1981, Elrod [12] developed a numerical algorithm that can automatically determine the meniscus position based on the given inlet oil layer thickness. Since then, the studies on the starvation effect in EHL problems were mostly based on such an algorithm. The influence of initial drop height on the dynamic response has been explored with constant impact mass in many studies [6, 7, 9], while more efforts should be devoted to investigate the effects of the impact momentum and energy.

The above studies were all focused on Newtonian lubricants, assuming that the lubricant viscosity was independent of the lubricant shear strain rate to simplify the problems. However, in reality, lubricants always possess non-Newtonian properties including the shear-thinning and limiting shear stress behaviours. Many researchers have paid attention to non-Newtonian EHL problems in tangential motion conditions. A generalized Reynolds equation was derived by Yang and Wen [13] based on the equivalent viscosity that coupled the non-Newtonian effects and thermal effects. By applying such a method, non-Newtonian lubricants with different rheological properties were studied [14, 15]. The equivalent viscosity and shear rate were solved using a redundant iteration loop controlled by velocity boundary conditions. Later, Yang et al. [16, 17] put forward a separation flow method to derive analytical lubricant shear rate and velocity, which assumed the lubricant flow as the superposition of the Poiseuille flow and Couette flow. Mezlane et al. [18] studied the viscosity wedge deduced by the thermal effect and squeeze effect of an infinitely long cylinder in a transient non-Newtonian EHL problem by applying the Carreau–Yasuda rheological model. Recently, some experimental studies were conducted to investigate the rheological behaviours of non-Newtonian lubricants under the impact load [19, 20]. Nevertheless, few theoretical studies have been carried out for the non-Newtonian lubrication problems in impact motion conditions.

In previous studies, the materials in contact were commonly assumed to be homogeneous. However, materials are naturally inhomogeneous, consisting of micro-defects such as inclusions and voids. Therefore, it is critical to build a reliable model considering the material inhomogeneity. Micro-defects and their interactions within materials have been intensively studied [21, 22]. So far, limited studies have been conducted for materials with near-surface defects under contact loading with or without the oil film. Zhou et al. [23, 24] employed Eshelby’s equivalent inclusion method (EIM) to model inhomogeneous inclusions and developed a semi-analytical solution for materials with multiple inclusions under contact loading. The EIM was further employed to model the EHL problem of inhomogeneous materials [25]. Slack et al. [26] proposed an approach for computing the pressure, fluid film thickness, and subsurface stresses by employing the discrete element method. Later, the effects of multiple inclusions on the pressure profile, film thickness profile, and subsurface elastic fields were investigated [27, 28]. However, the impact motion, lubrication starvation, and non-Newtonian properties were not considered in the mentioned solutions.

The present study develops a numerical model for starved thermal-EHL of inhomogeneous material in impact motion. A rigid ball bounces on an oil-covered plane surface of an elastic semi-infinite heterogeneous solid with inhomogeneous inclusions. The pressure, film thickness, starvation parameter, temperature, and subsurface elastic fields during the impact–rebound process are simulated. The EIM, Elrod algorithm, and separation flow method are adopted sequentially to deal with the effects of the inhomogeneous inclusions, lubrication starvation, and non-Newtonian properties. The dynamic response of the cases subjected to constant impact mass, momentum, and energy is discussed to reveal the influence of the initial drop height on the impact–rebound process.

2 Mathematical Model

The problem considered is modelled as a rigid ball with a radius R that moves downwards under gravity and hits onto the oil-covered plane surface of an elastic semi-infinite heterogeneous solid with inclusions, as sketched in Fig. 1. The half-space is bounded by the plane surface $z = 0$ in an $x-y-z$ Cartesian coordinate system whose origin is set at the initial dry contact point. The half-space and ball have the densities $\rho _{\mathrm {1}}$ and $\rho _{\mathrm {2}}$, Young’s moduli $E_{\mathrm {1}}$ and $E_{\mathrm {2}}$, Poisson’s ratios $\upsilon _{\mathrm {1}}$ and $\upsilon _{\mathrm {2}}$, thermal conductivity $k_{\mathrm {1}}$ and $k_{\mathrm {2}}$, and heat capacity $c_{\mathrm {1}}$ and $c_{\mathrm {2}}$, respectively. The ball freely falls from the initial drop height $h_{{0}}^{{*}}$ onto the oil film with initial layer thickness $h_{\mathrm {oil}}$, density $\rho$, isotropic viscosity $\eta$, thermal conductivity $k_{\mathrm {f}}$, and heat capacity $c_{\mathrm {f}}$. It is noted that the impact ball applied in this study is a rigid one with an infinitely large Young’s modulus. The half-space contains n arbitrarily shaped sub-domains $\Omega _{\psi }$ ($\psi = 1, 2, \ldots , n)$, each of which has different mechanical material constants but the same thermal material constants with the matrix. Moreover, the influences of thermal deformation and stresses are small enough to be neglected due to the limiting contact area. The elastic moduli of the matrix and the subdomains are denoted by $C_{ijkl}$ and $C_{ijkl}^{\psi }$ ($i,\, j,\, k,\, l = x,\, y,\, z)$, respectively (see Ref. [11] for more details about the elastic moduli, i.e. the fourth-order stiffness tensor). The inhomogeneous inclusions are considered as the subdomains containing the initial eigenstrain $\varepsilon _{ij}^{0}$.

Fig. 1

Schematic of starved thermal-EHL of the inhomogeneous material subjected to impact

By incorporating the Elrod algorithm [12], the general Reynolds equation of the impact motion under the starved condition can be rewritten as

$$\begin{aligned} \frac{\partial }{\partial x}\left( {\varepsilon ^{{*}}\frac{\partial p}{\partial x}} \right) +\frac{\partial }{\partial y}\left( {\varepsilon ^{{*}}\frac{\partial p}{\partial y}} \right) -\frac{\partial \theta \rho _{1} }{\partial t}=0 \end{aligned}$$

(1)

where

$$\begin{aligned} {\varepsilon ^{{*}}}= & {} {\left( {\frac{\rho }{\eta }} \right) _\mathrm{{1}}}, \quad {\left( {\frac{\rho }{\eta }} \right) _\mathrm{{1}}} = \frac{{{\eta _\mathrm{{1}}}\rho _\mathrm{{1}}^\prime }}{{\eta _\mathrm{{1}}^\prime }} - \rho _\mathrm{{1}}^{\prime \prime }\nonumber \\ \rho _{{1}}= & {} \int _0^h \rho \mathrm{d}z,\quad \rho _{1}^{\prime } =\int _0^h \rho \int _0^z {\frac{\mathrm{d}z^{\prime }}{\eta }} \mathrm{d}z,\quad \rho _{1}^{\prime \prime } =\int _0^h \rho \int _0^z {\frac{z^{\prime }\mathrm{d}z^{\prime }}{\eta }} \mathrm{d}z\nonumber \\ \frac{1}{\eta _{1} }= & {} \int _0^h {\frac{\mathrm{d}z}{\eta }} , \frac{1}{\eta _{1}^{\prime } }=\int _0^h {\frac{z\mathrm{d}z}{\eta }} \end{aligned}$$

(2)

The fraction film content $\theta$, defined as the ratio of the film thickness $h_{\mathrm {L}}$ to the geometry gap h, is introduced to deal with the lubrication starvation. Therefore, the complementary condition is written as

$$\begin{aligned} p(1-\theta )=0 \end{aligned}$$

(3)

where 0 $\le \theta \le 1$ should be satisfied. The boundary and initial conditions for Eq. (1) are given as

$$\begin{aligned} \left\{ {{\begin{array}{c} {p\left( {x_{\mathrm{in}} ,y,t} \right) =p\left( {x_{\mathrm{out}} ,y,t} \right) =p\left( {x,y_{\mathrm{in}} ,t} \right) =p\left( {x,y_{\mathrm{out}} ,t} \right) =0} \\ {p(x,y,t)\geqslant 0 \quad \left( {x_{\mathrm{in}}<x<x_{\mathrm{out}} ,y_{\mathrm{in}}<y<y_{\mathrm{out}} } \right) } \\ \end{array} }} \right. \end{aligned}$$

(4)

The geometry gap h is defined as

$$\begin{aligned} h(x,y,t)=h_{0} (t)+\frac{x^{2}\text{+ }y^{2}}{2R}+U_{z} (x,y,t) \end{aligned}$$

(5)

The normal elastic deformation $U_{\mathrm {z}}$ can be decomposed into two parts: the surface deformation $U_{z}^{\mathrm {e}}$ induced by the external pressure and the disturbed deformation $U_{z}^{*}$ due to the eigenstrains. The surface deformation $U_{z}^{\mathrm {e}}$ is solved by the integral equation as

$$\begin{aligned} U_{z}^{\text{ e }} (x,y)=\frac{2}{\pi E^{{*}}}\iint {\frac{p\left( {x^{\prime },y^{\prime }} \right) }{\sqrt{\left( {x-x^{\prime }} \right) ^{2}+\left( {y-y^{\prime }} \right) ^{2}} }} \text{ d }x^{\prime }\text{ d }y^{\prime } \end{aligned}$$

(6)

where E* is the composite Young’s Modulus defined as

$$\begin{aligned} E^{*}=\left( {\frac{1-\nu _{1}^{2} }{E_{1} }+\frac{1-\nu _{2}^{2} }{E_{2} }} \right) ^{\text{- }1} \end{aligned}$$

(7)

The disturbed deformation $U_{z}^{*}$ is calculated based on the EIM by homogenizing the subsurface inclusions with unknown equivalent eigenstrains $\varepsilon _{ij}^{{*}}$. Further detailed calculations of the eigenstrains and the disturbed deformation can be found in “Appendix.”

The rigid body approach $h_{{0}}$ is controlled by dynamic relations as

$$\begin{aligned} \left\{ {{\begin{array}{c} \displaystyle {h_{0} (k\Delta t+\Delta t)=h_{0} (k\Delta t)+v_{0} (k\Delta t)\cdot \Delta t+\frac{1}{2}a_{0} (k\Delta t)\cdot (\Delta t)^{2}} \\ \displaystyle {v_{0} (k\Delta t+\Delta t)=v_{0} (k\Delta t)+a_{0} (k\Delta t)\cdot \Delta t} \\ \displaystyle {a_{0} (k\Delta t+\Delta t)=\frac{1}{m_{0} }\int \!\!\!\int _ p (x,y,k\Delta t+\Delta t)\text{ d }x\text{ d }y-g} \\ \end{array} }} \right. \end{aligned}$$

(8)

where k is the number of the time interval in the calculation, and $v_{{0}}$ and $a_{{0}}$ are the velocity and acceleration of the ball, respectively. By assuming the ball falls freely at $t = 0$, the initial conditions of these dynamic relations are

$$\begin{aligned} \left\{ {{\begin{array}{c} {h_{0} (0)=h_{\mathrm{oil}} } \\ {v_{0} (0)=-\sqrt{2g\left( {h_{0}^{*} -h_{\mathrm{oil}} } \right) } } \\ {a_{0} (0)=-g} \\ \end{array} }} \right. \end{aligned}$$

(9)

The temperature variation of the oil film is controlled by the energy equation as

$$\begin{aligned}&c_{\mathrm{f}} \left( {\rho \frac{\partial T}{\partial t}+\rho u\frac{\partial T}{\partial x}+\rho v\frac{\partial T}{\partial y}-q\frac{\partial T}{\partial z}} \right) =k_{\mathrm{f}} \frac{\partial ^{2}T}{\partial z^{2}} \nonumber \\&\quad -\frac{T}{\rho }\frac{\partial \rho }{\partial T}\left( {\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}+v\frac{\partial p}{\partial y}} \right) +\eta \left[ {\left( {\frac{\partial u}{\partial z}} \right) ^{2}+\left( {\frac{\partial v}{\partial z}} \right) ^{2}} \right] \end{aligned}$$

(10)

where

$$\begin{aligned} q=\frac{\partial }{\partial x}\int _0^z \rho u\mathrm{d}z^{\prime }+\frac{\partial }{\partial y}\int _0^z \rho v\mathrm{d}z^{\prime }+\frac{\partial }{\partial t}\int _0^z \rho \mathrm{d}z^{\prime } \end{aligned}$$

(11)

The initial and boundary conditions of Eq. (10) read

$$\begin{aligned} \left. {T(x,y,z,t)} \right| _{t=0} =T_{0} \end{aligned}$$

(12)

and

$$\begin{aligned} \left\{ {{\begin{array}{*{20}c} {T\left( {x_{\mathrm{in}} ,y,z,t} \right) =0} \\ {T\left( {x,y_{\mathrm{in}} ,z,t} \right) =0} \\ \end{array} }} \right. \end{aligned}$$

(13)

respectively.

The energy equations for both solids are

$$\begin{aligned} \left\{ {{\begin{array}{c} {c_{1} \rho _{1} \frac{\partial T}{\partial t}=k_{1} \left( {\frac{\partial ^{2}T}{dx^{2}}+\frac{\partial ^{2}T}{\partial y^{2}}+\frac{\partial ^{2}T}{\partial z_{1}^{2} }} \right) } \\ {c_{2} \rho _{2} \frac{\partial T}{\partial t}=k_{2} \left( {\frac{\partial ^{2}T}{\partial x^{2}}+\frac{\partial ^{2}T}{\partial y^{2}}+\frac{\partial ^{2}T}{\partial z_{2}^{2} }} \right) } \\ \end{array} }} \right. \end{aligned}$$

(14)

with the initial conditions

$$\begin{aligned} \left\{ {{\begin{array}{c} {T(x,y,z_{1} ,t)\left| {_{t=0} =T_{0} } \right. } \\ {\left. {T(x,y,z_{2} ,t)} \right| _{t=0} =T_{0} } \\ \end{array} }} \right. \end{aligned}$$

(15)

and boundary conditions

$$\begin{aligned}&T(x_{\mathrm{in}} ,y,z_{1} ,t)=T_{0} , \,\, T(x,y_{\mathrm{in}} ,z_{1} ,t)=T_{0} ,\,\, T(x_{\mathrm{out}} ,y,z_{1} ,t)=T_{0} ,\,\, T(x,y_{\mathrm{out}} ,z_{1} ,t)=T_{0} \nonumber \\&T(x_{\mathrm{in}} ,y,z_{2} ,t)=T_{0} ,\,\, T(x,y_{\mathrm{in}} ,z_{2} ,t)=T_{0} ,\,\, T(x_{\mathrm{out}} ,y,z_{2} ,t)=T_{0},\,\,T(x,y_{\mathrm{out}} ,z_{2} ,t)=T_{0} \nonumber \\&T(x,y,-d,t)=T_{0} ,\,\, T(x,y_{\mathrm{out}} ,d,t)=T_{0} \end{aligned}$$

(16)

The symbol $d = 3.15a_{\mathrm {Hz}}$ represents the depth of the computational domain ($a_{\mathrm {Hz}}$ is Hertzian contact radius) [29], and $z_{\mathrm {1}}$ and $z_{\mathrm {2}}$ represent the z-axes of the two solids, respectively. On the two oil–solid interfaces, the heat flux continuity conditions given by

$$\begin{aligned} \left\{ {{\begin{array}{c} {\left. {k_{1} \frac{\partial T}{\partial z_{1} }} \right| _{z_{1} =0} =\left. {k\frac{\partial T}{\partial z}} \right| _{z=0} } \\ {\left. {k_{2} \frac{\partial T}{\partial z_{2} }} \right| _{z_{2} =0} =\left. {k\frac{\partial T}{\partial z}} \right| _{z=h} } \\ \end{array} }} \right. \end{aligned}$$

(17)

must be satisfied.

The lubricant density $\rho$ (p, T) is specified by the Dowson–Higginson relation [30] as

$$\begin{aligned} \rho (p,T)=\rho _{\mathrm{R}} \left[ {1+\frac{0.6\times 10^{-9}p}{1+1.7\times 10^{-9}p}-6.5\times 10^{-4}\left( {T-T_{R} } \right) } \right] \end{aligned}$$

(18)

where p is the pressure in GPa, $\rho _{\mathrm {R}}$ is the reference density at the ambient condition, and T is the temperature in K.

Numerous rheological models based on experimental data have been proposed to describe the shear-thinning and/or limiting shear stress behaviours, such as the Eyring model, Carreau–Yasuda model, Elsharkawy–Hamrock model, and Bair–Winer model, among which the Bair–Winer model has a broad application when comes to modelling the EHL problems with non-Newtonian lubricants [11]. The Bair–Winer model is adopted in this study and written as [31]

$$\begin{aligned} \frac{\tau }{\tau _{\mathrm{L}} }=1-\exp \left( {-\dot{{\gamma }}\frac{\eta _{\mathrm{R}\,\mathrm{oe}}}{\tau _{\mathrm{L}} }} \right) \end{aligned}$$

(19)

The limiting shear stress $\tau _{\mathrm {L}}$, the typical stress beyond which the limiting shear stress behaviour significantly affects the lubrication, can be obtained by

$$\begin{aligned} \tau _{L} =\tau _{L0} \left( {1+\frac{0.05}{\tau _{L0} }p} \right) \exp \left[ {\beta \left( {\frac{1}{T+273.15}-\frac{1}{T_{0} +273.15}} \right) } \right] \end{aligned}$$

(20)

The initial limiting shear stress $\tau _{\mathrm {L0}} = 2\, \hbox {MPa}$ and the temperature coefficient $\beta = 585\, \hbox {K}$ under ambient conditions [14]. Additionally, the Roelands viscosity $\eta _{\mathrm {Roe}}$ is employed to specify the viscosity–pressure–temperature relationship and written as

$$\begin{aligned} \eta _{\mathrm{R}\,\mathrm{oe}} =\eta _{\mathrm{R}} \exp \left\{ {\left( {\ln \eta _{\mathrm{R}} +9.67} \right) \left[ {-1+\left( {1+5.1\times 10^{-9}p} \right) ^{Z_{P} }\times \left( {\frac{T-138}{T_{\mathrm{R}} -138}} \right) ^{-s}} \right] } \right\} \end{aligned}$$

(21)

where

$$\begin{aligned} Z_{P} =\frac{2.19\times 1.98}{\ln \eta _{\mathrm{R}} +9.67},\,\, S=\frac{\beta \left( {T_{\mathrm{R}} -138} \right) }{\ln \eta _{\mathrm{R}} +9.67} \end{aligned}$$

(22)

and $\eta _{\mathrm {R}}$ is the reference lubricant viscosity at the ambient condition.

According to the separation flow method, the Couette flow is hardly influenced by the non-Newtonian properties, while the velocity distribution of the Poiseuille flow is symmetrical about the centreline and the value is equal to zero at boundaries [17]. Thus, the shear rate $\dot{\gamma }$ of the lubricant in the rheological model can be derived and written as

$$\begin{aligned} \left\{ {{\begin{array}{c} {\dot{{\gamma }}_{zx} {=}\frac{\partial u}{\partial z}=f\left( {\frac{\partial p}{\partial x}\left( {z-\frac{h}{2}} \right) ,\eta _{\mathrm{R}\,\mathrm{oe}} } \right) } \\ {\dot{{\gamma }}_{zy} {=}\frac{\partial v}{\partial z}=f\left( {\frac{\partial p}{\partial y}\left( {z-\frac{h}{2}} \right) ,\eta _{\mathrm{R}\,\mathrm{oe}} } \right) } \\ \end{array} }} \right. \end{aligned}$$

(23)

Subsequently, the shear stress $\tau$ is obtained by introducing the predicted shear rate $\dot{\gamma }$ into the rheological equations, and the lubricant viscosity $\eta$ is calculated using $\tau$ divided by $\dot{\gamma }$. The lubricant velocity is obtained by integrating the shear rate with respect to z. The resulted shear rate and velocity are written as

$$\begin{aligned} \left\{ {{\begin{array}{c} {\frac{\partial u}{\partial z}=\frac{\tau _{\mathrm{L}} }{\eta _{\mathrm{R}\,\mathrm{oe}} }\ln \left( {1-\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial x}\left( {z-\frac{h}{2}} \right) } \right) } \\ {\frac{\partial v}{\partial z}=\frac{\tau _{\mathrm{L}} }{\eta _{\mathrm{R}\,\mathrm{oe}} }\ln \left( {1-\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial y}\left( {z-\frac{h}{2}} \right) } \right) } \\ \end{array} }} \right. \end{aligned}$$

(24)

and

$$\begin{aligned} \left\{ {{\begin{array}{c} {u=-\frac{\tau _{\mathrm{L}} }{\eta _{\mathrm{R}\,\mathrm{oe}} }\frac{\tau _{\mathrm{E}} }{\frac{\partial p}{\partial x}}\left\{ {-\left( {1-\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial x}\left( {z-\frac{h}{2}} \right) } \right) \ln \left( {1-\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial x}\left( {z-\frac{h}{2}} \right) } \right) -\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial x}\left( {z-\frac{h}{2}} \right) } \right\} } \\ {v=-\frac{\tau _{\mathrm{L}} }{\eta _{\mathrm{R}\,\mathrm{oe}} }\frac{\tau _{\mathrm{E}} }{\frac{\partial p}{\partial y}}\left\{ {-\left( {1-\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial y}\left( {z-\frac{h}{2}} \right) } \right) \ln \left( {1-\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial y}\left( {z-\frac{h}{2}} \right) } \right) -\frac{1}{\tau _{\mathrm{L}} }\frac{\partial p}{\partial y}\left( {z-\frac{h}{2}} \right) } \right\} } \\ \end{array} }} \right. \end{aligned}$$

(25)

respectively.

3 Numerical Techniques

To formulate the Reynolds and elastic deformation equations in Sect. 2, a cubic computational domain is set at $x_{\mathrm {in}} \le x \le x_{\mathrm {out}}$, $y_{\mathrm {in}} \le y \le$ $y_{\mathrm {out}}$, and 0 $\le z \le z_{\mathrm {btm}}$ with $N_{x} \times N_{y} \times N_{z}$ cuboid elements having the same size of $\Delta _{x} \times \Delta _{y} \times \Delta _{z}$. Each element is indexed by a sequence of three integers [$\alpha _0 ,\, \beta _0 ,\, \gamma _0$] ($0 \le \alpha _0 \le N_{x} -1, \, 0 \le \beta _0 \le N_{y} -1, 0 \le \gamma _0 \le N_{z}-1$). The Poiseuille flow terms and squeeze term of the Reynolds equation are discretized using the second-order central scheme and first-order backward scheme, respectively. The temperature computation is implemented in a cubic domain with a depth range of $7.3a_{\mathrm {Hz}}$, which is discretized by $N_{\mathrm {zf}}$ equidistance nodes in the z-direction across the film, and $N_{\mathrm {zs}}$ non-equidistant nodes in the $z_{\mathrm {1}}$- and $z_{\mathrm {2}}$-directions within the half-space and ball, respectively. The element lengths along the x- and y-directions are also fixed at $\Delta _{x}$ and $\Delta _{y}$, respectively (see Ref. [29] for more details). The energy equations are discretized using the first-order backward scheme.

The entire problem couples the lubricated contact and the homogenization of heterogeneous materials. An iterative algorithm as shown in Fig. 2 is developed to obtain solutions of the contact pressure, film thickness, fractional film content, temperature, and subsurface stress field at each time instant. The calculational procedure is comprised of the computations of pressure, temperature, and disturbed deformation.

Fig. 2

Iterative algorithm for solving the problem of thermal-EHL of materials in impact motion at each time instant

First, pressure p, film thickness h, fractional film content $\theta$, and equivalent density $\rho _{\mathrm {1}}$ at the previous step and dynamic parameters $h_0,\, v_0$, and $a_0$ at the current step are read in as the initial values. The Reynolds equation, film thickness equation, and Elrod algorithm are solved sequentially by assuming the temperature and disturbed deformation are known. This iterative procedure is repeated until the relative difference between the summations of the pressure from two adjacent steps is smaller than 10$^{{-3}}$. Second, the new values of pressure and film thickness are taken into the temperature iteration to iteratively obtain the new temperature distribution until the relative difference of the temperature summations is smaller than $10^{{-3}}$. Third, the eigenstrain $\varepsilon _{ij}^{{*}}$ and eigenstress $\sigma _{ij}^{{*}}$ are calculated according to the EIM until the relative difference of the eigenstrain summations is smaller than $10^{{-6}}$. Finally, the disturbed deformation is calculated for the given pressure and eigenstrain. The global loop at the current time instant is terminated until the relative difference of the maximum deformation by the eigenstrain is smaller than $10^{\mathrm {-3}}$.

The differential Reynolds equation and energy equation are solved using the successive overrelaxation (SOR) method, while the governing equation of the equivalent homogeneous inclusions is solved using the conjugate gradient method (CGM). Discrete convolution and fast Fourier transform (DC-FFT) with duplicated padding technique are used to determine the surface deformation and stress caused by the eigenstrains (see Ref. [24] for more details). The multi-grid (MG) method, which can be considered as a hierarchical SOR algorithm, was widely used in many studies of EHL problems to improve the computational efficiency [9, 10]. Several discretization grids of different densities were used alternately following a certain cyclic pattern in order to accelerate the solution process and efficiently reduce the numerical errors of various wavelengths. However, it has been reported by Zhu [32] and Pu et al. [33] that the MG algorithm may not be desirable for EHL cases with thin film and mixed lubrication conditions. Furthermore, because the most time-consuming process of the calculation is the iterative determination of the eigenstrains in the inhomogeneity zones, the time consumption is almost at the same level whether employing the MG method or not. In this study, the computational grid is fixed at a single level, and the FFT is utilized to improve the calculation efficiency.

4 Results and Discussions

4.1 Operating Conditions

In this section, the contact pressure, film thickness, fractional film content, temperature, and subsurface elastic fields of a rigid ball that freely falls onto an elastic half-space with a single cuboidal stiff inclusion are simulated. The Reynolds equation and elastic deformation equations are implemented at the computational domain set at $-3.0a_{\mathrm {Hz}} \le x \le 3.0a_{\mathrm {Hz}}$, $-3.0a_{\mathrm {Hz}} \le y \le 3.0a_{\mathrm {Hz}}$, and $0 \le z \le 2.0a_{\mathrm {Hz}}$ with $129 \times 129 \times 32$ discretization grids equally spaced. In the temperature computation domain, the node numbers across the oil film and solids along the z-direction are 10 and 6, respectively. The time steep $\Delta t$ is $0.05\ \upmu$s in this study.

According to Ref. [34], a macro-inclusion with an edge length of 5–80 $\upmu$m may be formed in interstitial free steels and cause severe surface quality problems. A single inclusion with its edge length of $0.75a_{\mathrm {Hz}}$ (about 64.5 $\upmu$m) is focused on in this study, as shown in Fig. 3. It has been realized that the inclusion may influence the contact within a certain depth range; otherwise, the inclusion may not induce any disturbance to the contact. In this case, the inclusion is centred at (0, 0, $0.75a_{\mathrm {Hz}}$ ) to reveal its effect on the impact contact.

Fig. 3

Schematic of a cuboidal inclusion beneath a half-space surface

To discuss the effect of the initial drop height of the ball with constant impact mass, momentum, and energy on the dynamic response of the system, seven cases with different values of impact mass $m_0$ and initial drop height $h_{{0}}^{{*}}$ are studied as depicted in Table 1. The universal geometric and material parameters of all the cases are listed in Table 2.

Table 1 Impact parameters of different cases

Table 2 Geometric and material parameters

To analyze the results conveniently, the dimensionless coordinate system and temperature, respectively, normalized by the Hertzian contact radius $a_{\mathrm {Hz}}$ and ambient temperature $T_{\mathrm {R}}$ are presented in the following figures.

4.2 Model Verification

The central pressure, central film thickness, and central temperature obtained from the present model are compared with those obtained by Wang et al. [9] for the Newtonian lubricant, as shown in Fig. 4. The starvation condition is not considered in both models, meaning $\theta = 1$ in the simulations. The input parameters of the validation case have the same values as those summarized in Table 1, except that $E_{\mathrm {1}}= E_{\mathrm {2}} = 250$ GPa and $a_{\mathrm {Hz}} = 0.145\ \hbox {mm}$. It is observed that the results from the two models agree well with each other.

The results of the Bair–Winer lubricant are also displayed in Fig. 4, and they are consistent with those found in the study of the EHL problem in tangential motion conditions [13]. The non-Newtonian behaviours affect the pressure and thickness around the pressure spike region. The film thickness is increased by the large shear rate, which is induced by the non-Newtonian behaviours. The central temperature–time traces of Newtonian and non-Newtonian lubricants are quite different. Besides, the influence of the non-Newtonian behaviours is enhanced by the increase of initial drop height.

Fig. 4

Comparison between the results obtained by Wang et al. [9] and the present model: a central pressure, b central film thickness, and c dimensionless central temperature

4.3 Essential Features of Solutions

The pressure, film thickness, fractional film content, and mid-layer temperature profiles of case 2 at nine instants along the x-axis in the impact–rebound process are given in Figs. 5, 6, and  7. The results of the homogenous material under the same operating condition are also presented for comparison. At the first instant ($t = 0.05\ \upmu \mathrm{s}$) that the ball comes into contact with the oil, the pressure is not zero and fractional film content $\theta$ is 1.0 at that contact point. For the rest of the parts, the pressure is zero and $\theta < 1.0$. With the enlargement of the contact area, the pressure domain is expanded, while the area where $\theta = 1.0$ is enlarged accordingly. From 40 to $8 \ \upmu \mathrm{s}$, the contact pressure is rapidly increased and produces a significant elastic deformation to form the oil entrapment. The position of the minimum film thickness is located at the peripheral contact area. The pressured area is extended in the impact process until the majority of the $\theta$ profile achieves 1.0.

Fig. 5

Variations of the pressure and film thickness profiles in the x-direction: a $t = 0.05\,\upmu \mathrm{s}$, b $t = 10\,\upmu \mathrm{s}$, c $t = 40\ \upmu \mathrm{s}$, d $t = 80\ \upmu \mathrm{s}$, e $t = 100\ \upmu \mathrm{s}$, f $t = 115\ \upmu \mathrm{s}$, g $t = 130\ \upmu \mathrm{s}$, h $t = 135\ \upmu \mathrm{s}$, and i $t = 140\ \upmu \mathrm{s}$

As the rebound process starts, the oil at the periphery starts to flow out of the contact area so that the pressure at the corresponding position increases to balance the outflow. Specifically, the peripheral pressure is increased and there is a gap between the oil layer and the ball surface. As the process continues, $\theta$ is reduced gradually and the peripheral pressure is further increased. Outside the contact area, the gap between the ball surface and the oil surface is increased gradually. The elastic deformation at the periphery location is increased by the increase of the pressure so that a periphery entrapment is formed, as shown in Fig. 5g–i. As the rebound process continues, the minimum film thickness recovers to the contact centre.

Fig. 6

Variations of the fractional film content profiles in the x-direction: a $t = 0.05\ \upmu \mathrm{s}$, b $t = 10 \ \upmu \mathrm{s}$, c $t = 40 \ \upmu \mathrm{s}$, d $t = 80 \ \upmu \mathrm{s}$, e $t = 100 \ \upmu \mathrm{s}$, f $t = 115 \ \upmu \mathrm{s}$, g $t = 130 \ \upmu \mathrm{s}$, h $t = 135 \ \upmu \mathrm{s}$, and i $t = 140 \ \upmu \mathrm{s}$

Similar patterns can be found in the profile shapes of the pressure and mid-layer temperature distributions in Fig. 7, which is consistent with the solutions for the EHL of materials in tangential motion [13]. In the initial impact stage, the maximum mid-layer temperature occurs at the contact centre. In the subsequent stages, it shifts to the contact periphery in the rebound process and returns to the contact centre at last.

Fig. 7

Variations of the mid-layer temperature profiles in the x-direction: a $t = 0.05 \ \upmu \mathrm{s}$, b $t = 10 \ \upmu \mathrm{s}$, c $t = 40 \ \upmu \mathrm{s}$, d $t = 80 \ \upmu \mathrm{s}$, e $t = 100 \ \upmu \mathrm{s}$, f $t = 115 \ \upmu \mathrm{s}$, g $t = 130 \ \upmu \mathrm{s}$, h $t = 135 \ \upmu \mathrm{s}$, and i $t = 140 \ \upmu \mathrm{s}$

It is seen that the pressure, film thickness, fractional film content, and mid-layer temperature profiles of homogeneous and inhomogeneous materials in the x-direction are similar even if the pressure distributions around the inclusion are different. Moreover, the pressure of the inhomogeneous material shifts from the contact centre to the contact periphery at an earlier instant and the contact area is reduced due to the stiff inclusion at the same instant. It indicates that the impact duration is reduced by the stiff inclusion.

Fig. 8

Variations of the von Mises stress on the x–z plane and y–z plane: a $t = 0.05\ \upmu \mathrm{s}$, b $t = 10\ \upmu \mathrm{s}$, c $t = 40 \ \upmu \mathrm{s}$, d $t = 80\ \upmu \mathrm{s}$, e $t = 100\ \upmu \mathrm{s}$, f $t = 115\ \upmu \mathrm{s}$, g $t = 130\ \upmu \mathrm{s}$, h $t = 135\ \upmu \mathrm{s}$, and i $t = 140\ \upmu \mathrm{s}$

The inclusion-disturbed von Mises stress on the $x-z$ plane and $y-z$ plane is displayed in Fig. 8. The consistency of the stress on the $x-z$ plane and $y-z$ plane is due to the symmetry of this pure squeeze problem, which also demonstrates the reasonability of the developed numerical model. The stress concentrations are found along the inclusion edges where the crack nucleation is likely to start, and the von Mises stress inside is larger than that outside during the approaching process as the inclusion is stiffer than the substrate. The maximum stress moves from the inclusion centre to the top surface in the rebound process. At the end of the rebound process, the stress outside the inclusion is increased rapidly and beyond the stress value inside. The maximum von Mises stress occurs at the subsurface instead of the plane surface of the half-space. Consistent with the contact pressure, the maximum von Mises stress appears at the central vertical axis during the approaching process, then moves to the contact periphery during the rebound process, and finally returns to the central vertical axis abruptly.

The variations of the central pressure $p_{\mathrm {cen}}$, central film thickness $h_{\mathrm {cen}}$, minimum film thickness $h_{\mathrm {min}}$, central temperature $T_{\mathrm {cen}}$, and maximum von Mises stress $\sigma _{\mathrm {von}}$ with time are plotted in Fig. 9. The normal velocity $v_{{0}}$ of the ball and the hydrodynamic force $p_{z}$ throughout the impact–rebound process are also shown. The pressure spike in Fig. 9a for either the inhomogeneous material or the inhomogeneous material reaches a value of about 14 GPa, which is significantly higher than the yield strength of common metal materials and may produce surface damage in the rebound process. The central contact pressure is increased by the inclusion during the approaching process. According to energy conservation, the impact–rebound duration of the inhomogeneous material is shorter than that of the homogeneous material as shown in Fig. 9a.

Fig. 9

Dynamic response of homogeneous and inhomogeneous materials: a central pressure, b central film thickness, c minimum film thickness, d dimensionless central temperature, e maximum subsurface von Mises stress, and f ball velocity and hydrodynamic force

In Fig. 9b, c, both the central film thickness–time and minimum film thickness–time traces have a constriction formed to balance the pressure spike. The central film thickness $h_{\mathrm {cen}}$ of the inhomogeneous material is decreased gradually during the major part of the impact process (25–$125 \ \upmu \mathrm{s}$) rather than being remained constant as that of the homogeneous material. The contact pressure of the inhomogeneous material is concentrated around the region where the inclusion is embedded below, comparing with the solutions for the homogeneous material in Fig. 5a. Therefore, a part of the lubricant flows towards the contact periphery where the pressure is lower, making the central film thickness smaller, the minimum film thickness larger, and the pressure in the contact periphery smaller.

It is seen in Fig. 9d that the temperature is increased quickly during the initial approaching process, followed by a drop caused by the reduction of the film thickness before reaching the spike. The heat in the EHL of materials in impact motion is mainly produced by the compressive work that is highly related to the contact area. Since the contact area at the end of the rebound process is small, the spike of central temperature is not as significant as that of the central pressure.

As shown in Fig. 9e, the maximum von Mises stress–time trace of the inhomogeneous material has an obvious rise during the major part of the impact process (25 to $125 \ \upmu \mathrm{s}$) compared with that of the homogeneous material. The difference indicates that the inclusion mainly disturbs the stress field during the middle of the impact–rebound process when the contact pressure is relatively high. In the later part of the rebound process, a very small damping caused by the oscillation of the maximum contact pressure occurs, which is followed by a spike for both the inhomogeneous and homogeneous materials.

It is seen in Fig. 9f that the influences of the stiff inclusion on the ball velocity and hydrodynamic force are small enough to be neglected. The leaving velocity of the ball as it escapes from the 10-$\upmu \mathrm{m}$-thick layer of lubricant is estimated to be $0.280\ \mathrm{m}\cdot \mathrm{s}^{\mathrm {-1}}$, with the restitution coefficient of 0.893. The hydrodynamic force profiles in the approaching process and rebound process are approximately symmetrical about the rebound instant $t = 90 \ \upmu \mathrm{s}$.

4.4 Influence of the Initial Drop Height

4.4.1. Constant mass

Solutions for the drop heights of 2 mm (case 1), 5 mm (case 2), and 10 mm (case 3) are shown in Fig. 10. In the case of constant impact mass, with the increase of the initial drop height $h_{{0}}^{{*}}$, the peaks of the central pressure $p_{\mathrm {cen}}$, von Mises stress $\sigma _{\mathrm {von}}$, hydrodynamic force $p_{z}$, and central temperature $T_{\mathrm {cen}}$ increase. However, $h_{{0}}^{{*}}$ has a reverse effect on the spike of $p_{\mathrm {cen}}$. The case with a larger $h_{{0}}^{{*}}$ has a shorter impact duration and produces a thicker film layer and a larger film dimple. Consequently, less oil flows out from the contact centre to balance the pressure spike. The leaving velocities of the ball dropping from 2 and 10 mm are 0.160 and $0.409\ \mathrm{m}\cdot \mathrm{s}^{\mathrm {-1}}$, with the restitution coefficients of 0.808 and 0.924, respectively. It implies that the larger is the initial drop height, the less momentum and energy are lost in the impact–rebound process.

Figure 11 shows the distributions of pressure, film thickness, and fractional film thickness of cases 1–3 at $t = 40$, 115 and $140 \ \upmu \mathrm{s}$. Identical to the discussions of the essential features of case 2 in Sect. 4.3, a central dimple is formed in the impact process and a periphery dimple is formed in the rebound process. Both the maximum pressure and minimum film thickness occur at the contact periphery. The film with a greater central thickness is thinner at the contact periphery, and vice versa. Moreover, the difference of the pressure at the contact periphery among the three cases is similar to that of the central pressure. Since the result differences induced by the loading conditions can be revealed by the time traces, the solutions for cases with constant momentum, energy, and drop height at fixed times are not shown here.

Fig. 10

Dynamic response dependence on the impact mass: a central pressure, b central film thickness, c central temperature, d minimum film thickness, e maximum subsurface von Mises stress, and f ball velocity and hydrodynamic force

Fig. 11

Variations of the pressure and film thickness profiles in the x-direction: a $t = 40 \ \upmu \mathrm{s}$, b $t = 115 \ \upmu \mathrm{s}$, and c $t = 140 \ \upmu \mathrm{s}$

Fig. 12

Dynamic response dependence on the impact momentum: a central pressure, b central film thickness, c central temperature, d minimum film thickness, e maximum subsurface von Mises stress, and f ball velocity and hydrodynamic force

Fig. 13

Dynamic response dependence on the impact energy: a central pressure, b central film thickness, c central temperature, d minimum film thickness, e maximum subsurface von Mises stress, and f ball velocity and hydrodynamic force

4.4.2 Constant momentum

The dynamic response of cases with the same impact momentum (i.e. cases 2, 4, and 5) is compared in Fig. 12. The effects of the initial drop height $h_{{0}}^{{*}}$ on the contact pressure, film thickness, temperature, stress field, and hydrodynamic force are in accordance with the conclusions from the cases with constant impact mass. Nonetheless, the difference among the dynamic response of cases with different initial drop heights is reduced because of the change of impact mass. The leaving velocities of the ball dropping from 2 and 10 mm are 0.170 and $0.400 \mathrm{m}\cdot \mathrm{s}^{\mathrm {-1}}$, with the restitution coefficients of 0.858 and 0.903, respectively. It is demonstrated again that the larger is the initial drop height, the less momentum and energy are lost in the impact–rebound process.

4.4.3 Constant energy

In Fig. 13, the dynamic response of cases with the same impact energy (i.e. cases 2, 6, and 7) is depicted. It is consistent with the results in the above two subsections that the film thickness and the spike of the central pressure $p_{\mathrm {cen}}$ grow with the initial drop height $h_{{0}}^{{*}}$. However, the peaks of the central pressure and central temperature increase with $h_{{0}}^{{*}}$ more slowly than those of the cases with constant mass or momentum. Furthermore, the effects of $h_{{0}}^{{*}}$ on the peaks of the von Mises stress and hydrodynamic force $p_{z}$ are too small to be observed. Moreover, although the initial drop velocities of cases 2, 6, and 7 are different, their resultant restitution coefficients are all approximately 0.893. Therefore, it can be concluded that the impact energy is the decisive factor for the stress peaks, maximum hydrodynamic force, and restitution coefficient.

Based on the above analysis on the influence of the initial drop height $h_{{0}}^{{*}}$, it is found that the impact duration increases with the rise of $h_{{0}}^{{*}}$. Additionally, impact mass $m_0$ also affects the impact duration positively, as shown in Fig. 14. Therefore, the impact duration may be controlled by the coupling of $h_{{0}}^{{*}}$ and $m_0$ rather than the impact momentum or energy. Moreover, the dynamic response during the early approaching process remains unchanged with the variation of $h_{{0}}^{{*}}$, indicating that the dynamic response during the early approaching process is controlled by the initial drop height.

Fig. 14

Dynamic response dependence on the initial drop height: a central pressure, b central film thickness, c central temperature, d minimum film thickness, e maximum subsurface von Mises stress, and f ball velocity and hydrodynamic force

5 Concluding Remarks

In this paper, the material inhomogeneity is involved in the EHL problem in impact motion by means of Eshelby’s equivalent inclusion method. The developed model of the EHL problem for impact motion conditions integrates the non-Newtonian, thermal, starvation, and material inhomogeneity effects together in a strongly nonlinear equation system. By comparing the results of the homogeneous and inhomogeneous materials, it is found that the pattern of central pressure–time trace is similar to that of the pressure–distance trace at the steady-state rolling condition, and the inclusion disturbs the stress field and affects the dynamic response during the middle of the impact–rebound process when the contact pressure is high. Specifically, the impact duration and central film thickness are reduced by the stiff inclusion, while the pressure at the contact centre and the film thickness at the contact periphery are increased. However, the influences of the stiff inclusion on the ball velocity and hydrodynamic force are small enough to be neglected. Apart from the impact mass, the impact momentum and energy are concerned to investigate the influence of the initial drop height on the dynamic response of the EHL problems that were seldom discussed previously. The results show that the impact energy is the decisive factor for the stress peak, maximum hydrodynamic force, and restitution coefficient, while the dynamic response during the early approaching process is controlled by the initial drop height.

Some further research topics may be triggered by the present study: (i) the effect of randomly distributed inclusions on the contact conditions; (ii) the effect of the inlet oil layer thickness on the response of the EHL problems in impact motion conditions; (iii) the flow fields of various non-Newtonian lubricants in impact motion conditions.

Appendix

With Hooke’s law and stress superposition considered, the stresses in the equivalent homogeneous inclusions satisfy [24]

$$\begin{aligned}&C_{ijkl}^{\psi } C_{klmq}^{-1} (\sigma _{mq}^{0} +\sigma _{mq}^{p} +\sigma _{mq}^{*} )-\sigma _{ij}^{0} -\sigma _{ij}^{p} -\sigma _{ij}^{*} +C_{ijkl}^{\psi } \varepsilon _{kl}^{*} =0,\nonumber \\&\quad (\psi =1,\,2,\,\ldots ,n_{1} ;\,i,\,j,\,k,\,l,\,m,\,q=1,\,3)\, \hbox {within }\Omega _{\psi } \end{aligned}$$

(A1)

where $\sigma _{ij}^{0}$ is the eigenstress at a point within an inclusion caused by all the initial eigenstrains $\varepsilon _{ij}^{0}$, $\sigma _{ij}^{*}$ the eigenstress caused by all eigenstrains $\varepsilon _{ij}^{*}$ within the equivalent inclusions, and $\sigma _{ij}^{p}$ the stress caused by surface pressure.

According to Zhou et al. [24], the discretization expressions for $\sigma _{ij}^{0}$, $\sigma _{ij}^{*}$, and $\sigma _{ij}^{p}$ can be obtained with the assistance of the influence coefficients ${{B}}_{\alpha _{{0}}{-}\xi {,}\beta _{{0}}{-}\zeta {,}\gamma _{{0}},\varphi }$ and ${{M}}_{\alpha _{{0}}{-}\xi ,\beta _{{0}}{-}\zeta ,\gamma _{{0}}}$. ${{B}}_{\alpha _{{0}}{-}\xi {,}\beta _{{0}}{-}\zeta {,}\gamma _{{0}},\varphi }$ relates the eigenstresses $\sigma _{\alpha _{{0}},\beta _{{0}},\gamma _{{0}}}^{{0}}$ and $\sigma _{\alpha _{{0}},\beta _{{0}},\gamma _{{0}}}^{{*}}$ at the observation point in the cuboid [$\alpha _{{0}}$, $\beta _{{0}}$, $\gamma _{{0}}$] to the initial eigenstrains $\varepsilon _{\xi ,\zeta ,\varphi }^{{0}}$ and $\varepsilon _{\xi ,\zeta ,\varphi }^{{*}}$ in the cuboid [$\xi$, $\zeta$, $\varphi$], while ${{M}}_{\alpha _{{0}}{-}\xi ,\beta _{{0}}{-}\zeta ,\gamma _{{0}}}$ relates to the surface pressure $p_{\xi ,\zeta }$. Therefore, the governing equation can then be written as

$$\begin{aligned}&\left( {C_{\alpha 0,\beta _{0} ,\gamma _{0} } C^{-1}-I} \right) \left[ {\sum \limits _{\varphi =0}^{N_{Z} -1} {\sum \limits _{\varsigma =0}^{N_{Y} -1} {\sum \limits _{\xi =0}^{N_{X} -1} {{{ B}}_{\alpha _{0} -\xi ,\beta _{0} -\zeta ,\gamma _{0} ,\varphi } } } } \left( {\varepsilon _{\xi ,\zeta ,\varphi }^{P} +\varepsilon _{\xi ,\zeta ,\varphi }^{*} } \right) } \right. \nonumber \\&\quad \left. {+\sum \limits _{\varsigma =0}^{N_{Y} -1} {\sum \limits _{\xi =0}^{N_{X} -1} {{{M}}_{\alpha _{0} -\xi ,\beta _{0} -\zeta ,\gamma _{0} } } } p_{\xi ,\zeta } } \right] +C_{\alpha _{0} ,\beta _{0} ,\gamma _{0} } \varepsilon _{\xi ,\zeta ,\varphi }^{*} =0\nonumber \\&\quad \left( {0\le \alpha _{0} \le N_{x} -1,0\le \beta _{0} \le N_{y} -1,0\le \gamma _{0} \le N_{z} -1} \right) \hbox { within }\Omega _{\psi } \end{aligned}$$

(A2)

The unknown eigenstrains within the governing equation can be determined by a modified CGM with achievable convergence. The FFT algorithm is applied to improve the computational efficiency of the summation [24].

With the determined eigenstrains, the disturbed deformation $U_{\alpha _{{0}},\beta _{{0}}}^{{*}}$ at the patch [$\alpha _{{0}}, \beta _{{0}}$] can be obtained by

$$\begin{aligned} U_{\alpha _{0} ,\beta _{0} }^{*}= & {} \sum \limits _{\varphi =0}^{N_{Z} -1} {\sum \limits _{\varsigma =0}^{N_{Y} -1} {\sum \limits _{\xi =0}^{N_{X} -1} {\mathrm{{ S}}_{\alpha _{0} -\xi ,\beta _{0} -\zeta ,\varphi } } } } \left( {\varepsilon _{\xi ,\varsigma ,\varphi }^{0} +\varepsilon _{\xi ,\zeta ,\varphi }^{*} } \right) \nonumber \\&\left( {0\le \alpha _{0} \le N_{x} -1,0\le \beta _{0} \le N_{y} -1} \right) \end{aligned}$$

(A3)

where ${{S}}_{\alpha _{{0}}{-}\xi ,\beta _{{0}}{-}\zeta ,\gamma _{{0}},\varphi }$ relates the disturbed deformation $U_{\alpha _{{0}},\beta _{{0}}}^{{*}}$ due to the inclusions at point ($x_{\alpha _{{0}}}$, $y_{\beta _{{0}}})$ to the eigenstrains $\varepsilon _{\xi ,\zeta ,\varphi }^{{0}}$ and $\varepsilon _{\xi ,\zeta ,\varphi }^{{*}}$ at point $(x_{\xi }, \ y_{\zeta }, \ z_{\varphi })$.