Research article

Automatic detection method of abnormal vibration of engineering electric drive construction machinery

  • Aiming at the problem that the extraction effect of abnormal vibration characteristics of current engineering electric drive construction machinery is poor, an automatic detection method of abnormal vibration of engineering electric drive construction machinery is proposed. Firstly, the abnormal data of mechanical abnormal vibration are collected and identified, and based on the identification results, the dynamic characteristic model of engineering electric drive construction machinery is constructed. The empirical mode decomposition and Hilbert spectrum are used to decompose the abnormal vibration of machinery, calculate the response amplitude and time lag value generated by the operation of the engineering electric drive construction machinery to simplify the diagnosis steps of the abnormal vibration of the engineering electric drive construction machinery and realize the positioning and detection of the transverse and torsional vibration characteristics. Finally, through experiments, it was confirmed that the automatic detection method of the abnormal vibration of the engineering electric drive construction machinery has high accuracy, which can better ensure the healthy operation of mechanical equipment. This endeavor aims to establish scientific methodologies and standards for fault detection techniques in construction machinery, ultimately forging a versatile solution better suited for detecting and resolving issues across various categories of construction equipment.

    Citation: Jian Yuan, Hao Liu, Yang Zhang. Automatic detection method of abnormal vibration of engineering electric drive construction machinery[J]. Electronic Research Archive, 2023, 31(10): 6327-6346. doi: 10.3934/era.2023320

    Related Papers:

    [1] Heesung Shin, Jiang Zeng . More bijections for Entringer and Arnold families. Electronic Research Archive, 2021, 29(2): 2167-2185. doi: 10.3934/era.2020111
    [2] Bin Han . Some multivariate polynomials for doubled permutations. Electronic Research Archive, 2021, 29(2): 1925-1944. doi: 10.3934/era.2020098
    [3] Shishuo Fu, Zhicong Lin, Yaling Wang . Refined Wilf-equivalences by Comtet statistics. Electronic Research Archive, 2021, 29(5): 2877-2913. doi: 10.3934/era.2021018
    [4] Jiafan Zhang . On the distribution of primitive roots and Lehmer numbers. Electronic Research Archive, 2023, 31(11): 6913-6927. doi: 10.3934/era.2023350
    [5] Dmitry Krachun, Zhi-Wei Sun . On sums of four pentagonal numbers with coefficients. Electronic Research Archive, 2020, 28(1): 559-566. doi: 10.3934/era.2020029
    [6] Massimo Grossi . On the number of critical points of solutions of semilinear elliptic equations. Electronic Research Archive, 2021, 29(6): 4215-4228. doi: 10.3934/era.2021080
    [7] Ji-Cai Liu . Proof of Sun's conjectural supercongruence involving Catalan numbers. Electronic Research Archive, 2020, 28(2): 1023-1030. doi: 10.3934/era.2020054
    [8] Hongjian Li, Kaili Yang, Pingzhi Yuan . The asymptotic behavior of the reciprocal sum of generalized Fibonacci numbers. Electronic Research Archive, 2025, 33(1): 409-432. doi: 10.3934/era.2025020
    [9] Hanpeng Gao, Yunlong Zhou, Yuanfeng Zhang . Sincere wide τ-tilting modules. Electronic Research Archive, 2025, 33(4): 2275-2284. doi: 10.3934/era.2025099
    [10] Taboka Prince Chalebgwa, Sidney A. Morris . Number theoretic subsets of the real line of full or null measure. Electronic Research Archive, 2025, 33(2): 1037-1044. doi: 10.3934/era.2025046
  • Aiming at the problem that the extraction effect of abnormal vibration characteristics of current engineering electric drive construction machinery is poor, an automatic detection method of abnormal vibration of engineering electric drive construction machinery is proposed. Firstly, the abnormal data of mechanical abnormal vibration are collected and identified, and based on the identification results, the dynamic characteristic model of engineering electric drive construction machinery is constructed. The empirical mode decomposition and Hilbert spectrum are used to decompose the abnormal vibration of machinery, calculate the response amplitude and time lag value generated by the operation of the engineering electric drive construction machinery to simplify the diagnosis steps of the abnormal vibration of the engineering electric drive construction machinery and realize the positioning and detection of the transverse and torsional vibration characteristics. Finally, through experiments, it was confirmed that the automatic detection method of the abnormal vibration of the engineering electric drive construction machinery has high accuracy, which can better ensure the healthy operation of mechanical equipment. This endeavor aims to establish scientific methodologies and standards for fault detection techniques in construction machinery, ultimately forging a versatile solution better suited for detecting and resolving issues across various categories of construction equipment.



    Poroelastic materials play an important role in various fields such as geophysics [1,2,3], biomechanics, chemical engineering [4,5], and materials science [6,7]. The physical process in porous media is usually described by a poroelasticity model. In homogeneous isotropic linear elastic porous media, the governing equations for the quasi-static poroelasticity model are given as

    2divσ(u)+αp=fin ΩT:=Ω×(0,T)Rd×(0,T), (1.1)
    (c0p+αdivu)t+divvf=ϕin ΩT, (1.2)

    where

    σ(u):=με(u)+λtr(ε(u))I,ε(u):=12(u+Tu), (1.3)
    vf:=Kμfp. (1.4)

    Here u and p represent the displacement vector of solid and the pressure of solvent, σ(u) is the (effective) stress tensor, vf is the volumetric solvent flux, and Eq (1.4) is so-called Darcy's law. I denotes the d×d identity matrix and ε(u) is the Green strain tensor. f is the body force and ϕ is the source term. The permeability tensor K=K(x) is assumed to be symmetric and uniformly positive definite in the sense that there exist positive constants K1 and K2 such that K1|ζ|2K(x)ζζK2|ζ|2 for a.e. xΩ and ζRd. The parameters μf,α and c0 are the solvent viscosity, Biot-Willis constant and constrained specific storage coefficient, respectively. λ and μ (Note that we use μ instead of 2μ in (1.3), as often seen in the literature, to simplify the notation) are Lamé constants, which can be computed from the Young's modulus E and Poisson ratio ν as

    λ=Eν(1+ν)(12ν),μ=E2(1+ν).

    Much effort has been carried out for the above poroelasticity model. Biot [1] proposed the three-dimensional consolidation theory, providing important theoretical support for the development of the poroelasticity model. Showalter [8] showed the well-posedness of a weak solution of a linear poroelasticity model. Phillips and Wheeler proposed and analyzed a continuous-in-time linear poroelasticity model [9] and a discrete-in-time linear poroelasticity model [10] with the continuous Galerkin mixed finite element method (MFEM) in space. After that, a discontinuous Galerkin MFEM for the model is also developed [11]. The nonconforming FEM was studied for the linear poroelasticity model [12,13], and the MINI and stabilized P1P1 elements were employed to solve the linear poroelasticity model [14]. Ju et al. [15] developed parameter-robust numerical algorithms for the Biot model and applied the algorithms to simulations of brain edema. Yi [16] studied two kinds of locking phenomena (Poisson locking and pressure oscillation) in linear poroelasticity and proposed a new family of MFEM to overcome them. To reveal the underlying deformation and diffusion processes of the original model, Feng et al. reformulated the linear poroelasticity model to a fluid-fluid coupled system and proposed a multiphysics MFEM to overcome the locking phenomena [17]. One can refer to more numerical methods for solving the poroelasticity model, such as the Galerkin least square methods [18,19], two-field finite element method [20,21], four-field finite element method[22,23], finite volume method [24], and Weak Galerkin finite element method [25,26].

    In practical problems, the model parameters have different ranges, such as the Poisson ratio ν=0.10.3 and the permeability K=10141016m2 in geophysics [6]. This requires us to construct numerical methods to be robust for the selection of model parameters. Following the idea of [17], we first reformulate the original problem into a four-field by introducing some auxiliary variables. After that, we discretize the four-field problem by the backward Euler method in time and design the coupled and decoupled algorithms by using MFEM with P2P1P1P1 element pairs in space. The inf-sup condition is proved to be uniformly independent of the physical parameters, and the optimal convergence rates both in space and time are also obtained for the proposed algorithms. It is pointed out that, for the decoupled algorithm [17], there exists the constraint τ=O(h2) on the mesh size h and time step τ while performing the error estimates. But there is no constraint on h and τ for the decoupled algorithm arising from our four-field problem. The reason is that the decoupled algorithm ensures that the variables p,ξ, and η are in the same time layer, thereby avoiding the use of inverse inequality. To improve computational efficiency and achieve high-accuracy in time, we further design a two-step backward differentiation formula (BDF2) Galerkin algorithm, which combines the MFEM with P2P1P1 element pairs for the space variables. The corresponding stability and error estimates of the BDF2 algorithm are performed, and the optimal convergence rate in time is also obtained in Theorem 7. Numerical tests are also provided to verify that our proposed algorithms are robust for selection of ν and K, overcome the Poisson locking in stress field and nonphysical oscillation in the pressure field, and have the optimal convergence rates. In addition, our algorithms are also used to simulate the brain edema. The numerical results on pressure are consistent with the experimental results [27]. The influence of the Poisson rate ν on the displacement of brain tissue and maximum pressure is investigated, which shows that the values of ν have little influence on the maximum pressure, and the maximum tissue deformation becomes smaller as ν gets bigger.

    The remainder of this paper is organized as follows: In Section 2, the original problem is equivalently transformed into a four-field problem, and the four-field problem is also proved to satisfy the inf-sup condition. In Section 3, the coupled and decoupled algorithms with the Euler method are proposed for the four-field problem, and the error estimates are also presented. In Section 4, the BDF2 algorithm for the poroelasticity problem is also designed, and the stability analysis and error estimates are established. In Section 5, numerical examples are provided to verify that the algorithms are robust with respect to the key physical parameters, reach the optimal convergence rates, and overcome the locking phenomenon. In Section 6, the algorithms are applied to simulate the brain edema. Finally, we draw a conclusion to summarize the main results of this paper.

    The standard function space notations are adopted; their precise definitions can be found in [28,29,30]. In particular, (,) and , denote respectively the standard L2(Ω) and L2(Ω) inner products, where ΩRd(d=2,3) is a bounded polygonal domain with the boundary Ω. For any Banach space B, we let B=[B]d, and use B to denote its dual space. We also introduce the function spaces

    L20(Ω):={qL2(Ω);(q,1)=0},X:=H1(Ω).

    Define the space of infinitesimal rigid motions by RM={r:=a+b×x;a,b,xRd}. As in [28,30,31], RM is the kernel of the operator ε, i.e., rRM if and only if ε(r)=0. Hence, we have

    ε(r)=0,divr=0,rRM. (2.1)

    Denote L2(Ω) and H1(Ω) respectively by the subspaces of L2(Ω) and H1(Ω), which are orthogonal to RM, and are given as,

    H1(Ω):={vH1(Ω);(v,r)=0rRM},L2(Ω):={gL2(Ω);g,r=0rRM}.

    It follows from [17,32] that there exists a constant ck>0 such that

    vL2(Ω)=infrRMv+rL2(Ω)rL2(Ω)ckε(v)L2(Ω),vH1(Ω). (2.2)

    To close the above system (1.1)–(1.4), here we set the following boundary and initial conditions:

    2σ(u)nαpn=f1on ΩT:=Ω×(0,T), (2.3)
    1μfKpn=ϕ1on ΩT, (2.4)
    u=u0,p=p0in Ω×{t=0}. (2.5)

    By introducing the variables q=divu,η=c0p+αq,ξ=αpλq, one has

    p=κ1ξ+κ2η,q=κ1ηκ3ξ, (2.6)

    with

    κ1=αα2+λc0,κ2=λα2+λc0,κ3=c0α2+λc0.

    Thus, the problem (1.1)–(1.4) can be equivalently written as

    2μdivε(u)+ξ=fin ΩT, (2.7)
    κ3ξ+divu=κ1ηin ΩT, (2.8)
    κ1ξ+κ2ηp=0in ΩT, (2.9)
    ηt1μfdiv(Kp)=ϕin ΩT. (2.10)

    Remark 1. We point out that the following results with respect to the pure Neumann condition can be extended to the case of the pure Dirichlet condition or mixed Dirichlet-Neumann condition. If the boundary condition (2.3) is replaced by a pure Dirichlet condition or the mixed Dirichlet-Neumann condition, there is no need to introduce the space H1(Ω) for the variable u.

    We here consider an equidistant partition of the time interval [0,T] with a time step τ=T/N and tn=nτ (n=1,2,,N), where NN. The focus of this section is on the proof of the inf-sup condition, which will hold uniformly independently of physical parameters. To do that, we first discretize the time-dependent problem (2.7)–(2.10) with (2.3)–(2.5) by the backward Euler method and BDF2 scheme in time to make it a static problem in each time step. Note that the discussion on inf-sup condition is similar for the backward Euler method and BDF2 scheme. Below, we mainly consider the backward Euler method. By the backward Euler method, the weak formulation of problem (2.7)–(2.10) with (2.3)–(2.5) reads: (v,φ,ω,ψ)H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω), (u,ξ,η,p)H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω) such that

    2μ(ε(u),ε(v))(ξ,divv)=(f,v)+f1,v, (2.11)
    κ3(ξ,φ)+(divu,φ)=κ1(η,φ), (2.12)
    (κ1ξ,ω)+(κ2η,ω)(p,ω)=0, (2.13)
    (η,ψ)+τμf(Kp,ψ)=(˜ϕ,ψ)+˜ϕ1,ψ, (2.14)

    where ˜ϕ(x,tn)=τϕ(x,tn)+η(x,tn1) and ˜ϕ1=τϕ1.

    Define

    A(u,ξ,η,p)=(μdivε00divκ3κ100κ1κ21001τμfdivK)(uξηp).

    Taking Λ=H1(Ω)×L2(Ω)×L2(Ω)×H1(Ω) and ξ0 as the mean-value zero part of ξ, we define

    (u,ξ,η,p)Λ=μ12ε(u)L2(Ω)+κ123ξL2(Ω)+δ120ξ0L2(Ω)+κ122ηL2(Ω)+(τμf)12K12pL2(Ω), (2.15)

    where δ0 is a positive constant, uH1(Ω), ξL2(Ω), ηL2(Ω) and pH1(Ω).

    The following theorem shows that A is invertible and the inverse is a map from Λ to Λ.

    Theorem 1. There exists a constant β>0 for the problem (2.11)(2.14), independent of μ,κ2,κ3 and Kτ/μf, such that the following inf-sup condition holds:

    inf(u,ξ,η,p)Λsup(v,φ,ω,ψ)Λ(A(u,ξ,η,p),(v,φ,ω,ψ))(u,ξ,η,p)Λ(v,ϕ,ω,ψ)Λβ. (2.16)

    Proof. To prove (2.16), we need to prove the following inequalities: (0,0,0,0)(u,ξ,η,p)Λ, (v,φ,ω,ψ)Λ such that

    (v,φ,ω,ψ)Λβ1(u,ξ,η,p)Λ,(A(u,ξ,η,p),(v,φ,ω,ψ))β2(u,ξ,η,p)2Λ,

    where β1 and β2 are positive constants and independent of μ,κ1,κ2,κ3 and Kτ/μf.

    Suppose that (0,0,0,0)(u,ξ,η,p)Λ is given. It follows from the theory of the Stokes equation that there exists a constant β0>0, depending only on the domain Ω, and wH1(Ω) such that

    (divw,ξ)=ξ02L2(Ω),(ε(w),ε(w))β20ξ02L2(Ω). (2.17)

    Taking v=uδ0w, ϕ=ξ, ω=η, ψ=p in (A(u,ξ,η,p),(v,φ,ω,ψ)), we have

    (A(u,ξ,η,p),(v,φ,ω,ψ))=μ(ε(u),ε(u))δ0μ(ε(u),ε(w))(ξ,divu)+δ0(ξ,divw)+κ3(ξ,ξ)+(divu,ξ)κ1(η,ξ)+(κ1ξ+κ2η,η)(p,η)+(η,p)+τμfKp2L2(Ω)=με(u)2L2(Ω)δ0μ(ε(u),ε(w))+δ0(ξ,divw)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω). (2.18)

    Using Young's inequality and (2.17), we have

    (ε(u),ε(w))12θ0(ε(u),ε(u))+θ02(ε(w),ε(w))12θ0(ε(u),ε(u))+θ0β202(ξ0,ξ0),θ0>0. (2.19)

    Inserting (2.19) into (2.18), using (2.17) and choosing θ0=δ0=β20μ1, we obtain

    (A(u,ξ,η,p),(v,φ,ω,ψ))μ2ε(u)2L2(Ω)+δ02ξ02L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω). (2.20)

    Noting that

    (u,ξ,η,p)2Λ=(μ12ε(u)L2(Ω)+κ123ξL2(Ω)+κ122ηL2(Ω)+(τμf)12KpL2(Ω))24(με(u)2L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω))4(με(u)2L2(Ω)+δ0ξ02L2(Ω)+κ3ξ2L2(Ω)+κ2η2L2(Ω)+τμfKp2L2(Ω)), (2.21)

    and using (2.21) in (2.20), we have

    (A(u,ξ,η,p),(v,φ,ω,ψ))18(u,ξ,η,p)2Λ. (2.22)

    It is easy to check that

    (v,φ,ω,ψ)Λ(1+μ12δ120β0)(u,ξ,η,p)Λ. (2.23)

    Combining (2.22) and (2.23), the (2.16) holds and β=1/(8(1+μ12δ120β0)). The proof is complete.

    Theorem 1 shows that the bilinear form (A(),) is coercive. Moreover, it is easy to check that (A(),) is bounded. According to the Lax-Milgram theorem, there exists a unique solution (u,ξ,η,p)Λ for the linear system (2.11)–(2.14). Thus, there holds the following theorem:

    Theorem 2. There exists a unique weak solution to problem (2.7)(2.10) with (2.3)(2.5). Likewise, there exists a unique weak solution to problem (1.1)(1.4) and (2.3)(2.5).

    Let Th be a quasi-uniform triangulation or rectangular partition of Ω with maximum mesh size h, and ˉΩ=KThˉK. In this part, we use the backward Euler method in time discretization. For a smooth function v defined on [0,T], we denote dtvn=(vnvn1)/τ.

    A number of stable mixed finite element spaces (Xh,Mh) have been presented in the literature [33]. A typical example is the following so-called Taylor-Hood element (cf. [33,34]):

    Xh={vhC0(¯Ω);vh|KP2(K)KTh},Mh={φhC0(¯Ω);φh|KP1(K)KTh}.

    Finite element approximation space Wh for p and η variables can be chosen independently; any piecewise polynomial space is acceptable provided that WhMh. A convenient choice is Wh=Mh.

    Define

    Vh:={vhXh;(vh,r)=0,rRM}. (3.1)

    If the problem (2.11)–(2.14) is discretized by the finite element spaces Vh,Mh, and Wh, then the discrete counterpart is to find (uh,ξh,ηh,ph)Vh×Mh×Wh×Wh such that

    2μ(ε(uh),ε(vh))(ξh,divvh)=(f,vh)+f1,vh,vhVh, (3.2)
    κ3(ξh,φh)+(divuh,φh)=κ1(ηh,φh),φhMh, (3.3)
    κ1(ξh,ωh)+κ2(ηh,ωh)(p,ωh)=0,ωhWh, (3.4)
    (ηh,ψh)+τμf(Kph,ψh)=(˜ϕ,ψh)+˜ϕ1,ψh,ψhWh. (3.5)

    Set Λh=Vh×Mh×Wh×Wh, and define the operator Ah:ΛhΛh arisen from (3.2)–(3.5). Then we have the following theorem:

    Theorem 3. For the problem (3.2)(3.5), there exists a constant β>0, independent of μ,κ2,κ3, and Kτ/μf, such that the following inf-sup condition holds:

    inf(uh,ξh,ηh,ph)Λhsup(vh,φh,ωh,ψh)Λh(Ah(uh,ξh,ηh,ph),(vh,φh,ωh,ψh))(uh,ξh,ηh,ph)Λh(vh,φh,ωh,ψh)Λhβ. (3.6)

    We omit the proof of Theorem 3, since it is completely analogous to the proof of Theorem 1. Moreover, Theorem 3 implies that Ah is invertible and the inverse is a map from Λh to Λh with operator norm independent of the parameters.

    Algorithm 1 A coupled algorithm
    Input: Evaluate u0hVh,p0hWh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      for n=0,1,2,, do
        solve for (un+1h,ξn+1h,ηn+1h,pn+1h)Vh×Mh×Wh×Wh such that
                μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηn+1h,φh),φhMh,κ1(ξn+1h,ωh)+κ2(ηn+1h,ωh)(pn+1h,ωh)=0,ωhWh,(dtηn+1h,ψh)+1μf(Kpn+1h,ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh,
      end for

    Algorithm 2 A decoupled algorithm
    Input: Evaluate u0hVh,p0hWh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      for n=0,1,2,, do
        (i) solve for (un+1h,ξn+1h)Vh×Mh such that
                μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηnh,φh),φhMh.
        (ii) Using (un+1h,ξn+1h) obtained in (i), solve for (ηn+1h,pn+1h)Wh×Wh by
                κ1(ξn+1h,ωh)+κ2(ηn+1h,ωh)(pn+1h,ωh)=0,ωhWh,(dtηn+1h,ψh)+1μf(Kpn+1h,ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh,
      end for

    Below, we assume f,f1, ϕ, and ϕ1 all are independent of t. Note that all the results can be easily extended to the case of time-dependent source functions.

    Now, we consider the fully discrete MFEM for the problem (2.7)–(2.10) with (2.3)–(2.5) based on the backward Euler method, which involves the coupled (i.e., Algorithm 1) and decoupled (i.e., Algorithm 2) algorithms. The coupled algorithm solves four unknowns simultaneously based on problem (2.7)–(2.10). The decoupled algorithm only solves two subproblems. Noting that problem (2.7) and (2.8) can be regarded as a generalized Stokes problem for a given η and problem (2.9) and (2.10) is a coupled problem associated with the diffusion problem, we can solve a generalized Stokes problem for u and ξ by using the solution of η at the previous time-step. After that, we solve the diffusion problem for η and p.

    To derive the error estimates of the MFEM, we define L2(Ω)-projection Qh:L2(Ω)Xkh by

    (Qhφ,ψh)=(φ,ψh),ψhXkh and φL2(Ω), (3.7)

    where Xkh:={ψhC0;ψh|KPk(K)KTh}, k is the degree of piecewise polynomial on K.

    Next, for any φH1(Ω), we define the elliptic projection Sh:H1(Ω)Xkh as

    (KShφ,φh)=(Kφ,φh),φhXkh, (3.8)
    (Shφ,1)=(φ,1), (3.9)

    and for any vH1(Ω), we define its elliptic projection Rh:H1(Ω)Vkh as

    (ε(Rhv),ε(wh))=(ε(v),ε(wh))whVkh, (3.10)

    where Vkh:={vhC0;vh|KPk(K),(vh,r)=0rRM}.

    It follows from [28] that for any 0sk, Qh,Sh and Rh satisfy

    QhφφL2(Ω)+h(Qhφφ)L2(Ω)Chs+1φHs+1(Ω),φHs+1(Ω), (3.11)
    ShφφL2(Ω)+h(Shφφ)L2(Ω)Chs+1φHs+1(Ω),φHs+1(Ω), (3.12)
    RhvvL2(Ω)+h(Rhvv)L2(Ω)Chs+1vHs+1(Ω),vHs+1(Ω). (3.13)

    Furthermore, we introduce the following error notations:

    Enu=u(tn)unh,Enξ=ξ(tn)ξnh,Enη=η(tn)ηnh,Enp=p(tn)pnh.

    Also, we denote

    Enu=u(tn)Rh(u(tn))+Rh(u(tn))unh:=Ynu+Znu,Eni=i(tn)Sh(i(tn))+Sh(i(tn))inh:=Yni+Zni,i=ξ,η,p,Eni=i(tn)Qh(i(tn))+Qh(i(tn))inh:=Fni+Gni,i=ξ,η,p.

    Then we present the error estimates for the coupled and decoupled algorithms by the following theorems.

    Theorem 4. Let {(unh,ξnh,ηnh,pnh)}n0 be solutions of the coupled algorithm (i.e., Algorithm 1). Then it holds

    max0nN[μ(u(tn)unh)L2(Ω)+κ2η(tn)ηnhL2(Ω)+κ3ξ(tn)ξnhL2(Ω)]+(τNn=01μfK12(p(tn)pnh)2L2(Ω))12C1τ+C2h2. (3.14)

    Theorem 5. Let {(unh,ξnh,ηnh,pnh)}n0 be solutions of the decoupled algorithm (i.e., Algorithm 2). Then for the case c0>0 and λ<, it holds

    max0nN[μ(u(tn)unh)L2(Ω)+κ2η(tn)ηnhL2(Ω)+κ3ξ(tn)ξnhL2(Ω)]+(τNn=01μfK12(p(tn)pnh)2L2(Ω))12ˉC1τ+ˉC2h2. (3.15)

    The proof of the aforementioned theorems follows a similar approach to that of [17, Theorem 3.7]; thus, further details are omitted here. Our primary focus is on the constraint condition τ=O(h2) in the error estimates for the algorithm with θ=0 [17]. This constraint arises because the choice of the test function ψh results in the variables p,ξ, and η being in different time layers, which in turn leads to the application of an inverse inequality in the estimation of the right-hand side term. To address this issue, a straightforward solution is to select an appropriate test function ψh that ensures the variables p,ξ, and η reside in the same time layer. To do that, we first examine the following error equations:

    2κ1(Gn+1ξ,ωh)+κ2(Gn+1η,ωh)(Gn+1p,ωh)=0,ωhWh, (3.16)
    (dtGn+1η,ψh)+1μf(KZn+1p,ψh)=(Rn+11,η,ψh),ψhWh, (3.17)

    which are derived from the fourth equation of Algorithm 2 and the continuous problem (2.10) the fourth equation of Algorithm 2 and the continuous problem (2.10). Based on Algorithm 2, we can define η1h by

    κ1(η1h,φh)=κ3(ξ0h,φh)+(divu0h,φh).

    Setting ωh=dtGnη,ψh=Znp in (3.16) and (3.17) separately after lowering the super index n+1 to n, applying the summation operator τln=0 to both sides of (3.17) and utilizing the fact that G1η=0, there holds

    κ22Glη2L2(Ω)+τln=0[κ1(dtGnη,Gn+1ξ)+κ2τ2dtGnη2L2(Ω)+1μfK12Znp2L2(Ω)]=τln=0[(dtGnη,YnpFnp)+(dtη(tn+1)ηt(tn+1),Znp)+κ1τ(dtGn+1ξ,dtGnη)]. (3.18)

    Note that the choice of ψh=Znp is critical. Unlike the setting in [17], where ψh=ˆZnp:=κ1Zn+1ξ+κ2Znη is employed, selecting ψh=Znp eliminates the need to apply the inverse inequality to bound a gradient term (as seen in [17, estimate (3.52)]). Instead, we only need to apply the Cauchy-Schwarz and Young inequalities to the last term on the right-hand side of (3.18) to derive

    κ1τ2ln=0(dtGn+1ξ,dtGnη)κ1τ2ln=0[12ϵdtGn+1ξ2L2(Ω)+ϵ2dtGnη)2L2(Ω)], (3.19)

    where ϵ(2κ2κ3,+) is a positive constant from Young inequality.

    Additionally, we point out that:

    (i) The aim of ϵ(2κ2κ3,+) is to ensure the signs of terms κ2κ1ϵ2τ2ln=0dtGnη2L2(Ω) and (κ34κ12ϵ)τ2ln=0dtGnη2L2(Ω) to be positive. In addition, for the case of λ+, there does not exist the term κ1τ2ln=0(dtGn+1ξ,dtGnη) in error estimates. For the case of c0=0, it has a similar discussion for the term κ1τ2ln=0(dtGn+1ξ,dtGnη).

    (ii) For the decoupled algorithm, it does not exist the constraint on the time step τ and mesh size h.

    We now consider the BDF2 scheme in time discretization for the problem (2.7)–(2.10) with (2.3)–(2.5). For a smooth function v defined on [0,T], we denote

    vn+12:=(vn+1+vn)/2andDtvn=3vn4vn1+vn22τ. (4.1)

    Based on (4.1), we introduce the BDF2 fully discrete scheme in Algorithm 3.

    Algorithm 3 BDF2 algorithm
    Input: Evaluate u0hVh,p0hMh, ξ0hMh by ξ0h=αp0hλdivu0h and η0hWh by η0h=c0p0h+αdivu0h.
      Step 1: Solve for (u1h,ξ1h,η1h)Vh×Mh×Wh such that
            2μ(ε(u12h),ε(vh))(ξ12h,divvh)=(f,vh)+f1,vh,vhVh,(4.2)
            κ3(ξ12h,φh)+(divu12h,φh)=κ1(η12h,φh),φhMh,(4.3)
            (dtη1h,ψh)+1μf(K(κ1ξ12h+κ2η12h),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh.(4.4)
      Step 2:
      for n=1,2,, do
        (i) Solve for (un+1h,ξn+1h,ηn+1h)Vh×Mh×Wh such that
            2μ(ε(un+1h),ε(vh))(ξn+1h,divvh)=(f,vh)+f1,vh,vhVh,(4.5)
            κ3(ξn+1h,φh)+(divun+1h,φh)=κ1(ηn+1h,φh),φhMh,(4.6)
            (Dtηn+1h,ψh)+1μf(K(κ1ξn+1h+κ2ηn+1h),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh.(4.7)
        (ii) Update pn+1h and qn+1h by
            pn+1h=κ1ξn+1h+κ2ηn+1h,qn+1h=κ1ηn+1hκ3ξn+1h,(4.8)
      end for

    The next lemma establishes a discrete energy inequality for the BDF2 algorithm.

    Lemma 4.1. Let {(unh,ξnh,ηnh)}n0 be defined by the BDF2 algorithm. Then it holds

    μ8ε(ul+1h)2L2(Ω)+κ34ξl+1h2L2(Ω)+κ24ηl+1h2L2(Ω)+τ2μfln=1K12(κ1ξn+1h+κ2ηn+1h)2L2(Ω)C3(f2L2(Ω)+f12L2(Ω)+ε(u0h)2L2(Ω)+ξ0h2L2(Ω)+η0h2L2(Ω)+ϕ2L2(Ω)+ϕ12L2(Ω)),l1, (4.9)

    where C3=C3(ck,ca,cb,cp,μ,κ3,κ2,T,K,μf).

    Proof. Setting vh=dtu1h in (4.2), φh=dtξ1h in (4.3) and ψh=p12h in (4.4), we get

    2μ(ε(u12h),ε(dtu1h))(ξ12h,divdtu1h)=(f,dtu1h)+f1,dtu1h, (4.10)
    κ3(ξ12h,dtξ1h)+(divu12h,dtξ1h)=κ1(η12h,dtξ1h), (4.11)
    (dtη1h,p12h)+1μf(K(κ1ξ12h+κ2η12h),p12h)=(ϕ,p12h)+ϕ1,p12h. (4.12)

    Adding (4.10)–(4.12), there holds

    μ2τ(ε(u1h)2L2(Ω)ε(u0h)2L2(Ω))+κ32τ(ξ1h2L2(Ω)ξ0h2L2(Ω))+κ22τ(η1h2L2(Ω)η0h2L2(Ω))+1μfK12(κ1ξ12h+κ2η12h)2L2(Ω)=(f,dtu1h)+f1,dtu1h+(ϕ,p12h)+ϕ1,p12h+1τ((ξ0h,divu1h)(ξ1h,divu0h))+κ1τ((η0h,ξ1h)(η1h,ξ0h)). (4.13)

    Using the Cauchy-Schwarz inequality, triangle inequality, (2.2), (4.32), trace inequality, Poincaré inequality, and Young inequality, we can bound the right-side terms of (4.13) as follows:

    (f,dtu1h)1τ(5c2kμf2L2(Ω)+μ16ε(u1h)2L2(Ω)+μ4ε(u0h)2L2(Ω)), (4.14)
    f1,dtu1h1τ(5c2ac2kμf12L2(Ω)+μ16ε(u1h)2L2(Ω)+μ4ε(u0h)2L2(Ω)), (4.15)
    (ϕ,p12h)c2pμfK12ϕ2L2(Ω)+14μfK12p12h2L2(Ω), (4.16)
    ϕ1,p12hc2ac2pμfK12ϕ12L2(Ω)+14μfK12p12h2L2(Ω), (4.17)
    (ξ0h,divu1h)(ξ1h,divu0h)2c2bμξ0h2L2(Ω)+μ8ε(u1h)2L2(Ω)+κ38ξ1h2L2(Ω) (4.18)
    +2c2bκ3ε(u0h)2L2(Ω),κ1(η0h,ξ1h)κ1(η1h,ξ0h)2κ21κ3η0h2L2(Ω)+κ38ξ1h2L2(Ω)+κ24η1h2L2(Ω)+κ21κ2ξ0h2L2(Ω), (4.19)

    where ca is a positive constant related to Trace inequality.

    Substituting (4.14)–(4.19) into (4.13), we have

    μ4ε(u1h)2L2(Ω)+κ34ξ1h2L2(Ω)+κ24η1h2L2(Ω)+τ2μfK12(κ1ξ12h+κ2η12h)2L2(Ω)5c2kμf2L2(Ω)+5c2ac2kμf12L2(Ω)+c2pμfτK12ϕ2L2(Ω)+c2ac2pμfτK12ϕ12L2(Ω)+(μ+2c2bκ3)ε(u0h)2L2(Ω)+(κ32+2c2bμ+κ21κ2)ξ0h2L2(Ω)+(κ22+2κ21κ3)η0h2L2(Ω). (4.20)

    Next, setting vh=Dtun+1h,φh=ξn+1h and ψh=pn+1h in (4.5)–(4.7), we obtain

    2μ(ε(un+1h),ε(Dtun+1h))(ξn+1h,divDtun+1h)=(f,Dtun+1h)+f1,Dtun+1h, (4.21)
    κ3(Dtξn+1h,ξn+1h)+(divDtun+1h,ξn+1h)=κ1(Dtηn+1h,ξn+1h), (4.22)
    (Dtηn+1h,pn+1h)+1μf(K(κ1ξn+1h+κ2ηn+1h),pn+1h)=(ϕ,pn+1h)+ϕ1,pn+1h. (4.23)

    Adding (4.21)–(4.23) and using the fact of

    (Dtvn+1h,vn+1h)=14τ(vn+1h2L2(Ω)vnh2L2(Ω)+2vn+1hvnh2L2(Ω)2vnhvn1h2L2(Ω)+vn+1h2vnh+vn1h2L2(Ω)), (4.24)

    we have

    μ4τ(ε(un+1h)2L2(Ω)ε(unh)2L2(Ω)+2ε(un+1h)ε(unh)2L2(Ω)2ε(unh)ε(un1h)2L2(Ω)+ε(un+1h)2ε(unh)+ε(un1h)2L2(Ω))+κ34τ(ξn+1h2L2(Ω)ξnh2L2(Ω)+2ξn+1hξnh2L2(Ω)2ξnhξn1h2L2(Ω)+ξn+1h2ξnh+ξn1h2L2(Ω))+κ24τ(ηn+1h2L2(Ω)ηnh2L2(Ω)+2ηn+1hηnh2L2(Ω)2ηnhηn1h2L2(Ω)+ηn+1h2ηnh+ηn1h2L2(Ω))+1μfK12(κ1ξn+1h+κ2ηn+1h)2L2(Ω)=(f,Dtun+1h)+f1,Dtun+1h+(ϕ,pn+1h)+ϕ1,pn+1h. (4.25)

    Applying the summation operator τln=1 to both sides of (4.25), there holds

    μ4(ε(ul+1h)2L2(Ω)+2ε(ul+1h)ε(ulh)2L2(Ω))+κ34(ξl+1h2L2(Ω)+2ξl+1hξlh2L2(Ω))+κ24(ηl+1h2L2(Ω)+2ηl+1hηlh2L2(Ω))+ln=1(μ2ε(un+1h)2ε(unh)+ε(un1h)2L2(Ω)+κ34ξn+1h2ξnh+ξn1h2L2(Ω)+κ24ηn+1h2ηnh+ηn1h2L2(Ω)+τμfK12(κ1ξn+1h+κ2ηn+1h)2L2(Ω))=(f,3ul+1hulh3u1h+u0h2)+f1,3ul+1hulh3u1h+u0h2+τln=1[(ϕ,pn+1h)+ϕ1,pn+1h]+μ4(ε(u1h)2L2(Ω)+2ε(u1h)ε(u0h)2L2(Ω))+κ34(ξ1h2L2(Ω)+2ξ1hξ0h2L2(Ω))+κ24(η1h2L2(Ω)+2η1hη0h2L2(Ω)). (4.26)

    Similar to (4.13), we can bound the first four terms of the right side of (4.26) as follows:

    (f,3ul+1hulh3u1h+u0h2)4c2kμf2L2(Ω)+μ162ε(ul+1h)ε(ulh)2L2(Ω) (4.27)
    +μ16ε(ul+1h)2L2(Ω)+3μ8ε(u1h)2L2(Ω)+μ8ε(u0h)2L2(Ω),f1,3ul+1hulh3u1h+u0h24c2ac2kμf12L2(Ω)+μ162ε(ul+1h)ε(ulh)2L2(Ω) (4.28)
    +μ16ε(ul+1h)2L2(Ω)+3μ8ε(u1h)2L2(Ω)+μ8ε(u0h)2L2(Ω),τln=1(ϕ,pn+1h)τln=1(c2pμfK12ϕ2L2(Ω)+14μfK12pn+1h2L2(Ω)), (4.29)
    τln=1ϕ1,p12hτln=1(c2ac2pμfK12ϕ12L2(Ω)+14μfK12p12h2L2(Ω)). (4.30)

    Substituting (4.27)–(4.30) into (4.26) and using (4.20) for the resulting inequality, we can deduce that (4.9) holds. The proof is complete.

    In this subsection, we consider the optimal convergence analysis of the BDF2 algorithm. Before that, we have the following estimate for the BDF2 algorithm.

    Theorem 6. Let {(unh,ξnh,ηnh)}n0 be defined by the BDF2 algorithm. Then there exists a positive constant C4 independent of mesh size h and time step τ such that

    μ8ε(Zl+1u)2L2(Ω)+κ34Gl+1ξ2L2(Ω)+κ28Gl+1η2L2(Ω)+τ2μfln=1K12(κ1Zn+1ξ+κ2Zn+1η)2L2(Ω)C4(h4+τ4). (4.31)

    Proof. First, it is easy to check that there exists a positive constant cb such that

    divwL2(Ω)cbε(w)L2(Ω),wH1(Ω). (4.32)

    For n1, at t=tn+1, we have for (2.7)–(2.10) that

    2μ(ε(u(tn+1)),ε(vh))(ξ(tn+1),divvh)=(f,vh)+f1,vh,vhVh, (4.33)
    κ3(ξ(tn+1),φh)+(divu(tn+1),φh)=κ1(η(tn+1),φh),φhMh, (4.34)
    (ηt(tn+1),ψh)+1μf(K(κ1ξ(tn+1)+κ2η(tn+1)),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh. (4.35)

    Moreover, we have from (2.7)–(2.10) at t=t0 and t=t1 that

    2μ(ε(u(t12)),ε(vh))(ξ(t12),divvh)=(f,vh)+f1,vh,vhVh, (4.36)
    κ3(ξ(t12),φh)+(divu(t12),φh)=κ1(η(t12),φh),φhMh, (4.37)
    (ηt(t12),ψh)+1μf(K(κ1ξ(t12)+κ2η(t12)),ψh)=(ϕ,ψh)+ϕ1,ψh,ψhWh, (4.38)

    where v(t12)=(v(t0)+v(t1))/2.

    Subtracting (4.2)–(4.4) from (4.36)–(4.38), there holds

    2μ(ε(E12u),ε(vh))(E12ξ,divvh)=0,vhVh, (4.39)
    κ3(E12ξ,φh)+(divE12u,φh)=κ1(E12η,φh),φhMh, (4.40)
    (dtE1η,ψh)+1μf(K(κ1E12ξ+κ2E12η),ψh)=(R1,η,ψh)+(R2,η,ψh),ψhWh, (4.41)

    where

    R1,η:=12τt1t0Γ2(t)ηtttdt,R2,η:=12t1t0Γ(t)ηtttdt,Γ(t):={t1twhentt0+t12,tt0whent<t0+t12.

    Using the definitions of the projection operators Qh,Sh,Rh for (4.39)–(4.41), we have

    2μ(ε(Z12u),ε(vh))(G12ξ,divvh)=(F12ξ,divvh),vhVh, (4.42)
    κ3(G12ξ,φh)+(divZ12u,φh)=κ1(G12η,φh)(divY12u,φh),φhMh, (4.43)
    (dtG1η,ψh)+1μf(K(κ1Z12ξ+κ2Z12η),ψh)=(R1,η,ψh)+(R2,η,ψh),ψhWh. (4.44)

    Setting vh=dtZ1u in (4.42), φh=dtG1ξ in (4.43) and ψh=Z12p in (4.44), we get

    2μ2τ(ε(Z1u)2L2(Ω)ε(Z0u)2L2(Ω))(G12ξ,divdtZ1u)=(F12ξ,divdtZ1u), (4.45)
    κ32τ(G1ξ2L2(Ω)G0ξ2L2(Ω))+(divZ12u,dtG1ξ)=κ1(G12η,dtG1ξ)(divY12u,dtG1ξ), (4.46)
    κ22τ(G1η2L2(Ω)G0η2L2(Ω))+κ1(dtG12η,G12ξ)+1μfK12(κ1Z12ξ+κ2Z12η)2L2(Ω)=(R1,η,Z12p)+(R2,η,Z12p)+(dtG1η,Y12pF12p). (4.47)

    Adding (4.45)–(4.47), there holds

    μ2τε(Z1u)2L2(Ω)+κ32τG1ξ2L2(Ω)+κ22τG1η2L2(Ω)+1μfK12(κ1Z12ξ+κ2Z12η)2L2(Ω)=μ2τε(Z0u)2L2(Ω)+κ32τG0ξ2L2(Ω)+κ22τG0η2L2(Ω)+(F12ξ,divdtZ1u)(divY12u,dtG1ξ)+(R1,η,Z12p)+(R2,η,Z12p)+(dtG1η,Y12pF12p)+1τ((G0ξ,divZ1u)(G1ξ,divZ0u))+κ1τ((G0η,G1ξ)(G1η,G0ξ)). (4.48)

    Then, we bound the fourth to tenth terms on the right side of (4.48). Using the Cauchy-Schwarz inequality, triangle equality, and (4.32), Poincaré and Young inequalities, we have

    (F12ξ,divdtZ1u)12τ(2c2bμF1ξ2L2(Ω)+2c2bμF0ξ2L2(Ω)+μ2ε(Z1u)2L2(Ω) (4.49)
    +μ2ε(Z0u)2L2(Ω)),(divY12u,dtG1ξ)12τ(4c2bκ3ε(Y1u)2L2(Ω)+4c2bκ3ε(Y0u)2L2(Ω)+κ34G1ξ2L2(Ω) (4.50)
    +κ34G0ξ2L2(Ω)),(R1,η,Z12p)c2pμfτ3640K1ηttt2L2(t0,t1;L2(Ω))+14μfK12Z12p2L2(Ω), (4.51)
    (R2,η,Z12p)c2pμfτ396K1ηttt2L2(t0,t1;L2(Ω))+14μfK12Z12p2L2(Ω), (4.52)
    (dtG1η,Y12pF12p)12τ(κ22G1η2L2(Ω)+κ22G0η2L2(Ω)+4κ2Y1p2L2(Ω) (4.53)
    +4κ2Y0p2L2(Ω)+4κ2F1p2L2(Ω)+4κ2F0p2L2(Ω)),(G0ξ,divZ1u)(G1ξ,divZ0u)2c2bμG0ξ2L2(Ω)+μ8ε(Z1u)2L2(Ω)+κ38G1ξ2L2(Ω) (4.54)
    +2c2bκ3ε(Z0u)2L2(Ω),κ1(G0η,G1ξ)κ1(G1η,G0ξ)2κ21κ3G0η2L2(Ω)+κ38G1ξ2L2(Ω)+κ24G1η2L2(Ω)+κ21κ2G0ξ2L2(Ω). (4.55)

    Substituting (4.49)–(4.55) into (4.48) and using (3.11)–(3.13), the facts of Z0u=0,G0ξ=0, and G0η=0, we find there exists a positive constant ˜C=˜C(uL(t0,t1;H2(Ω)),ξL(t0,t1;H2(Ω)),ηL(t0,t1;H2(Ω)),ηtttL2(t0,t1;L2(Ω)),cb,μ,κ3,cp,μf,K1,κ2,κ1) such that

    με(Z1u)2L2(Ω)+κ3G1ξ2L2(Ω)+κ2G1η2L2(Ω)+τμfK(κ1Z12ξ+κ2Z12η)2L2(Ω)˜C(h4+τ4). (4.56)

    Next, we show the error estimate at t=tn+1. Subtracting (4.5)–(4.7) from (4.33)–(4.35), it holds

    2μ(ε(Zn+1u),ε(vh))(Gn+1ξ,divvh)=(Fn+1ξ,divvh),vhVh, (4.57)
    (4.58)
    (4.59)

    where

    Setting in (4.57), in (4.58) after using the operator to both sides of (4.58), in (4.59), using (4.24) and adding the result equation, we get

    (4.60)

    Applying the summation operator to both sides of (4.60), there holds

    (4.61)

    Now, we begin to estimate the first four terms on the right side of (4.61). For the first term, it is easy to check that

    (4.62)

    Applying the Cauchy-Schwarz inequality, (4.32), triangle and Young inequalities to (4.62), we have

    (4.63)

    For the second term, using the Canchy-Schwarz and Young inequalities, we obtain

    (4.64)

    Applying the Canchy-Schwarz inequality, Poincaré inequality, Young inequality, and the fact of

    to the third term, there holds

    (4.65)

    Similar to the first term, we can bound the fourth term as follows:

    (4.66)

    Substituting (4.63)–(4.66) into (4.61), using the facts of , (4.56), (3.11)–(3.13), there exists a positive constant such that

    (4.67)

    Applying the Gronwall inequality to (4.67), one has that (4.31) holds. The proof is complete.

    Furthermore, we state the main result of this part by a theorem.

    Theorem 7. There exists a positive constant independent of and such that the solution of the BDF2 algorithm satisfies the following error estimates:

    (4.68)
    (4.69)

    Proof. Using (3.11)–(3.13), the relation , Theorem 6, (4.56), and the triangle inequality on

    one can see that (4.68) and (4.69) hold. The proof is complete.

    Here the element pairs are employed for the coupled algorithm and decoupled algorithms, and element pairs are employed for the BDF2 algorithm. To verify the robustness of our algorithms with respect to parameters, we compute the spatial errors and convergence rates of the variables and for different values of Poisson rates and permeability tensor separately. In the meantime, we show the effectiveness of our algorithms in overcoming the Poisson locking as . After that, to illustrate that there is no nonphysical oscillation in the pressure field, we use the decoupled algorithm to test a Barry-Mercer problem. Finally, we use the following examples to compute the convergence rates for the above proposed algorithms in subsections 5.1, 5.2, and 5.4.

    Example 1. Let , , , , and . The boundary and initial conditions are

    The exact solutions of (1.1)–(1.4) and (2.3)–(2.5) are and

    We first test the robustness of the coupled and decoupled algorithms on the Poisson ratio by fixing , , , e-7, and for and , respectively.

    Tables 1 and 2 show the spatial errors and convergence rates of the coupled algorithm for and , respectively. Numerical results show that the convergence rates of and the convergence rates of are all approaching . The convergence rates of are around . The convergence rates of are around . It can be concluded that the convergence rates of and are optimal, which is consistent with Theorem 4.

    Table 1.  Spatial errors and convergence rates of and of the coupled algorithm .
    Rate Rate Rate Rate
    4.7630e-2 1.6047e-1 8.5755e-2 4.2637e-1
    5.6649e-3 3.07 4.4075e-2 1.86 1.7262e-2 2.31 2.0452e-1 1.06
    6.5936e-4 3.10 1.1349e-2 1.96 3.8013e-3 2.18 9.9924e-2 1.03
    7.9656e-5 3.05 2.8614e-3 1.99 9.8977e-4 1.94 4.9462e-2 1.01

     | Show Table
    DownLoad: CSV
    Table 2.  Spatial errors and convergence rates of and of the coupled algorithm for .
    Rate Rate Rate Rate
    4.7651e-2 1.6055e-1 8.5700e-2 4.2574e-1
    5.6669e-3 3.07 4.4078e-2 1.86 1.7255e-2 2.31 2.0459e-1 1.06
    6.5948e-4 3.10 1.1349e-2 1.96 3.7745e-3 2.19 9.9988e-2 1.03
    7.9662e-5 3.05 2.8614e-3 1.99 8.7440e-4 2.11 4.9501e-2 1.01

     | Show Table
    DownLoad: CSV

    For the decoupled algorithm, we also take and separately to observe the convergence rates of and . Tables 3 and 4 show that the convergence rates of all variables are optimal.

    Table 3.  Spatial errors and convergence rates of and of the decoupled algorithm for .
    Rate Rate Rate Rate
    4.7630e-2 1.6047e-1 8.5755e-2 4.2637e-1
    5.6649e-3 3.07 4.4075e-2 1.86 1.7262e-2 2.31 2.0452e-1 1.06
    6.5935e-4 3.10 1.1349e-2 1.96 3.8013e-3 2.18 9.9924e-2 1.03
    7.9656e-5 3.05 2.8614e-3 1.99 9.8977e-4 1.94 4.9462e-2 1.01

     | Show Table
    DownLoad: CSV
    Table 4.  Spatial errors and convergence rates of and of the decoupled algorithm for .
    Rate Rate Rate Rate
    4.7651e-2 1.6055e-1 8.5700e-2 4.2574e-1
    5.6669e-3 3.07 4.4078e-2 1.86 1.7255e-2 2.31 2.0459e-1 1.06
    6.5948e-4 3.10 1.1349e-2 1.96 3.7745e-3 2.19 9.9988e-2 1.03
    7.9662e-5 3.05 2.8614e-3 1.99 8.7440e-4 2.11 4.9501e-2 1.01

     | Show Table
    DownLoad: CSV

    Thus, Tables 14 imply the coupled and decoupled algorithms are robust on the Poisson rate .

    Table 5 reports the CPU time for both the coupled and decoupled algorithms. The results indicate that the decoupled algorithm significantly outperforms the coupled one in terms of runtime, demonstrating its superior efficiency.

    Table 5.  The CPU time of coupled and decoupled algorithms for .
    the CPU time of coupled algorithm the CPU time of decoupled algorithm
    0.19 s 0.13 s
    2.61 s 1.63 s
    51.82 s 30.32 s
    915.23 s 503.59 s

     | Show Table
    DownLoad: CSV

    We now test the robustness of the algorithms with respect to the permeability tensor by fixing , , , , and for and e-, respectively.

    For the cases of and e-, Tables 6 and 7 exhibit the spatial errors and convergence rates of and for the coupled algorithm. The convergence rates of with respect to norm and norm are all around 3 and 2 separately. The convergence rates of with respect to norm and norm are all around 2 and 1, respectively. This implies that the coupled algorithm is robust for .

    Table 6.  Spatial errors and convergence rates of and of the coupled algorithm for .
    Rate Rate Rate Rate
    4.7630e-2 1.6047e-1 1.3174e-1 3.7941e-1
    5.6649e-3 3.07 4.4075e-2 1.86 3.5460e-2 1.89 1.9461e-1 0.96
    6.5936e-4 3.10 1.1349e-2 1.96 9.0399e-3 1.97 9.7954e-2 0.99
    7.9656e-5 3.05 2.8614e-3 1.99 2.2712e-3 1.99 4.9060e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 7.  Spatial errors and convergence rates of and of the coupled algorithm for e-.
    Rate Rate Rate Rate
    4.7630e-2 1.6047e-1 8.5753e-2 4.2449e-1
    5.6649e-3 3.07 4.4075e-2 1.86 1.7287e-2 2.31 2.0376e-1 1.06
    6.5936e-4 3.10 1.1349e-2 1.96 3.8309e-3 2.17 9.9511e-2 1.03
    7.9656e-5 3.05 2.8614e-3 1.99 1.0018e-3 1.93 4.9279e-2 1.01

     | Show Table
    DownLoad: CSV

    Again, for the cases of and e-, we compute the convergence rates of and by the decoupled algorithm. Tables 8 and 9 show that the convergence rates of and are optimal, which is consistent with the results in Theorem 5 and shows that the decoupled algorithm is robust for .

    Table 8.  Spatial errors and convergence rates of and of the decoupled algorithm for .
    Rate Rate Rate Rate
    4.7630e-2 1.6047e-1 1.3174e-1 3.7941e-1
    5.6649e-3 3.07 4.4075e-2 1.86 3.5460e-2 1.89 1.9461e-1 0.96
    6.5935e-4 3.10 1.1349e-2 1.96 9.0399e-3 1.97 9.7954e-2 0.99
    7.9656e-5 3.05 2.8614e-3 1.99 2.2712e-3 1.99 4.9060e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 9.  Spatial errors and convergence rates of and of the decoupled algorithm for e-.
    Rate Rate Rate Rate
    4.7630e-2 1.6047e-1 8.5753e-2 4.2449e-1
    5.6649e-3 3.07 4.4075e-2 1.86 1.7287e-2 2.31 2.0376e-1 1.06
    6.5935e-4 3.10 1.1349e-2 1.96 3.8309e-3 2.17 9.9511e-2 1.03
    7.9656e-5 3.05 2.8614e-3 1.99 1.0018e-3 1.94 4.9279e-2 1.01

     | Show Table
    DownLoad: CSV

    We now verify that there is no nonphysical oscillation in the pressure field for the proposed algorithms. For brevity, we only consider the decoupled algorithm. The benchmark problem for the poroelasticity model is given as follows:

    Example 2. (Barry-Mercer problem) Again we take and . The boundary segments are defined in Example 1. The source functions and , and the boundary and initial conditions are

    where

    Set , , e, , e-, and . For the case of , low permeability , and small time step , it will lead to for the discrete form of Eq (1.2). If using element pairs to solve the original problem, there exists the nonphysical oscillation in the pressure field at a very short time see Figure 1. But Figure 2 shows that there is no nonphysical oscillation in the pressure field for the decoupled algorithm with element pairs. In addition, the arrows near the boundary match very well with the inner arrows in Figure 3. This implies that the decoupled algorithm is locking-free and stable.

    Figure 1.  The numerical solution by using the element pairs for the original problem at .
    Figure 2.  The numerical solution by using the decoupled algorithm with element pairs at .
    Figure 3.  Arrow plot of numerical solution by using the decoupled algorithm with element pairs at .

    We here test the convergence rates for the BDF2 algorithm by fixing , , e, e- and for and , respectively.

    Tables 10 and 11 show the spatial errors and convergence rates of the BDF2 algorithm for and , respectively. The numerical results show that the convergence rates of and the convergence rates of are all approaching . The convergence rates of are around . The convergence rates of are around . This implies that the BDF2 algorithm has the optimal convergence rates for the different values of .

    Table 10.  Spatial errors and convergence rates of and of the BDF2 algorithm for .
    Rate Rate Rate Rate
    4.8482e-2 1.6013e-1 5.7104e-2 4.0384e-1
    5.7075e-3 3.09 4.4041e-2 1.86 1.3284e-2 2.10 1.9872e-1 1.02
    6.6098e-4 3.11 1.1349e-2 1.96 3.2745e-3 2.02 9.8508e-2 1.01
    7.9759e-5 3.05 2.8615e-3 1.99 9.3353e-4 1.81 4.9128e-2 1.00

     | Show Table
    DownLoad: CSV
    Table 11.  Spatial errors and convergence rates of and of the BDF2 algorithm for .
    Rate Rate Rate Rate
    4.8493e-2 1.6012e-1 5.7040e-2 4.0369e-1
    5.7096e-3 3.09 4.4040e-2 1.86 1.3273e-2 2.10 1.9874e-1 1.02
    6.6116e-4 3.11 1.1349e-2 1.96 3.2409e-3 2.03 9.8517e-2 1.01
    7.9771e-5 3.05 2.8615e-3 1.99 8.0517e-4 2.01 4.9130e-2 1.00

     | Show Table
    DownLoad: CSV

    In this section, the coupled and decoupled algorithms with element pairs and BDF2 algorithm with element pairs are used to simulate brain swelling caused by brain injury. Since the coupled and decoupled algorithms present almost the same result, we only report the simulation result based on the coupled algorithm. For more details about , source functions, boundary conditions, and the values of parameters, one can see [15]. Here we list the required information, which is obtained from [15]. The is a 2D cross-section of a 3D model of the brain (see Figure 4), and the values of parameters are shown in Table 12. The injured area is marked with a yellow line in the right-hand figure of Figure 4.

    Figure 4.  A slice of the magnetic resonance imaging for a human brain (left) and the finite-element mesh (right).
    Table 12.  Values of parameters.
    Parameters Values Parameters Values
    4.5e-7 K 1.4e-9
    3e-5 mm/min/Pa 1
    1070 Pa 0.35
    1.48e-5 E 9010 Pa

     | Show Table
    DownLoad: CSV

    Moreover, the source functions are e-, and the boundary conditions are

    Figures 5 and 6 show the simulation results of the coupled and BDF2 algorithms for an injured brain with respect to . From Figure 5, one can see that the maximum pressures of the coupled algorithm are Pa and Pa separately, which are consistent with that in [27]. And Figure 6 illustrates the similar facts for the BDF2 algorithm by the values of maximum pressures ( Pa and Pa). In addition, the brain tissue can deform due to the influence of the total stress. In our simulation, the values of maximum tissue deformation are mm and mm for the coupled algorithm and mm and mm for the BDF2 algorithm, respectively. It is easy to find that the change of has little influence on the maximum pressure, and the maximum tissue deformation becomes smaller as gets bigger. When the Poisson ratio is approaching , the brain tissue is nearly incompressible such that the tissue deformation is very small. Moreover, from Table 13, the execution time of the CPU of the coupled algorithm is longer than the BDF2 algorithm; thus the latter has higher computational efficiency.

    Figure 5.  Simulation results of the coupled algorithm with element pairs.
    Figure 6.  Simulation results of the BDF2 algorithm with element pairs.
    Table 13.  Execution time of CPU.
    Coupled algorithm BDF2 algorithm
    0.35 1547.65s 1239.13s
    0.496 126.97s 87.21s

     | Show Table
    DownLoad: CSV

    In this paper, we transform the poroelasticity model into a four-field problem by introducing the new variables and discuss the inf-sup condition about the four-field problem. Besides, we design the coupled and decoupled algorithms for the four-field problem and perform the error estimates. Optimal convergence is obtained in space and time. It can be observed that there is no constraint on mesh size and time step in error estimates. Then we introduce the BDF2 algorithm and perform stability analysis and error estimates. Furthermore, some numerical examples are presented to verify the theoretical results. More importantly, we successfully apply the algorithms to simulate the brain edema and observe the influence on the maximum tissue deformation and maximum pressure.

    Peizhen Wang: Conceptualization, investigation, writing-original draft; Wenlong He: Validation; Yajuan Di: Writing-review, editing. All authors have read and approved the final version of the manuscript for publication.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work was supported by the National Natural Science Foundation of China (Grant No.11901197) and the Young Teacher Foundation of Henan Province (Grant No.2021GGJS080). The numerical simulations have been done on the supercomputing system in the Supercomputing Center of Wuhan University.

    The authors declare that they have no competing interests.



    [1] J. Lian, S. Fang, Y. Zhou, Model predictive control of the fuel cell cathode system based on state quantity estimation, Comput. Simul., 37 (2020), 119-122.
    [2] H. M. Numanoğlu, H. Ersoy, B. Akgöz, O. Civalek, A new eigenvalue problem solver for thermos-mechanical vibration of Timoshenko nanobeams by an innovative nonlocal finite element method, Math. Methods Appl. Sci., 45 (2022), 2592-2614. https://doi.org/10.1002/mma.7942 doi: 10.1002/mma.7942
    [3] V. Chaturvedi, T. Talapaneni, Effect of mechanical vibration and grain refiner on microstructure and mechanical properties of AZ91Mg alloy during solidification, J. Mater. Eng. Perform., 30 (2021), 3187-3202. https://doi.org/10.1007/s11665-021-05471-3 doi: 10.1007/s11665-021-05471-3
    [4] W. Booyse, D. N. Wilke, S. Heyns, Deep digital twins for detection, diagnostics and prognostics, Mech. Syst. Signal Process., 140 (2020), 106612. https://doi.org/10.1016/j.ymssp.2019.106612 doi: 10.1016/j.ymssp.2019.106612
    [5] F. Tao, X. Sun, J. Cheng, Y. Zhu, W. Liu, Y. Wang, et al., 2023, MakeTwin: a reference architecture for digital twin software platform, Chin. J. Aeronaut., in press, 2023. https://doi.org/10.1016/j.cja.2023.05.002
    [6] Q. Qi, F. Tao, T. Hu, N. Anwer, A. Liu, Y. Wei, et al., Enabling technologies and tools for digital twin, J. Manuf. Syst., 58 (2021), 3-21. https://doi.org/10.1016/j.jmsy.2019.10.001 doi: 10.1016/j.jmsy.2019.10.001
    [7] S. Liu, Y. Lu, P. Zheng, H. Shen, J. Bao, Adaptive reconstruction of digital twins for machining systems: a transfer learning approach, Rob. Comput. Integr. Manuf., 78 (2022), 102390. https://doi.org/10.1016/j.rcim.2022.102390 doi: 10.1016/j.rcim.2022.102390
    [8] C. Gao, H. Park, A. Easwaran, An anomaly detection framework for digital twin driven cyber-physical systems, in ICCPS '21: Proceedings of the ACM/IEEE 12th International Conference on Cyber-Physical Systems, (2021), 44-54. https://doi.org/10.1145/3450267.3450533
    [9] X. Wang, Y. Wang, F. Tao, A. Liu, New paradigm of data-driven smart customisation through digital twin, J. Manuf. Syst., 58 (2021), 270-280. https://doi.org/10.1016/j.jmsy.2020.07.023 doi: 10.1016/j.jmsy.2020.07.023
    [10] V. S. Vishnu, K. G. Varghese, B. Gurumoorthy, A data-driven digital twin of CNC machining processes for predicting surface roughness, Procedia CIRP, 104 (2021), 1065-1070. https://doi.org/10.1016/j.procir.2021.11.179 doi: 10.1016/j.procir.2021.11.179
    [11] C. Zhang, G. Zhou, J. He, Z. Li, W. Cheng, A data- and knowledge-driven framework for digital twin manufacturing cell, Procedia CIRP, 83 (2019), 345-350. https://doi.org/10.1016/j.procir.2019.04.084 doi: 10.1016/j.procir.2019.04.084
    [12] Y. Sun, Y. Lu, J. Bao, F. Tao, Prognostics and health management via long short-term digital twins, J. Manuf. Syst., 68 (2023), 560-575. https://doi.org/10.1016/j.jmsy.2023.05.023 doi: 10.1016/j.jmsy.2023.05.023
    [13] K. Feng, J. C. Ji, Q. Ni, Y. Li, W. Mao, L. Liu, A novel vibration-based prognostic scheme for gear health management in surface wear progression of the intelligent manufacturing system, Wear, 522 (2023), 204697. https://doi.org/10.1016/j.wear.2023.204697 doi: 10.1016/j.wear.2023.204697
    [14] L. Ma, B. Jiang, L. Xiao, N. Lu, Digital twin-assisted enhanced meta-transfer learning for rolling bearing fault diagnosis, Mech. Syst. Signal Process., 200 (2023), 110490. https://doi.org/10.1016/j.ymssp.2023.110490 doi: 10.1016/j.ymssp.2023.110490
    [15] L. Li, Y. Ren, Q. Jin, Electro-mechanical vibration and stress field of piezoelectric nanobeam with symmetrical FGM core under the low-velocity impact, Eur. Phys. J. Plus, 137 (2022), 1-20. https://doi.org/10.1140/epjp/s13360-022-02934-x doi: 10.1140/epjp/s13360-022-02934-x
    [16] M. Rigacci, R. Sato, K. Shirase, Power consumption simulation of servo motors focusing on the influence of mechanical vibration on motor efficiency, Int. J. Autom. Technol., 16 (2022), 104-116. https://doi.org/10.20965/ijat.2022.p0104 doi: 10.20965/ijat.2022.p0104
    [17] P. Ewert, C. T. Kowalski, M. Jaworski, Comparison of the effectiveness of selected vibration signal analysis methods in the rotor unbalance detection of PMSM drive system, Electronics, 11 (2022), 1748. https://doi.org/10.3390/electronics11111748 doi: 10.3390/electronics11111748
    [18] Y. W. Zhang, G. L. She, Wave propagation and vibration of FG pipes conveying hot fluid, Steel Compos. Struct., 42 (2022), 397-405.
    [19] Y. Kumar, A. Gupta, A. Tounsi, Size-dependent vibration response of porous graded nanostructure with FEM and nonlocal continuum model, Adv. Nano Res., 11 (2021), 1-17. https://doi.org/10.12989/anr.2021.11.1.001 doi: 10.12989/anr.2021.11.1.001
    [20] S. K. Barman, M. Mishra, D. K. Maiti, D. Maity, Vibration-based damage detection of structures employing Bayesian data fusion coupled with TLBO optimization algorithm, Struct. Multidiscip. Optim., 64 (2021), 2243-2266. https://doi.org/10.1007/s00158-021-02980-6 doi: 10.1007/s00158-021-02980-6
    [21] F. L. Zhang, C. W. Kim, Y. Goi, Efficient Bayesian FFT method for damage detection using ambient vibration data with consideration of uncertainty, Struct. Control Health Monit., 28 (2021), e2659. https://doi.org/10.1002/stc.2659 doi: 10.1002/stc.2659
    [22] A. Turnbull, J. Carroll, A. McDonald, Combining SCADA and vibration data into a single anomaly detection model to predict wind turbine component failure, Wind Energy, 24 (2021), 197-211. https://doi.org/10.1002/we.2567 doi: 10.1002/we.2567
    [23] S. K. Barman, D. K. Maiti, D. Maity, Vibration-based delamination detection in composite structures employing mixed unified particle swarm optimization, AIAA J., 59 (2021), 386-399. https://doi.org/10.2514/1.J059176 doi: 10.2514/1.J059176
    [24] C. Tarawneh, J. Montalvo, B. Wilson, Defect detection in freight railcar tapered-roller bearings using vibration techniques, Railway Eng. Sci., 29 (2021), 42-58. https://doi.org/10.1007/s40534-020-00230-x doi: 10.1007/s40534-020-00230-x
    [25] Z. Mousavi, S. Varahram, M. M. Ettefagh, H. M. Sadeghi, N. S. Razavi, Deep neural networks-based damage detection using vibration signals of finite element model and real intact state: An evaluation via a lab-scale offshore jacket structure, Struct. Health Monit., 20 (2021), 379-405. https://doi.org/10.1177/1475921720932614 doi: 10.1177/1475921720932614
    [26] N. Wu, S. Haruyama, The 20k samples-per-second real time detection of acoustic vibration based on displacement estimation of one-dimensional laser speckle images, Sensors, 21 (2021), 2938. https://doi.org/10.3390/s21092938 doi: 10.3390/s21092938
    [27] M. H. M. Ghazali, W. Rahiman, Vibration-based fault detection in drone using artificial intelligence, IEEE Sensors J., 22 (2022), 8439-8448. https://doi.org/10.1109/JSEN.2022.3163401 doi: 10.1109/JSEN.2022.3163401
    [28] B. R. F. Rende, A. A. Cavalini, I. F. Santos, Fault detection using vibration data-driven models—a simple and well-controlled experimental example, J. Braz. Soc. Mech. Sci. Eng., 44 (2022), 1-11. https://doi.org/10.1007/s40430-022-03462-6 doi: 10.1007/s40430-022-03462-6
    [29] X. Huang, Q. Huang, H. Cao, W. Yan, L. Cao, Q. Zhang, Optimal design for improving operation performance of electric construction machinery collaborative system: Method and application, Energy, 263 (2023), 125629. https://doi.org/10.1016/j.energy.2022.125629 doi: 10.1016/j.energy.2022.125629
    [30] J. L. Conradi Hoffmann, L. P. Horstmann, M. Martínez Lucena, M. G. de Araujo, A. A. Fröhlich, H. M. Napoli Nishioka, Anomaly detection on wind turbines based on a deep learning analysis of vibration signals, Appl. Artif. Intell., 35 (2021), 893-913. https://doi.org/10.1080/08839514.2021.1966879
    [31] Y. Zhu, F. Li, W. Bao, Fatigue crack detection under the vibration condition based on ultrasonic guided waves, Struct. Health Monit., 20 (2021), 931-941. https://doi.org/10.1177/1475921719860772 doi: 10.1177/1475921719860772
  • This article has been cited by:

    1. Shishuo Fu, Jiaxi Lu, Yuanzhe Ding, A skeleton model to enumerate standard puzzle sequences, 2021, 30, 2688-1594, 179, 10.3934/era.2022010
    2. Kanasottu Anil Naik, Rayappa David Amar Raj, Chepuri Venkateswara Rao, Thanikanti Sudhakar Babu, Generalized cryptographic image processing approaches using integer-series transformation for solar power optimization under partial shading, 2022, 272, 01968904, 116376, 10.1016/j.enconman.2022.116376
    3. Sen-Peng Eu, Tung-Shan Fu, Springer Numbers and Arnold Families Revisited, 2024, 10, 2199-6792, 125, 10.1007/s40598-023-00230-9
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1618) PDF downloads(69) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog