Non-Fick mechanical-diffusion with space-dependent diffusivity and time-domain finite element method for transient impact responses analysis

Abstract

The mechanical-diffusion coupling behavior between molar concentration and strain fields have aroused great interest in biosensors, artificial muscles and actuators of adaptive structures and so on. In such nonuniform concentration environment, the space-dependent diffusivity verified in experimental observation is still not considered. Present work aims to establish the non-Fick mechanical-diffusion model with space-dependent diffusivity. To numerically solve the nonlinear governing equations, the time-domain finite element method is developed based on the principle of virtual work. The newly established model and numerical approach are applied to investigate impact responses of a thick circular plate with space-dependent diffusivity under transient chemical shock loadings. With the increase in the space-dependent diffusivity parameter, the dimensionless results reveal that the diffusive wave propagation is proportionally accelerated. And the nonlinear mechanical/chemical responses are maximally enhanced.

Introduction

Diffusion process refers to the random movement of material particles from regions of higher concentration to lower concentration, driven by concentration gradients. It is known that the diffusion of atoms can generate the chemical stresses (diffusion-induced stresses). This concept can be dated back to the experimental observation of the doping effect on the mechanical deformation of Si [1]. In the practical engineering, the coupling between diffusion and mechanics is of great significance in domains of robotic artificial skin, adaptive structural actuators, artificial muscles and so on. Mechanical deformation of solids in nonuniform concentration field has attracted significant research interests due to its extensive engineering applications in micro/nano-electromechanical systems. Notable examples include polymer electrolyte [2], viscoelastic beams [3], comb-like structure [4], thin thermoelastic plate [5], and porous micro-electrodes [6]. Currently, the interaction between molar concentration and strain fields has become a hot topic of numerous theoretical and experimental investigations in geophysics, microelectronics and environmental mechanics [7, 8]. As multifunctional materials rapidly develop, the impact of chemical stress on mechanical behavior has become increasingly important. This effect has been seen in nano-sized circular elastic matrices [9], Cu–Cr nano-layered films [10], metal matrix composites [11], and lithium-ion batteries (LIBs) [12]. To better understand the interaction between strain and concentration fields, much research has focused on theoretical models for coupled diffusion-mechanics. Early studies established the related constitutive relationships [13]. So far, there are still extensive works focused on developing coupled models that integrate diffusion and elasticity.

However, when the relevant time scale is as short as the femtosecond or picosecond, the applicability of the above-mentioned models becomes questionable. For instances, during rapid solidification of micro/nano-scale binary mixtures, ultra-fast mass transfer often occurs. At the moving solid–liquid interface, significant concentration gradients arise under highly non-equilibrium conditions due to phase transitions. This indicates that conventional equilibrium-based models are inadequate, and a non-equilibrium approach is required. Instead of the classical Fick mass diffusion law, the local non-equilibrium diffusion equation of hyperbolic type (i.e., non-Fick mass diffusion law) is more reasonable to characterize the transport feature of mass diffusion at extreme small time-scale. The non-Fick model incorporates the relaxation effect of mass diffusion, meaning that a sudden change in concentration at the boundary does not instantaneously influence distant points. This suggests that the mass diffusion occurs at a finite speed. A detailed understanding of these phenomena can be found in [14]. Currently, the non-Fick diffusion model has been widely applied in fields such as rapid solidification of metallic alloys [15] and colloidal crystallization [16]. Inspired by this, researchers have developed coupled non-Fick diffusion–elasticity theories that account for the relaxation behavior of mass transfer. These frameworks introduce a diffusion relaxation time parameter into the classical formulation, incorporating both the diffusion flux and its time derivative. As a result, the unrealistic assumption of infinite diffusion speed in both uncoupled and coupled classical models is eliminated, and the wave-like nature of mass diffusion is more accurately represented. One of the most representative models [17] is built upon the concept of inertial concentration or chemical potential [18]. Based on the above non-Fourier diffusion models, the recent progresses have been made in studying the diffusion-mechanical coupling problems with various effects, such as concentration dependency [19], nonlocal effect [20], lagging behavior [21] and fractional time dependency [22, 23]. The various studies have focused on the transient shock responses in non-Fick diffusion-mechanical systems using analytical and numerical techniques, including the Laplace transform (LT) method [24, 25], the meshless local Petrov–Galerkin (MLPG) method [26], and the mesh-free generalized finite difference (MGFD) method [27]. These techniques often involve Laplace transforms when addressing one- or two-dimensional problems. However, the process of numerical inversion introduces discretization and truncation errors. Alternatively, hybrid methods such as the Laplace transform combined with the finite element method (FEM) have also been employed for solving multi-field coupling problems [28]. When the concerned problems involve nonlinear terms (higher-order and coupled terms), the nonlinear solutions in Laplace domain are even hardly obtained. To overcome these challenges, direct time-domain finite element method [29] offers a more practical solution. This method bypasses the complexity of numerical inversion and solves both linear and nonlinear governing equations directly in the time domain.

Diffusivity can be defined as the proportionality parameter between the diffusion flux and the concentration gradient of a substance (e.g., an atom, ion or molecule). It can be used to measure how fast the mass diffuses in solids/liquids. Until now, the investigations of transient shock behavior of non-Fick mechanical-diffusion problems are mostly focused on the structure composed of space-independent material properties. However, existing theoretical and experimental studies [30, 31, 32–33] have shown that the diffusivity is characterized by spatial dependence, unlike the traditional assumption that the diffusivity of a material is uniformly distributed (constant value). Space-dependent diffusivity describes the spatial variation of a material's mass propagation properties. It can more realistically reflect the non-homogeneity that exists in a material or medium in practical applications, thus improving the prediction of concentration distribution, mass diffusion behavior, and chemical stress response. Space-dependent diffusivity mainly stems from the inherent spatial differences caused by the non-homogeneity of the material itself and the spatial correlation produced by the originally homogeneous material under specific external conditions. Space-dependent diffusivity is mainly found in non-homogeneous materials in the normal state such as composites or functional gradient materials. The diffusivity of these materials shows the spatial dependence on the macroscopic scale due to the different mass transfer properties in each region within the material. In these non-homogeneous material systems, the differences in the internal components or microstructure of the material result in inconsistent diffusivity across regions, which manifests itself as a spatial correlation. In the case of homogeneous materials such as most metallic materials, material properties such as diffusivity and thermal conductivity can also vary considerably due to natural or artificially induced inhomogeneities in the spatial system in practical problems [34]. Especially under the excitation of localized impact loadings (e.g., laser heating, chemical doping, etc.), the material structure undergoes changes (e.g., phase transitions, crystal rearrangements, carbonization, etc.), which further cause spatial variations in diffusivity. Current research on space-dependent diffusivity has focused on numerical inversion, i.e., inverting the numerical form of space-dependent diffusivities from concentration measurements, which is usually approximated as a linear or exponential dependence equation [35, 36, 37–38].

As previously discussed, the various factors of concentration-dependent diffusivity/elastic constant, stress gradient theory and spatial nonlocal/phase lagging/temporal memory dependent effects of mass transport on the non-Fick mechanical-diffusion coupling behavior have been investigated. Nevertheless, the inherent space-dependent diffusivity in the nonuniform molar concentration environment is still not considered on this topic and its influences on the transient diffusion-mechanical responses are unclear. To address the deficiencies, present work aims to establish a non-Fick mechanical-diffusion model with space-dependent diffusivity. The nonlinear time-domain finite element method is developed to directly solve the nonlinear governing equations. Following the newly established theoretical model and numerical approach, the transient shock responses of a two-dimensional thick circular plate subjected to chemical shock loadings are analyzed. Numerical results are presented and illustrated graphically. Parametric studies are conducted to evaluate the effects of space-dependent diffusivity on the nonlinear transient impact responses.

Non-Fick mechanical-diffusion model with space-dependent diffusivity

To develop a physically sound and mathematically manageable non-Fick mechanical-diffusion model with space-dependent diffusivity, the solid is assumed to be linearly elastic and homogeneous materials. The physical properties of the material (such as elastic modulus and density) are regarded as constants. In the present study, the space-dependent diffusivity is considered in the non-Fick mechanical-diffusion models [17, 24, 25], and the basically equations are given as follows:

The motion of equation and diffusion equation can be expressed as follows (in the absence of body force and diffusion source):

1

$${\sigma }_{ij,j}=\rho {\ddot{u}}_i,$$

2

$$\dot{\mu }+{\dot{\mu }}^{\left(a\right)}=-\beta J_{i,i},$$where ${\sigma }_{\mathit{ij}}$ represent the components of stress tensor, $\rho$ is the mass density, $u_i$ are the components of displacement vector. $\mu$, ${\mu }^{\left(a\right)}$ and $J_i$ are, respectively, the components of chemical potential, inertial chemical potential and diffusion flux.$$\beta =\frac{{\beta }^{\prime }}{c_0},{\beta }^{\prime }=RT,{\mu }^{\left(a\right)}=\gamma \dot{c},$$where $\beta$, $\gamma$, $R$ and $T$ are, respectively, the chemical potential constant, proportionality coefficient, the universe gas constant and absolute temperature. $c$ is the concentration of diffusing substance, $c_0$ represents the reference concentration. And it is worth noting that the superscript dot $\left(·\right)$ denotes the differentiation with respect to time, and the subscript $\left(,ij\right)$ represents the differential with respect to the material coordinates.

The stress and chemical potential constitutive relation [24, 25]:

3

$${\sigma }_{\mathit{ij}}=C_{\mathit{ijkl}}{\epsilon }_{\mathit{kl}}-{\alpha }_{\mathit{ij}}c,$$

4

$$\mu ={\alpha }_{\mathit{ij}}{\epsilon }_{\mathit{ij}}+\beta c,$$where Eq. (3) is the stress constitutive equation for the general case of anisotropic material. ${\epsilon }_{\mathit{ij}}$, $C_{\mathit{ijkl}}$, and ${\alpha }_{\mathit{ij}}$ are, respectively, the components of strain tensor, the elastic modulus and the mechanical diffusion coefficient.

The relation between the strain and displacement is given as:

5

$${\epsilon }_{\mathit{ij}}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right).$$

The Fick’s law for isotropic and homogeneous material with space-dependent diffusivity is:

6

$$J_i=-D\left(\mathbf{x}\right)c_{,i},$$where $D\left(\mathbf{x}\right)$ represents the space-dependent diffusivity. Then, the governing equations of displacement field and concentration field can be obtained, respectively:

7

$$C_{\mathit{ijkl}}{u}_{k,lj}-{\alpha }_{\mathit{ij}}{c}_{,j}=\rho {\ddot{u}}_i,$$

8

$${\left[D\left(\mathbf{x}\right)\right]}_{,i}{c}_{,i}+D\left(\mathbf{x}\right)c_{,ii}=\frac{{\alpha }_{\mathit{ij}}}{\beta }{\dot{\epsilon }}_{\mathit{ij}}+\dot{c}+\tau \ddot{c},$$where $\tau =\gamma \left./\phantom{\gamma \beta }\right)\beta$ is the diffusion relaxation time. Note that if the coupled term in Eq. (8) is ignored, the generalized mass diffusion equation with space-dependent diffusivity reads as:

9

$$\dot{c}+\tau \ddot{c}={\left[D\left(\mathbf{x}\right)\right]}_{,i}{c}_{,i}+D\left(\mathbf{x}\right)c_{,ii},$$which can degenerate into following limiting cases:

1. The generalized mass diffusion model is obtained when the diffusivity is assumed to be constant $D_m$. In such situation, the propagation speed of concentration can be evaluated as $v^{\prime }=\sqrt{D_m\left./\phantom{D_m\tau }\right)\tau }$ [14]. Evidently, the generalized diffusion equation is similar to hyperbolic heat conduction equation, which has been elaborately discussed in detail [39, 40]. Nowadays, the hyperbolic mass transfer model has been widely applied in the theoretical analysis of colloidal crystallization [16], diffusion-limited solidification in metal alloys [41] and coupled diffusion–elasticity in binary mixtures [42] and so on.

2. The classical Fick’s mass transfer model with space-dependent diffusivity can be derived when the diffusion relaxation time is ignored.

Equation (8) implies that the concentration diffuses at a finite speed of $v=\sqrt{D\left(\mathbf{x}\right)\left./\phantom{D\left(\mathbf{x}\right)\tau }\right)\tau }$. Thus far, the governing equations of the coupled non-Fick mechanical-diffusion problems based on space-dependent diffusivity are formulated.

Time-domain finite element method

When the space-dependent diffusivity is involved, the nonlinear governing equations of displacement and concentration fields are obtained (see Eqs. 7 and 8). Due to the presence of nonlinear terms in the governing equations, the mostly used numerical methods, such as LT, MLPG and MGFD, may be not applicable. To address the nonlinear problem formulated in Sect. 2, time-domain finite element method is applied to solve the nonlinear governing equations of non-Fick mechanical-diffusion problems based on space-dependent diffusivity in time domain directly. This method avoids the tedious processes and precision losses in the applications of numerical methods of LT, MLPG and MGFD. For example, in the application of MLPG method, the local integral equations (LIEs) are initially derived by using a unit step function as the test functions in the local weak-form and then the geometry domain is subdivided into a finite number of regions or elements with a circular shape. In the next step, the radial basis functions and Laplace transformation are, respectively, utilized to the treatments of the spatial and temporal variations of field variables. Finally, the inverse Laplace inverse Laplace transformation is employed to obtain the solutions in time domain.

This section is mainly devoted to establish the finite element formulations of the non-Fick mechanical-diffusion problems with space-dependent diffusivity. The constitutive Eqs. (3) and (4) can be written as the following matrix forms:

10

$$\left\{\sigma \right\}=\left[C\right]\left\{\epsilon \right\}-\left\{\alpha \right\}\left\{c\right\},$$

11

$$\left\{\mu \right\}={\left\{\alpha \right\}}^{\text{T}}\left\{\epsilon \right\}+\left\{\beta \right\}\left\{c\right\}.$$

The matrix form of Fick law of mass diffusion equation is given as:

12

$$\left\{J\right\}=-\left[D\left(\mathbf{x}\right)\right]\left\{p\right\},$$where $p=c_{,i}$. For a given element, the displacement $\left\{u\right\}$ and concentration $\left\{c\right\}$ can be expressed with two sets of shape functions $\left[N_1^e\right]$ and $\left[N_2^e\right]$:

13

$$\left\{u\right\}=\left[N_1^e\right]\left\{u^e\right\},\left\{c\right\}={\left[N_2^e\right]}^{\text{T}}\left\{c^e\right\},$$

14

$$\begin{array}{c}\left[N_1^e\right]=\left[\begin{array}{ccccccc}{N}_1& 0& N_2& 0& \cdots & N_n& 0\\ 0& N_1& 0& N_2& 0& \cdots & N_n\end{array}\right],\\ {\left[N_2^e\right]}^T=\left\{\begin{array}{cccc}{N}_1& N_2& \cdots & N_n\end{array}\right\}.\\ \end{array}$$

In Eq. (14), $n$ is the number of nodes. In terms of the nodal values of displacement $\left\{u^e\right\}$ and concentration $\left\{c^e\right\}$, the elastic strain $\left\{\epsilon \right\}$ and concentration gradient $\left\{p\right\}$ can be, respectively, expressed as:

15

$$\left\{\epsilon \right\}=\left[B_1\right]\left\{u^e\right\},\left\{p\right\}=\left[B_2\right]\left\{c^e\right\},$$where $\left[B_1\right]$ and $\left[B_2\right]$ are the strain matrix and concentration gradient matrix. The variational form of Eq. (15) can be written as:

16

$$\mathrm{\delta }\left\{\epsilon \right\}=\left[B_1\right]\mathrm{\delta }\left\{u^e\right\},\mathrm{\delta }\left\{p\right\}=\left[B_2\right]\mathrm{\delta }\left\{c^e\right\}.$$

The principle of virtual work for the non-Fick mechanical-diffusion problems with space-dependent diffusivity has the following form [29]:

17

$$\begin{array}{c}{\int }_V\left[\delta {\left\{\epsilon \right\}}^T\left\{\sigma \right\}+\left(\left\{\dot{\mu }\right\}+{\left\{\dot{\mu }\right\}}^{\left(a\right)}\right)\delta \left\{c\right\}+\delta {\left\{p\right\}}^T\left\{J\right\}\right]\text{d}V\\ =-{\int }_V\delta {\left\{u\right\}}^T\rho \left\{\ddot{u}\right\}dV+{\int }_{A_{\sigma }}\delta {\left\{u\right\}}^T\left\{\tilde{f}\right\}dA+{\int }_{A_c}\delta {\left\{c\right\}}^T\left\{\tilde{J}\right\}\text{d}A,\\ \end{array}$$where $\left\{\tilde{f}\right\}$ is the traction vector on surface $A_{\sigma }$, and $\left\{\tilde{J}\right\}$ is the mass flux vector on surface $A_c$. Substitution of Eqs. (10–16) into Eq. (17), the finite element governing equations of the non-Fick mechanical-diffusion problems with space-dependent diffusivity have the form:

18

$$\sum _{i=1}^{\text{ne}}\left({\mathbf{M}}^e{\ddot{X}}^e+{\mathbf{C}}^e{\dot{\mathbf{X}}}^e+{\mathbf{K}}^e{\mathbf{X}}^e={\mathbf{F}}_{\text{external}}^e\right),$$where $\text{ne}$ is the total number of elements. ${\mathbf{M}}^e$, ${\mathbf{C}}^{\text{e}}$, ${\mathbf{K}}^{\text{e}}$, ${\mathbf{X}}^{\text{e}}$ and ${\mathbf{F}}_{\text{external}}^{\text{e}}$ are the mass matrices, the damping matrices, the stiffness matrices, the vector of unknown constitutive variables and the external force vector, respectively, which can be represented as:$${\mathbf{M}}^{\text{e}}=\left[\begin{array}{cc}{M}_{\text{mm}}^e& 0\\ 0& M_{\text{cc}}^e\end{array}\right],{\mathbf{C}}^{\text{e}}=\left[\begin{array}{cc}\text{0}& 0\\ C_{\text{cm}}^e& C_{\text{cc}}^e\end{array}\right],{\mathbf{K}}^{\text{e}}=\left[\begin{array}{cc}{K}_{\text{mm}}^e& -K_{\text{mc}}^e\\ 0& K_{\text{cc}}^e\end{array}\right],$$$${\mathbf{X}}^{\text{e}}=\left\{\begin{array}{c}{u}^e\\ c^e\\ \end{array}\right\},{\mathbf{F}}_{\text{external}}^{\text{e}}=\left\{\begin{array}{c}{f}_m^e\\ c_c^e\\ \end{array}\right\},$$$$\left[M_{\text{mm}}^e\right]={\int }_V{\left[N_1^e\right]}^{\text{T}}\rho \left[N_1^e\right]\text{d}V,\left[M_{\text{cc}}^e\right]={\int }_V\tau \left[N_2^e\right]\left\{\beta \right\}{\left[N_2^e\right]}^{\text{T}}\text{d}V,$$$$\left[C_{\text{cm}}^e\right]={\int }_V\left[N_2^e\right]{\left\{\alpha \right\}}^{\text{T}}\left[B_1\right]\text{d}V,\left[C_{\text{cc}}^e\right]={\int }_V\left[N_2^e\right]\left\{\beta \right\}{\left[N_2^e\right]}^{\text{T}}\text{d}V,$$$$\left[K_{\text{mm}}^e\right]={\int }_V{\left[B_1\right]}^{\text{T}}\left[C\right]\left[B_1\right]\text{d}V,\left[K_{\text{mc}}^e\right]={\int }_V{\left[B_1\right]}^{\text{T}}\left\{\alpha \right\}{\left[N_2^e\right]}^{\text{T}}\text{d}V,$$$$\left[K_{\mathit{cc}}^e\right]={\int }_V{\left[B_2\right]}^{\text{T}}\left[D\left(\mathbf{x}\right)\right]\left[B_2\right]\text{d}V,\left\{f_m^e\right\}={\int }_{A_{\sigma }}{\left[N_1^e\right]}^{\text{T}}\left\{\tilde{f}\right\}\text{d}A,$$$$\left\{c_c^e\right\}={\int }_{A_{\eta }}\left[N_2^e\right]\left\{\tilde{J}\right\}\text{d}A.$$

Thus far, the finite element formulations for non-Fick mechanical-diffusion problems involving space-dependent diffusivity are given. The nonlinear finite element governing Eq. (18) can be solved directly in time domain by associating with appropriate initial and boundary conditions.

To verify the validity of time-domain finite element method before performing numerical calculations, the dynamic mechanical-diffusion response is investigated. The time-domain finite element method is applied to compare with the Laplace transform technique [21]. The initial and boundary conditions are the same as those used in [21]. The concentration curves for dimensionless time $\overline{t}=0.05$ are represented graphically. As illustrated in Fig. 1, the high degree of consistency is observed. This ensures the validity and accuracy of the method.

[See PDF for image]

Fig. 1

Verification of time-domain finite element method

Numerical examples: results and discussion

In this section, nonlinear transient responses of non-Fick mechanical-diffusion problems based on space-dependent diffusivity for an isotropic homogeneous elastic solid is studied. Hence, the elastic modulus and material coefficient ${\alpha }_{\mathit{ij}}$ have the following forms:

19

$$C_{\mathit{ijkl}}=\lambda {\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}}+\upsilon {\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+\upsilon {\delta }_{\mathit{il}}{\delta }_{\mathit{jk}},{\alpha }_{\mathit{ij}}=\alpha {\delta }_{\mathit{ij}},$$where $\lambda$ and $\upsilon$ are Lame’s constants. The governing equations can be written as:

20

$${\sigma }_{ij,j}=\rho {\ddot{u}}_i,\dot{\mu }+{\dot{\mu }}^{\left(a\right)}=-\beta J_{i,i},$$

21

$${\sigma }_{\mathit{ij}}=\lambda {\epsilon }_{\mathit{kk}}{\delta }_{\mathit{ij}}+2\upsilon {\epsilon }_{\mathit{ij}}-\alpha c{\delta }_{\mathit{ij}},\mu =\alpha {\epsilon }_{\mathit{kk}}+\beta c,$$

22

$${\epsilon }_{\mathit{ij}}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right),J_i=-D\left(\mathbf{x}\right)c_{,i}.$$

Substitution of Eqs. (21)₁ and (22)₁ into Eq. (20)₁ leads to:

23

$$\upsilon u_{i,jj}+\left(\lambda +\upsilon \right)u_{j,ij}-\alpha c_{,i}=\rho {\ddot{u}}_i.$$

Substitution of Eqs. (21)₂ and (22)₂ into (20)₂ yields:

24

$${\left[D\left(\mathbf{x}\right)\right]}_{,i}{c}_{,i}+D\left(\mathbf{x}\right)c_{,ii}=\frac{\alpha }{\beta }{\dot{\epsilon }}_{\mathit{kk}}+\dot{c}+\tau \ddot{c}.$$

In the governing equations of motion and concentration (i.e., Eqs. 23 and 24), it needs to be pointed out that the space-dependent diffusivity (i.e., $D\left(\mathbf{x}\right)$) is assumed as the following form in the two-dimensional case [31]:

25

$$D\left(\mathbf{x}\right)=D\left(x,y\right)=D_m\left[1+{\chi }_1\left(\frac{x}{x_0}+\frac{y}{y_0}\right)+{\chi }_2\left(\frac{x^2}{x_0^2}+\frac{\mathit{xy}}{x_0{y}_0}+\frac{y^2}{y_0^2}\right)\right],$$where $D_m$ is reference constant diffusivity. The parameters ${\chi }_1$ and ${\chi }_2$ are the parameters describing the variation of the diffusivity.

For numerical evaluations, the material parameters are shown in Table 1.

Table 1. Material constant in numerical simulation [21]

$\alpha$	$\phi$	$g$	$\zeta$	$E$
$1.871$	$1.978\times {10}^{-4}$	$0.583$	$0.3$	$2.3\times {10}^9\text{Pa}$

Verification

To check the stability of the established non-Fick mechanical-diffusion model, the transient shock behavior of a thin plate in the absence of space-dependent material properties is investigated. In such case, the parameters ${\chi }_1={\chi }_2=0$. The schematic model of a thin plate with the thickness of $h$ is displayed in Fig. 2. The medium is initially quiescent. At the upper surface ($z=0$), the thin plate is under the transient chemical shock loadings in the form of a Heaviside unit step function of time. It is assumed that the dimension of x-axis and y-axis is much larger than the thickness of the plate (i.e., z-axis). During the analysis, loading time is limited so that neither elastic wave nor diffusive wave reaches to the lower surface of the plate ($z=h$). Hence, the problem considered can be treated as a one-dimensional case. All the field variables depend only on $z$ and $t$. The components of the displacement and concentration are:

26

$$u_x=0,u_y=0,u_z=w\left(z,t\right),c=c\left(z,t\right).$$

[See PDF for image]

Fig. 2

The schematic model of the thin plate

Initial conditions are:

27

$$w\left(z,0\right)=\frac{\partial w\left(z,0\right)}{\partial t}=0,c\left(z,0\right)=\frac{\partial c\left(z,0\right)}{\partial t}=0\text{.}$$

Boundary conditions are:

28

$${\sigma }_{\mathit{zz}}\left(0,t\right)=0,c\left(0,t\right)=c_{0}H\left(t\right),$$

29

$$u\left(h,t\right)=0,c\left(h,t\right)=0.$$

Then, the governing Eqs. (20–22) have the forms:

30

$$\frac{\partial {\sigma }_{\mathit{zz}}}{\partial z}=\rho \frac{{\partial }^{2}w}{\partial t^2},\frac{\partial \mu }{\partial t}+\frac{\partial {\mu }^{\left(a\right)}}{\partial t}=-\beta \frac{\partial J_z}{\partial z},$$

31

$${\sigma }_{\mathit{zz}}=\left(\lambda +2\upsilon \right){\epsilon }_{\mathit{zz}}-\alpha c,{\sigma }_{\mathit{xx}}={\sigma }_{\mathit{yy}}=\lambda {\epsilon }_{\mathit{zz}}-\alpha c,$$

32

$$\mu =\alpha {\epsilon }_{\mathit{zz}}+\beta c,$$

33

$${\epsilon }_{\mathit{zz}}=\frac{\partial w}{\partial z},J_z=-D_m\frac{\partial c}{\partial z}.$$

For convenience, the following dimensionless quantities are introduced:

34

$$\begin{array}{c}\left({\overline{x}}_i,{\overline{u}}_i\right)=\sqrt{\frac{\lambda +2\upsilon }{\rho }}\frac{1}{D_m}\left(x_i,u_i\right),\\ \left(\overline{t},\overline{\tau }\right)=\frac{\lambda +2\upsilon }{\rho }\frac{1}{D_m}\left(t,\tau \right),{\overline{\sigma }}_{\mathit{ij}}=\frac{{\sigma }_{\mathit{ij}}}{\upsilon },\overline{c}=\frac{c}{c_0}.\\ \end{array}$$

In terms of dimensionless quantities, the dimensionless governing equations of motion and concentration can be written as:

35

$${\varsigma }^2\frac{{\partial }^2\overline{w}}{\partial {\overline{z}}^2}-\phi \frac{\partial \overline{c}}{\partial \overline{z}}={\varsigma }^2\frac{{\partial }^2\overline{w}}{\partial {\overline{t}}^2},$$

36

$$\frac{{\partial }^2\overline{c}}{\partial {\overline{z}}^2}=g\frac{{\partial }^2\overline{w}}{\partial \overline{z}\partial \overline{t}}+\frac{\partial \overline{c}}{\partial \overline{t}}+\overline{\tau }\frac{{\partial }^2\overline{c}}{\partial {\overline{t}}^2},$$where ${\varsigma }^2=\frac{\lambda +2\upsilon }{\upsilon }$, $\phi =\frac{\alpha c_0}{\upsilon }$ and $g=\frac{\alpha }{\beta c_0}$.

The governing Eqs. (35) and (36) can be transformed into Laplace domain and solved, then the results in time domain are obtained via the inverse Laplace transform [24, 25]. Nevertheless, numerical evaluation of inverse transform is time-consuming and the high-resolution results for transient shock problems are difficult to be captured. As a result, some delicate features expected in non-Fick mechanical-diffusion problems (e.g., the finite propagation speed of concentration) could not be faithfully predicted by using the Laplace transformation method. In terms of wave Eq. (36), the dimensionless velocity of mass diffusion can be calculated as ${\overline{\nu }}_{\text{molar concentration}}^{\prime }=\sqrt{1\left./\phantom{1\overline{\tau }}\right)\overline{\tau }}$. When the dimensionless diffusion relaxation time is chosen as $\overline{\tau }=0.05$, the propagation velocity of the diffusive wave is 4.472 and the jump of concentration at the wave front of concentration should be at about ${\overline{z}}_{\text{molar concentration}}^{\prime }=0.224$ at $\overline{t}=0.05$. This suggests that the response of concentration vanishes identically when $\overline{z}>0.224$. However, the value of concentration approaches to zero until when $\overline{z}>0.53$ in available literature [25]. Evidently, the finite propagation feature of concentration is not depicted exactly. To capture the transient responses of the problem more precisely, time-domain finite element method is used to solve the above one-dimensional governing Eqs. (35) and (36) in time domain directly associated with the initial and boundary conditions (27–29). In the calculation, the dimensionless thickness of the thin plate is selected as $\overline{h}=1.5$. The dimensionless concentration and displacement are illustrated graphically in Figs. 3 and 4, respectively.

[See PDF for image]

Fig. 3

The distribution of concentration at different times (a) and diffusion relaxation times (b)

[See PDF for image]

Fig. 4

The distribution of displacement for different times with $\overline{\tau }=0.05$

As displayed in Fig. 3a, the distribution of dimensionless concentration along the thickness direction of the plate is presented for the instants $\overline{t}=0.05$, 0.08 and 0.10 when $\overline{\tau }$ is 0.05. The corresponding wave fronts of diffusive wave are $\overline{z}\left|_{\overline{t}=0.05}^{}\right)=0.224$, $\overline{z}\left|_{\overline{t}=0.08}^{}\right)=0.358$ and $\overline{z}\left|_{\overline{t}=0.10}^{}\right)=0.447$ (see Eq. (36)). It is observed that the evaluated wave fronts of concentration from finite element analysis are slightly off that obtained from theoretical predictions. This implies that the finite element results are reliable. The feature is also clearly identified from the curves shown in Fig. 3a, namely, the sharp jumps of concentration are around at the diffusive wave fronts. It is found that the response region of concentration is essentially undisturbed when $\overline{z}>0.224/0.358/0.447$. On the other hand, it can be concluded that concentration diffuses at a finite speed in the generalized diffusion–elasticity theory, instead of an infinite velocity as implied in the classical theory of diffusion–elasticity. From Fig. 3a, it is also seen that the magnitude of the concentration jump reduces dynamically with the passage of time. It means that the jump will vanish after a long time.

To evaluate the effect of diffusion relaxation time on the distribution of concentration, Fig. 3b displays the distribution of concentration for $\overline{\tau }=0.02$, 0.05 and 0.08 at $\overline{t}=0.05$. The corresponding wave fronts are ${\overline{z}}^{\prime }\left|_{\overline{\tau }=0.02}^{}\right)=0.353$, ${\overline{z}}^{\prime }\left|_{\overline{\tau }=0.05}^{}\right)=0.224$ and ${\overline{z}}^{\prime }\left|_{\overline{\tau }=0.08}^{}\right)=0.177$. As shown in Fig. 3b, it is also obtained that the diffusive wave travels slower when the diffusion relaxation time is larger. Figure 4 presents the variations of displacement along the thickness direction of the plate at $\overline{t}=0.05$, 0.08 and 0.10 when the diffusion relaxation time $\overline{\tau }$ is set as 0.05. Because the lower surface of the thin plate is constrained, the upper surface of the thin plate moves upward (negative displacement) when it is exposed to the shock loading of concentration. From Figs. 3 and 4, it can be concluded that with the propagation of elastic wave and diffusive wave, the magnitudes of displacement at the upper surface of the thin plate increases steadily with the absorption of energy. Nevertheless, the magnitude of the jump of concentration is decreasing with the dissipation of energy. This is consistent with the essence that mass diffuses at a finite velocity in the generalized diffusion–elasticity theory.

Impact responses of transient chemical shock affected thick circular plate with space-dependent diffusivity

Formulation of the problem

To investigate the influences of space-dependent diffusivity on the transient response of non-Fick mechanical-diffusion problems, a thick circular plate of thickness of $h$ with space-dependent diffusivity is taken as an example (see Fig. 5). The plate occupies the region $\mathrm{\Omega }=\left\{\left(r,\theta ,z\right)\left|0\le r<\infty ,0\le \theta \le 2\pi ,0\le z\le h\right)\right\}$. The cylindrical system of coordinates $\left(r,\theta ,z\right)$ with the z-axis coinciding with the axis of the cylinder is chosen. The axis of symmetry is selected as z-axis and the origin of the cylindrical coordinates systems is fixed at the upper surface of the plate. The initial state of the plate is assumed to be quiescent. The upper surface of the plate is taken to be traction free and its central parts are subjected to a zonal time-dependent concentration disturbance. The issue can be treated as an axisymmetric problem. The components of displacement and concentration are simplified as:

37

$$u_r=u_r\left(r,z,t\right),u_{\theta }=0,u_z=w\left(r,z,t\right),$$

38

$$c=c\left(r,z,t\right).$$

[See PDF for image]

Fig. 5

The schematic model of the thick circular plate

The initial and boundary conditions can be expressed as follows:

Initial conditions are:

39

$$u_r\left(r,z,0\right)=w\left(r,z,0\right)=0,\frac{\partial u_r\left(r,z,0\right)}{\partial t}=\frac{\partial w\left(r,z,0\right)}{\partial t}=0,$$

40

$$c\left(r,z,0\right)=0,\frac{\partial c\left(r,z,0\right)}{\partial t}=0.$$

Boundary conditions are:

41

$${\sigma }_{\mathit{zz}}\left(r,0,t\right)={\sigma }_{\mathit{rz}}\left(r,0,t\right)=0,$$

42

$$u_r\left(r,h,t\right)=w\left(r,h,t\right)=0,$$

43

$$c\left(r,0,t\right)=c_{0}H\left(t\right)H\left(a-\left|r\right|\right).$$

The components of strain in Eq. (22)₁ can be written as:

44

$${\epsilon }_{\mathit{rr}}=\frac{\partial u_r}{\partial r},{\epsilon }_{\theta \theta }=\frac{u_r}{r},{\epsilon }_{\mathit{zz}}=\frac{\partial w}{\partial z},{\epsilon }_{\mathit{rz}}=\frac{\partial u_r}{\partial z}+\frac{\partial w}{\partial r},{\epsilon }_{r\theta }={\epsilon }_{\theta z}=0.$$

Thus, Eqs. (20) and (21) have the forms as:

45

$$\frac{\partial {\sigma }_{\mathit{rr}}}{\partial r}+\frac{\partial {\sigma }_{\mathit{rz}}}{\partial z}+\frac{{\sigma }_{\mathit{rr}}-{\sigma }_{\theta \theta }}{r}=\rho \frac{{\partial }^2{u}_r}{\partial t^2},\frac{\partial {\sigma }_{\mathit{zz}}}{\partial z}+\frac{\partial {\sigma }_{\mathit{rz}}}{\partial r}+\frac{{\sigma }_{\mathit{rz}}}{r}=\rho \frac{{\partial }^{2}w}{\partial t^2},$$

46

$${\sigma }_{\mathit{rr}}=\lambda e+2\upsilon {\epsilon }_{\mathit{rr}}-\alpha c,{\sigma }_{\theta \theta }=\lambda e+2\upsilon {\epsilon }_{\theta \theta }-\alpha c,$$

47

$${\sigma }_{\mathit{zz}}=\lambda e+2\upsilon {\epsilon }_{\mathit{zz}}-\alpha c,{\sigma }_{\mathit{rz}}=\upsilon {\epsilon }_{\mathit{rz}},{\sigma }_{r\theta }={\sigma }_{z\theta }=0,$$

48

$$\mu =\alpha e+\beta c.$$

Then, the governing equations of motion and concentration are obtained as follows:

49

$$\upsilon {\mathrm{\nabla }}^2{u}_r-\upsilon \frac{u_r}{r^2}+\left(\lambda +\upsilon \right)\frac{\partial e}{\partial r}-\alpha \frac{\partial c}{\partial r}=\rho \frac{{\partial }^2{u}_r}{\partial t^2},$$

50

$$\upsilon {\mathrm{\nabla }}^{2}w+\left(\lambda +\upsilon \right)\frac{\partial e}{\partial z}-\alpha \frac{\partial c}{\partial z}=\rho \frac{{\partial }^{2}w}{\partial t^2},$$

51

$$D\left(r,z\right){\mathrm{\nabla }}^{2}c+\left[\frac{\partial D\left(r,z\right)}{\partial r}+\frac{\partial D\left(r,z\right)}{\partial z}\right]\left(\frac{\partial c}{\partial r}+\frac{\partial c}{\partial z}\right)=\frac{\alpha }{\beta }\frac{\partial e}{\partial t}+\frac{\partial c}{\partial t}+\tau \frac{{\partial }^{2}c}{\partial t^2},$$where $e$ is the cubical dilatation:

52

$$e=\frac{1}{r}\frac{\partial }{\partial r}\left(u_{r}r\right)+\frac{\partial w}{\partial z}=\frac{u_r}{r}+\frac{\partial u_r}{\partial r}+\frac{\partial w}{\partial z}.$$

The Laplace operator is given by

53

$${\mathrm{\nabla }}^2=\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{{\partial }^2}{\partial z^2}.$$

In terms of dimensionless quantities (see Eq. 34), Eqs. (49–51) can be written as:

54

$${\mathrm{\nabla }}^2{\overline{u}}_r-\frac{{\overline{u}}_r}{{\overline{r}}^2}+\left({\varsigma }^2-1\right)\frac{\partial \overline{e}}{\partial \overline{r}}-\phi \frac{\partial \overline{c}}{\partial \overline{r}}={\varsigma }^2\frac{{\partial }^2{\overline{u}}_r}{\partial {\overline{t}}^2},$$

55

$${\mathrm{\nabla }}^2\overline{w}+\left({\varsigma }^2-1\right)\frac{\partial \overline{e}}{\partial \overline{z}}-\phi \frac{\partial \overline{c}}{\partial \overline{z}}={\varsigma }^2\frac{{\partial }^2\overline{w}}{\partial {\overline{t}}^2},$$

56

$${\xi }_1\left(\overline{r},\overline{z}\right){\mathrm{\nabla }}^2\overline{c}+{\xi }_2\left(\overline{r},\overline{z}\right)\left(\frac{\partial \overline{c}}{\partial \overline{r}}+\frac{\partial \overline{c}}{\partial \overline{z}}\right)=g\frac{\partial \overline{e}}{\partial \overline{t}}+\frac{\partial \overline{c}}{\partial \overline{t}}+\overline{\tau }\frac{{\partial }^2\overline{c}}{\partial {\overline{t}}^2},$$where$$\begin{array}{c}{\xi }_1\left(\overline{r},\overline{z}\right)=\left[1+{\chi }_1\left(\frac{\overline{r}}{{\overline{r}}_0}+\frac{\overline{z}}{{\overline{z}}_0}\right)+{\chi }_2\left(\frac{{\overline{r}}^2}{{\overline{r}}_0^2}+\frac{\overline{r}\overline{z}}{{\overline{r}}_0{\overline{z}}_0}+\frac{{\overline{z}}^2}{{\overline{z}}_0^2}\right)\right],\\ {\xi }_2\left(\overline{r},\overline{z}\right)=\left[{\chi }_1\left(\frac{1}{{\overline{r}}_0}+\frac{1}{{\overline{z}}_0}\right)+{\chi }_2\left(\frac{2\overline{r}}{{\overline{r}}_0^2}+\frac{\overline{r}+\overline{z}}{{\overline{r}}_0{\overline{z}}_0}+\frac{2\overline{z}}{{\overline{z}}_0^2}\right)\right].\\ \end{array}$$

By using time-domain finite element method, the two-dimensional nonlinear governing Eqs. (54–56) can be solved in time domain directly by associating with the initial and boundary conditions (39–43). In the calculation, the finite element analysis model is sketched as shown in Fig. 6. The constants used for numerical calculation are taken as:$$\overline{t}=0.05,\overline{\tau }=0.02,\overline{a}=0.5,\overline{h}=1.5.$$

[See PDF for image]

Fig. 6

Finite element analysis model

The effects of space-dependent diffusivity parameters ${\chi }_1$

To investigate the influences of different space-dependent diffusivity parameters ${\chi }_1$ on the nonlinear transient responses, four cases are discussed (Table 2): in which ${\chi }_1>0$ (i.e., ${\chi }_1=0.2$, 0.5 and 0.8) represents the variation of the first-order space-dependent diffusivity.

Table 2. Simulation cases with different parameters ${\chi }_1$

Case	1	2	3	4
${\chi }_1$	0.0	0.2	0.5	0.8
${\chi }_2$	0.0	0.0	0.0	0.0

Figure 7 displays the evolution of the dimensionless concentration contours $\overline{c}$ in the thick circular plate under different values of ${\chi }_1$, with ${\chi }_2=0$. It is observed that as ${\chi }_1$ increases from 0.0 to 0.8, the response region of concentration is limited in a finite area and the distribution of concentration is almost unchanged beyond this area. This means that the mass diffuses at a finite speed due to the relaxation time $\overline{\tau }$. The region of the concentration field expands in both the $\overline{r}$ and $\overline{z}$ directions, indicates that the enhanced diffusivity gradient promotes broader mass transport. This trend lies in the form of the space-dependent diffusivity function $D\left(\overline{r},\overline{z}\right)$, which increases with ${\chi }_1$. This implies that the diffusivity of material increases radially and axially with distance from the origin. To further discuss the influence of space-dependent diffusivity on the concentration, the distributions of concentration along $\overline{z}$-axis and $\overline{r}$-axis are presented in Fig. 8a, b, respectively. As shown in Fig. 8a, the jump of concentration along z-axis at the diffusive wave front moves forward with the increases in ${\chi }_1$. Reviewing the governing equation of concentration (56), the propagation velocity of concentration and the corresponding wave front can be evaluated as:

57

$${\overline{v}}_{\text{Diffusive wave}}=\sqrt{{\xi }_1\left(\overline{r},\overline{z}\right)\left./\phantom{{\xi }_1\left(\overline{r},\overline{z}\right)\overline{\tau }}\right)\overline{\tau }},{\overline{z}}_{\text{Diffusive wave front}}=\sqrt{{\xi }_1\left(\overline{r},\overline{z}\right)\left./\phantom{{\xi }_1\left(\overline{r},\overline{z}\right)\overline{\tau }}\right)\overline{\tau }}\times \overline{t},$$which indicates that the propagation velocity of concentration depends on the space distribution of diffusivity. The diffusive wave front shifts from $\overline{z}\approx 0.35$ to $\overline{z}\approx 0.45$ with ${\chi }_1$ increases from 0.0 to 0.8. And the value of $\overline{c}$ at the same position increases due to the enlargement of ${\chi }_1$. This signifies that the enhanced ${\chi }_1$ has decreased the concentration gradient over a broader region. Figure 8b further supports this observation by showing the variation of $\overline{c}$ along the $\overline{r}$-axis. The diffusive wave front along the $\overline{r}$-axis moves from $\overline{r}\approx 0.85$ to $\overline{r}\approx 1.0$. In addition, the concentration keeps constant within $0\le \overline{r}\le 0.5$, which is consistent with the boundary conditions at the upper surface of the thick circular plate. The outer layers of the plate allow faster mass transport due to the increase in space-dependent diffusivity.

[See PDF for image]

Fig. 7

The contours of concentration in the thick circular plate

[See PDF for image]

Fig. 8

The distribution of concentration along a$\overline{z}$-axis b$\overline{r}$-axis

Figure 9 depicts the dimensionless axial stress contours ${\overline{\sigma }}_{\mathit{zz}}$ in the thick circular plate under different values of the parameter ${\chi }_1$, with ${\chi }_2=0$. As shown in the figures, the elastic wave propagates at a finite speed and the stress response is confined within a finite region. The compressive stress primarily concentrated in the center of the high-stress zone. As ${\chi }_1$ increases, significant changes are observed in both the stress distribution and the extent of the affected region. Specifically, the stress field exhibits a sharp localized concentration near the chemically loaded surface (around $\overline{z}=0.1$), which forms a pronounced compressive stress zone. As ${\chi }_1$ increases to 0.80, the high-stress region expands gradually along both the radial and axial directions. While ${\chi }_1$ increases from 0.0 to 0.8, the axial stress response region shifts from $\overline{z}=0.4$ to $\overline{z}=0.5$. The radial stress response region enhances from $\overline{r}\approx 0.9$ to $\overline{r}\approx 1.0$. The contour patterns become smoother, which indicates a broader stress diffusion region and the stress gradient gradually decreases. Notably, under ${\chi }_1=0.5$ and 0.8, the high-stress zone expands outward. It indicates that the increase of parameter ${\chi }_1$ under chemical impact loading enhances the stress propagation in elastic media. Furthermore, the distributions of ${\overline{\sigma }}_{\mathit{zz}}$ along the axial $\overline{z}$-direction and ${\overline{\sigma }}_{\mathit{rr}}$ along the radial r-direction are shown in Fig. 10a, b, respectively. Figure 10a demonstrates that the peak value of ${\overline{\sigma }}_{\mathit{zz}}$ initially reaches of $-1.6\times {10}^{-4}$ at $\overline{z}=0.1$. As ${\chi }_1$ improves, the absolute value of the peak compressive stress also increases and the value of ${\overline{\sigma }}_{\mathit{zz}}$ rises within $0.1\le \overline{r}\le 0.5$. The stress response range expands rightward along the z-axis, reflecting an enhancement in wave transmission due to space-dependent diffusivity. Figure 10b presents the radial stress ${\overline{\sigma }}_{\mathit{rr}}$ distribution under varying ${\chi }_1$. Due to boundary constraints, ${\overline{\sigma }}_{\mathit{rr}}$ shows no significant oscillation within the range $0\le \overline{r}\le 0.5$. In summary, the enhancement of ${\chi }_1$ improves the stress propagation and mitigates the stress concentration.

[See PDF for image]

Fig. 9

The contours of stress in the thick circular plate

[See PDF for image]

Fig. 10

The distribution of stress along a$\overline{z}$-axis b$\overline{r}$-axis

Figure 11a, b show the distributions of axial displacement $\overline{w}$ along the $\overline{z}$-axis and radial displacement ${\overline{u}}_r$ along the $\overline{r}$-axis, respectively. In Fig. 11a, the axial displacement curve exhibits a distinct tensile peak value ($2.5\times {10}^{-7}$) near $\overline{z}=0.05$ prior to the arrival of the elastic wave. This peak pattern corresponds to the arrival and propagation of the elastic wave triggered by the chemical load. With the increase of ${\chi }_1$, the position of the elastic wave front remains unchanged. This indicates that the elastic wave propagation speed is not affected by the space-dependent diffusivity, which agrees with the numerical results. The value peak of displacement slightly decreases, and the slope of the descending segment becomes gentler, resulting in a smoother curve. This suggests that stronger space-dependent diffusivity can suppress excessive local displacement induced by the chemical stress, thereby making the axial deformation of the structure more uniform and stable. Figure 11b presents a sharp peak of ${\overline{u}}_r$ approximately $3.5\times {10}^{-7}$ occurs at $\overline{r}\approx 0.6$. As ${\chi }_1$ increases, the amplitude of peak slightly decreases and the curve transitions become more gradual. This trend reflects the coupling mechanism of the mass diffusion field on the mechanical response. The enhancement of ${\chi }_1$ leads to a more uniform stress distribution, which causes a lower displacement peak value. Although the space-dependent diffusivity does not directly change the elastic wave velocity, it significantly influences the magnitude of structural deformation. This can be interpreted as a damping effect induced by space-dependent diffusivity under chemical impact loading.

[See PDF for image]

Fig. 11

The distribution of displacement along a$\overline{z}$-axis b$\overline{r}$-axis

The effects of space-dependent diffusivity parameters ${\chi }_2$

To investigate the influences of different space-dependent diffusivity parameters ${\chi }_2$ on the nonlinear transient responses, another four cases are investigated (Table 3): in which ${\chi }_2>0$ (i.e., ${\chi }_2=0.5$ and 0.8) represents the variation of higher-order space-dependent diffusivity. The value of space-dependent diffusivity is increased when the higher-order parameter ${\chi }_2$ is considered. This suggests the capability of mass transfer is enhanced.

Table 3. Simulation cases with different parameters ${\chi }_2$

Case	5	6	7	8
${\chi }_1$	0.0	0.2	0.2	0.2
${\chi }_2$	0.0	0.2	0.5	0.8

Figure 12 presents the two-dimensional contours of dimensionless concentration $\overline{c}$ in the thick circular plate, where the space-dependent diffusivity parameter ${\chi }_2$ varies from 0 to 0.8, while ${\chi }_1=0.2$ remains fixed. It is observed that with the increase of ${\chi }_2$, the region of elevated concentration expands significantly along both the axial and radial directions. For instance, the outermost contour line shifts from $\overline{z}\approx 0.4$ (at ${\chi }_2=0$) to $\overline{z}\approx 0.5$ (at ${\chi }_2=0.8$), which reflects a faster diffusive wave front advancement. This trend is attributed to the enhanced nonlinearity in the diffusivity function $D\left(\overline{r},\overline{z}\right)$ due to the higher-order term parameter ${\chi }_2$, which amplifies the spatial variation of the diffusivity. Figure 13a, b present the distributions of concentration along $\overline{z}$-axis and $\overline{r}$-axis, respectively. Figure 13a shows the axial concentration profiles, where the diffusive wave front shifts forward and the gradient near the front decreases with increasing ${\chi }_2$. The concentration curve jump occurs around $\overline{z}\approx 0.35$ for ${\chi }_2=0$. As ${\chi }_2=0.8$, the jump location delays to $\overline{z}\approx 0.48$. This reflects a higher diffusion velocity associated with increased space-dependent diffusivity. When ${\chi }_2>0$, the diffusive wave velocity is accelerated and the corresponding wave front moves forward with increasing ${\chi }_2$. When ${\chi }_2<0$, the diffusive wave velocity is slowed down and the corresponding wave front moves backward with decreasing ${\chi }_2$. Similarly, Fig. 13b demonstrates that the radial concentration curve remains relatively flat for $\overline{r}\le 0.5$, followed by a smooth decay. And the higher ${\chi }_2$ leads to a more gradual concentration decay, consistent with the smoother diffusive wave front. It can be concluded that the space-dependent effect on diffusivity is a significant factor that cannot be ignored in determining the responses of concentration in transient shock problems.

[See PDF for image]

Fig. 12

The contours of concentration in the thick circular plate

[See PDF for image]

Fig. 13

The distribution of concentration along a$\overline{z}$-axis b$\overline{r}$-axis

Figure 14 displays the contours of dimensionless axial stress ${\overline{\sigma }}_{\mathit{zz}}$ under various parameters ${\chi }_2$, with ${\chi }_1=0.2$. It is observed that as ${\chi }_2$ increases from 0 to 0.8, the region of high stress expands both radially and axially. The outer contour line moves from $\overline{z}\approx 0.4$ to $\overline{z}\approx 0.5$, which indicates a broader stress influence zone. Moreover, the spacing between adjacent contours increases, especially for ${\chi }_2>0.5$, reflecting a reduction in stress gradients and a smoother stress field. Although the high-stress region remains relatively concentrated, it gradually extends outward as ${\chi }_2$ increases. The stress concentration inside the structure is further mitigated. This reflects the fact that higher-order space-dependent diffusivity provides more effective dispersion and buffering effects on local impact responses. Figure 15a further confirms this trend by illustrating axial stress curves. The compressive stress valley becomes slightly deeper, reaching approximately ${\overline{\sigma }}_{\mathit{zz}}\approx 1.65\times {10}^{-4}$ for ${\chi }_2=0.8$. And the stress recovery slope becomes significantly more gradual compared to the sharp transition observed at ${\chi }_2=0$. This indicates improved value of stress due to the introduction of ${\chi }_2$. Figure 15b presents the radial stress ${\overline{\sigma }}_{\mathit{rr}}$. The radial stress gradually rises from ${\overline{\sigma }}_{\mathit{rr}}\approx 1.2\times {10}^{-4}$ to 0. As ${\chi }_2$ increases from 0 to 0.8, the value of ${\overline{\sigma }}_{\mathit{rr}}$ also increases within $0.6\le \overline{r}\le 1.0$. It shows a similar trend with ${\overline{\sigma }}_{\mathit{zz}}$. However, the effect of ${\chi }_2$ on radial stress is relatively weaker and mainly reflected in a slight change in the curvature of the recovery phase.

[See PDF for image]

Fig. 14

The contours of stress in the thick circular plate

[See PDF for image]

Fig. 15

The distribution of stress along a$\overline{z}$-axis b$\overline{r}$-axis

Figure 16a, b show the axial and radial distributions of the dimensionless displacement under different space-dependent diffusivity parameters ${\chi }_2$, with ${\chi }_1=0.2$. In Fig. 16a, the axial displacement curve initially reaches a tensile peak of approximately $2.6\times {10}^{-7}$ near $\overline{z}\approx 0.1$. As ${\chi }_2$ increases from 0 to 0.8, the peak amplitude slightly decreases about $3\%$-$5\%$. Similarly, Fig. 16b presents the radial displacement evolution. The maximum occurs at $\overline{r}\approx 0.5$, with the peak amplitude of roughly $3.6\times {10}^{-7}$. As ${\chi }_2$ increases, the peak magnitude reduces slightly and the curvature near the crest flattens. The influence of higher-order space-dependent diffusivity on displacement is weaker than that on stress. This is attributed to the fact that displacement is an integrated response quantity, less sensitive to local variations in diffusion compared to stress, which depends directly on the propagation of chemical concentration.

[See PDF for image]

Fig. 16

The distribution of displacement along a$\overline{z}$-axis b$\overline{r}$-axis

Conclusions

The main goal of this paper is to establish the non-Fick mechanical-diffusion model with space-dependent diffusivity to investigate the impact responses of a thick circular plate under transient chemical shock loadings by time-domain finite element method. From numerical results, the following conclusions can be reached:

1. The space-dependent diffusivity parameters ${\chi }_1$ and ${\chi }_2$ significantly influence the nonlinear transient response of the concentration field. The diffusive wave front shifts from $\overline{z}\approx 0.35$ to $\overline{z}\approx 0.45$ with ${\chi }_1$ increases from 0 to 0.8. As ${\chi }_1$ and ${\chi }_2$ increase, the value of $\overline{c}$ increases before the diffusive wave front and the velocity of diffusive wave increases accordingly.

2. The stress contours and curves demonstrate that space-dependent diffusivity modulates stress fields through mechanical-diffusion coupling. For ${\chi }_1$ increases from 0 to 0.8, the stress response region extends from $\overline{z}\approx 0.4$ to $\overline{z}\approx 0.5$ and the peak value of ${\overline{\sigma }}_{\mathit{zz}}$ increases from approximately $-1.6\times {10}^{-4}$ to $-1.7\times {10}^{-4}$. With ${\chi }_2$ increases, the contour spacing becomes significantly wider and the stress gradient is reduced.

3. Compared to stress and concentration, the displacement fields are less sensitive to space-dependent diffusivity variations. In the $\overline{z}$-axis, the maximum displacement $\overline{w}$ peak slightly reduces from $2.6\times {10}^{-7}$ to $2.5\times {10}^{-7}$ as ${\chi }_2$ increases from 0 to 0.8. Radial displacement ${\overline{u}}_r$ shows a peak at $\overline{r}\approx 0.5$ with amplitude decreasing from $3.6\times {10}^{-7}$ to $3.3\times {10}^{-7}$. This can be interpreted as a damping effect induced by space-dependent diffusivity under chemical impact loading.

Although the results provide valuable insights, the model has some limitations, and future research directions warrant attention. The current model assumes a quadratic functional space-dependent diffusivity to establish a clear theoretical framework and evaluate the intrinsic effects of space-dependent diffusivity. However, this abstraction has not yet explained all the complexities of real-world scenarios, such as space-dependent elastic modulus. And the present diffusivity function cannot adapt to all computational models. The future works can be contributed efforts to the non-Fick mechanical-diffusion modeling on such topics.

Acknowledgements

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Natural Science Foundation of China (12362014，12102158), the Natural Science Foundation of Gansu Province (21JR1RA241), the Opening Project from the State Key Laboratory for Strength and Vibration of Mechanical Structures (SV2019-KF-30, SV2021-KF-20), and Tianyou Youth Talent Lift Program of Lanzhou Jiaotong University.

Ushasi Roy

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Prussin, S. Generation and distribution of dislocations by solute diffusion. J Appl Phys; 1961; [DOI: https://dx.doi.org/10.1063/1.1728256]

2. Kim, H; Kim, J; Kim, D; Kim, GH; Kwon, O; Cha, H; Choi, H; Yoo, H; Park, T. Mass diffusion characteristics on performance of polymer electrolyte membrane fuel cells with serpentine channels of different width. Int J Heat Mass Transfer; 2022; 183, 122106. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.122106]

3. Yang, FQ. Diffusion-induced bending of viscoelastic beams. Int J Mech Sci; 2017; 131, pp. 137-145. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2017.06.055]

4. Liu, L; Zheng, LC; Zhang, XX. Fractional anomalous diffusion with Cattaneo–Christov flux effects in a comb-like structure. Appl Math Model; 2016; 40, pp. 6663-6675.3504900 [DOI: https://dx.doi.org/10.1016/j.apm.2016.02.013]

5. Aouadi, M; Passarella, F; Tibullo, V. A bending theory of thermoelastic diffusion plates based on Green–Naghdi theory. Eur J Mech A/Solids; 2017; 65, pp. 123-135.3667151 [DOI: https://dx.doi.org/10.1016/j.euromechsol.2017.03.001]

6. Le, TD; Lasseux, D; Nguyen, XP; Vignoles, G; Mano, N; Kuhn, A. Multi-scale modeling of diffusion and electrochemical reactions in porous micro-electrodes. Chem Eng Sci; 2017; 173, pp. 153-167. [DOI: https://dx.doi.org/10.1016/j.ces.2017.07.039]

7. Velichko, O. Simulation of the diffusion of point defects in structures with local elastic stresses. Appl Math Model; 2011; 35, pp. 1134-1141.2740855 [DOI: https://dx.doi.org/10.1016/j.apm.2010.08.002]

8. Wang, X; Wang, CY. Transient response of bilayered elastic strips with interfacial diffusion and sliding in cylindrical bending. Math Mech Solids; 2018; 23, pp. 748-774.3799708 [DOI: https://dx.doi.org/10.1177/1081286516689383]

9. Wang, X; Wang, CY; Schiavone, P. In-plane deformations of nano-sized circular inhomogeneity with interface slip and diffusion. Int J Eng Sci; 2016; 108, pp. 9-15.3546239 [DOI: https://dx.doi.org/10.1016/j.ijengsci.2016.08.006]

10. Raghavan, R; Wheeler, JM; Harzer, TP; Chawla, V; Djaziri, S; Thomas, K; Philippi, B; Kirchelechner, C; Jaya, BN; Wehrs, J; Michler, J; Dehm, G. Transition from shear to stress-assisted diffusion of copper-chromium nanolayered thin films at elevated temperatures. Acta Mater; 2015; 100, pp. 73-80. [DOI: https://dx.doi.org/10.1016/j.actamat.2015.08.016]

11. Xu, BB; Wei, AR; Ye, JY; Sha, WH; Guo, FL. Transverse creep induced by interfacial diffusion in fiber-reinforced composites: a micromechanics model and computational validation. Int J Mech Sci; 2022; 213, 106877. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2021.106877]

12. Li, Y; Zhang, K; Zheng, BL; Yang, FQ. Effect of local deformation on the coupling between diffusion and stress in lithium-ion battery. Int J Solids Struct; 2016; 87, pp. 81-89. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2016.02.029]

13. Yang, FQ. Interaction between diffusion and chemical stresses. Mater Sci Eng A; 2005; 409, pp. 153-159. [DOI: https://dx.doi.org/10.1016/j.msea.2005.05.117]

14. Sobolev, SL. Diffusion-stress coupling in liquid phase during rapid solidification of binary mixtures. Phys Lett A; 2014; 378, pp. 475-479. [DOI: https://dx.doi.org/10.1016/j.physleta.2013.12.031]

15. Sobolev, SL. Local non-equilibrium diffusion model for solute trapping during rapid solidification. Acta Mater; 2012; 60, pp. 2711-2718. [DOI: https://dx.doi.org/10.1016/j.actamat.2012.01.036]

16. Sobolev, SL. Rapid colloidal solidifications under local nonequilibrium diffusion conditions. Phys Lett A; 2012; 376, pp. 3563-3566. [DOI: https://dx.doi.org/10.1016/j.physleta.2012.10.031]

17. Suo, YH; Shen, SP. Dynamical theoretical model and variational principles for coupled temperature-diffusion-mechanics. Acta Mech; 2012; 223, pp. 29-41.2875877 [DOI: https://dx.doi.org/10.1007/s00707-011-0545-4]

18. Kuang, ZB. Energy and entropy equations in coupled nonequilibrium thermal mechanical diffusive chemical heterogeneous system. Sci Bull; 2015; 60, pp. 952-957. [DOI: https://dx.doi.org/10.1007/s11434-015-0784-1]

19. Li, CL; Guo, HL; Tian, XG; He, TH. Time-domain finite element method to generalized diffusion–elasticity problems with the concentration-dependent elastic constants and the diffusivity. Appl Math Model; 2020; 87, pp. 55-76.4109981 [DOI: https://dx.doi.org/10.1016/j.apm.2020.05.004]

20. Li, CL; Lu, YN; He, TH. An investigation into size-dependent dynamic mechanical-diffusion responses of multi-layered laminated sandwich-like nanocomposites under shock loadings of molar concentration for vibration control based on nonlocal diffusion–elasticity theory. Mech Adv Mater Struct; 2022; 30, pp. 647-660. [DOI: https://dx.doi.org/10.1080/15376494.2021.2021462]

21. Li, CL; Lu, YN; Guo, HL; He, TH; Tian, XG. Non-fick diffusion–elasticity based on a new nonlocal dual-phase-lag diffusion model and its application in structural transient dynamic responses. Acta Mech; 2023; 234, pp. 2745-2761.4599008 [DOI: https://dx.doi.org/10.1007/s00707-023-03519-0]

22. Lu, YN; Li, CL; He, TH. Fractional-order non-Fick mechanical-diffusion coupling model based on new fractional derivatives and structural transient dynamic responses of multilayered composite laminates. Arch Appl Mech; 2024; 94, pp. 239-259. [DOI: https://dx.doi.org/10.1007/s00419-023-02518-w]

23. Alhassan, Y; Alsubhi, M; Abouelregal, AE. A modified Moore–Gibson–Thompson fractional model for mass diffusion and thermal behavior in an infinite elastic medium with a cylindrical cavity. AIMS Math; 2024; 9, pp. 21860-21889.4771843 [DOI: https://dx.doi.org/10.3934/math.20241063]

24. Suo, YH; Shen, SP. Analytical solution for one-dimensional coupled non-Fick diffusion mechanics. Arch Appl Mech; 2013; 83, pp. 397-411. [DOI: https://dx.doi.org/10.1007/s00419-012-0687-4]

25. Hosseini, SA; Abolbashari, MH; Hosseini, SM. Shock-induced molar concentration wave propagation and coupled non-Fick diffusion–elasticity analysis using an analytical method. Acta Mech; 2013; 225, pp. 3591-3599.3279571 [DOI: https://dx.doi.org/10.1007/s00707-014-1161-x]

26. Hosseini, SM; Sladek, J; Sladek, V. Application of meshless local integral equations to two-dimensional analysis of coupled non-Fick diffusion–elasticity. Eng Anal Bound Elem; 2013; 37, pp. 603-615.3036600 [DOI: https://dx.doi.org/10.1016/j.enganabound.2013.01.010]

27. Hosseini, SM. Shock-induced two dimensional coupled non-Fickian diffusion–elasticity analysis using meshless generalized finite difference (GFD) method. Eng Anal Bound Elem; 2015; 61, pp. 232-240.3400010 [DOI: https://dx.doi.org/10.1016/j.enganabound.2015.07.019]

28. Aouadi, M. Hybrid laplace transform-finite element method to a generalized electromagneto-thermoelastic problem. Appl Math Model; 2007; 31, pp. 712-726. [DOI: https://dx.doi.org/10.1016/j.apm.2005.11.030]

29. Tian, XG; Shen, YP; Chen, CQ; He, TH. A direct finite element method study of generalized thermoelastic problems. Int J Solids Struct; 2006; 43, pp. 2050-2063. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2005.06.071]

30. Pankratov, EL. Influence of the spatial, temporal, and concentrational dependence of the diffusion coefficient on dopant dynamics: optimization of annealing time. Phys Rev B; 2005; 72, 075201. [DOI: https://dx.doi.org/10.1103/PhysRevB.72.075201]

31. Li, Y; Mei, R; Xu, Y; Kurths, J; Duan, J; Metzler, R. Particle dynamics and transport enhancement in a confined channel with position-dependent diffusivity. New J Phys; 2020; 22, 053016.4127583 [DOI: https://dx.doi.org/10.1088/1367-2630/ab81b9]

32. Garberoglio, G; Vallauri, R. Space-dependent diffusion of hydrogen in carbon nanotubes. J Mol Liq; 2005; 117, pp. 43-47. [DOI: https://dx.doi.org/10.1016/j.molliq.2004.08.029]

33. Lançon, P; Batrouni, G; Lobry, L; Ostrowsky, N. Brownian walker in a confined geometry leading to a space-dependent diffusion coefficient. Physica A; 2002; 304, pp. 65-76. [DOI: https://dx.doi.org/10.1016/S0378-4371(01)00510-6]

34. Saleh, M; Kovács, E; Barna, IF; Mátyás, L. New analytical results and comparison of 14 numerical schemes for the diffusion equation with space-dependent diffusion coefficient. Mathematics; 2022; 10, 2813. [DOI: https://dx.doi.org/10.3390/math10152813]

35. Li, GS; Gu, WJ; Jia, XZ. Numerical inversions for space-dependent diffusion coefficient in the time fractional diffusion equation. J Invers Ill-Posed Probl; 2012; 20, pp. 339-366.2984493 [DOI: https://dx.doi.org/10.1515/jip-2011-0012]

36. Li, GS; Zhang, DL; Jia, XZ; Yamamoto, M. Simultaneous inversion for the space-dependent diffusion coefficient and the fractional order in the time-fractional diffusion equation. Invers Probl; 2013; 29, 065014.3066390 [DOI: https://dx.doi.org/10.1088/0266-5611/29/6/065014]

37. Chen, J; Cui, BT; Chen, YQ. Backstepping-based boundary control design for a fractional reaction diffusion system with a space-dependent diffusion coefficient. ISA Trans; 2018; 80, pp. 203-211. [DOI: https://dx.doi.org/10.1016/j.isatra.2018.04.013]

38. Smyshlyaev, A; Krstic, M. On control design for PDEs with space-dependent diffusivity or time-dependent reactivity. Automatica; 2005; 41, pp. 1601-1608.2161123 [DOI: https://dx.doi.org/10.1016/j.automatica.2005.04.006]

39. Joseph, DD; Preziosi, L. Heat waves. Rev Mod Phys; 1989; 61, pp. 41-73.977943 [DOI: https://dx.doi.org/10.1103/RevModPhys.61.41]

40. Wang, GX; Prasad, V. Microscale heat and mass transfer and non-equilibrium phase change in rapid solidification. Mater Sci Eng A; 2000; 292, pp. 142-148. [DOI: https://dx.doi.org/10.1016/S0921-5093(00)01003-0]

41. Sobolev, SL. On the transition from diffusion-limited to kinetic-limited regimes of alloy solidification. Acta Mater; 2013; 61, pp. 7881-7888. [DOI: https://dx.doi.org/10.1016/j.actamat.2013.09.027]

42. Sobolev, SL. Nonlocal diffusion models: application to rapid solidification of binary mixtures. Int J Heat Mass Transf; 2014; 71, pp. 295-302. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.12.048]