Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article

A strategy for hepatitis diagnosis by using spherical q-linear Diophantine fuzzy Dombi aggregation information and the VIKOR method

  • Hepatitis is an infectious disease typified by inflammation in internal organ tissues, and it is caused by infection or inflammation of the liver. Hepatitis is often feared as a fatal illness, especially in developing countries, mostly due to contaminated water, poor sanitation, and risky blood transfusion practices. Although viruses are typically blamed, other potential causes of this kind of liver infection include autoimmune disorders, toxins, medicines, opioids, and alcohol. Viral hepatitis may be diagnosed using a variety of methods, including a physical exam, liver surgery (biopsy), imaging investigations like an ultrasound or CT scan, blood tests, a viral serology panel, a DNA test, and viral antibody testing. Our study proposes a new decision-support system for hepatitis diagnosis based on spherical q-linear Diophantine fuzzy sets (Sq-LDFS). Sq-LDFS form the generalized structure of all existing notions of fuzzy sets. Furthermore, a list of novel Einstein aggregation operators is developed under Sq-LDF information. Also, an improved VIKOR method is presented to address the uncertainty in analyzing the viral hepatitis categories demonstration. Interesting and useful properties of the proposed operators are given. The core of this research is the proposed algorithm based on the proposed Einstein aggregation operators and improved VIKOR approach to address uncertain information in decision support problems. Finally, a hepatitis diagnosis case study is examined to show how the suggested approach works in practice. Additionally, a comparison is provided to demonstrate the superiority and efficacy of the suggested decision technique.

    Citation: Huzaira Razzaque, Shahzaib Ashraf, Wajdi Kallel, Muhammad Naeem, Muhammad Sohail. A strategy for hepatitis diagnosis by using spherical q-linear Diophantine fuzzy Dombi aggregation information and the VIKOR method[J]. AIMS Mathematics, 2023, 8(6): 14362-14398. doi: 10.3934/math.2023735

    Related Papers:

    [1] Wenxiang Fang, Tao Xie, Biwen Li . Robustness analysis of fuzzy BAM cellular neural network with time-varying delays and stochastic disturbances. AIMS Mathematics, 2023, 8(4): 9365-9384. doi: 10.3934/math.2023471
    [2] N. Mohamed Thoiyab, Mostafa Fazly, R. Vadivel, Nallappan Gunasekaran . Stability analysis for bidirectional associative memory neural networks: A new global asymptotic approach. AIMS Mathematics, 2025, 10(2): 3910-3929. doi: 10.3934/math.2025182
    [3] R. Sriraman, R. Samidurai, V. C. Amritha, G. Rachakit, Prasanalakshmi Balaji . System decomposition-based stability criteria for Takagi-Sugeno fuzzy uncertain stochastic delayed neural networks in quaternion field. AIMS Mathematics, 2023, 8(5): 11589-11616. doi: 10.3934/math.2023587
    [4] Pratap Anbalagan, Evren Hincal, Raja Ramachandran, Dumitru Baleanu, Jinde Cao, Chuangxia Huang, Michal Niezabitowski . Delay-coupled fractional order complex Cohen-Grossberg neural networks under parameter uncertainty: Synchronization stability criteria. AIMS Mathematics, 2021, 6(3): 2844-2873. doi: 10.3934/math.2021172
    [5] Yunlong Ma, Tao Xie, Yijia Zhang . Robustness analysis of neutral fuzzy cellular neural networks with stochastic disturbances and time delays. AIMS Mathematics, 2024, 9(10): 29556-29572. doi: 10.3934/math.20241431
    [6] Tao Xie, Wenqing Zheng . Robustness analysis of Cohen-Grossberg neural network with piecewise constant argument and stochastic disturbances. AIMS Mathematics, 2024, 9(2): 3097-3125. doi: 10.3934/math.2024151
    [7] Ardak Kashkynbayev, Moldir Koptileuova, Alfarabi Issakhanov, Jinde Cao . Almost periodic solutions of fuzzy shunting inhibitory CNNs with delays. AIMS Mathematics, 2022, 7(7): 11813-11828. doi: 10.3934/math.2022659
    [8] R. Sriraman, P. Vignesh, V. C. Amritha, G. Rachakit, Prasanalakshmi Balaji . Direct quaternion method-based stability criteria for quaternion-valued Takagi-Sugeno fuzzy BAM delayed neural networks using quaternion-valued Wirtinger-based integral inequality. AIMS Mathematics, 2023, 8(5): 10486-10512. doi: 10.3934/math.2023532
    [9] Abdulaziz M. Alanazi, R. Sriraman, R. Gurusamy, S. Athithan, P. Vignesh, Zaid Bassfar, Adel R. Alharbi, Amer Aljaedi . System decomposition method-based global stability criteria for T-S fuzzy Clifford-valued delayed neural networks with impulses and leakage term. AIMS Mathematics, 2023, 8(7): 15166-15188. doi: 10.3934/math.2023774
    [10] Yijia Zhang, Tao Xie, Yunlong Ma . Robustness analysis of exponential stability of Cohen-Grossberg neural network with neutral terms. AIMS Mathematics, 2025, 10(3): 4938-4954. doi: 10.3934/math.2025226
  • Hepatitis is an infectious disease typified by inflammation in internal organ tissues, and it is caused by infection or inflammation of the liver. Hepatitis is often feared as a fatal illness, especially in developing countries, mostly due to contaminated water, poor sanitation, and risky blood transfusion practices. Although viruses are typically blamed, other potential causes of this kind of liver infection include autoimmune disorders, toxins, medicines, opioids, and alcohol. Viral hepatitis may be diagnosed using a variety of methods, including a physical exam, liver surgery (biopsy), imaging investigations like an ultrasound or CT scan, blood tests, a viral serology panel, a DNA test, and viral antibody testing. Our study proposes a new decision-support system for hepatitis diagnosis based on spherical q-linear Diophantine fuzzy sets (Sq-LDFS). Sq-LDFS form the generalized structure of all existing notions of fuzzy sets. Furthermore, a list of novel Einstein aggregation operators is developed under Sq-LDF information. Also, an improved VIKOR method is presented to address the uncertainty in analyzing the viral hepatitis categories demonstration. Interesting and useful properties of the proposed operators are given. The core of this research is the proposed algorithm based on the proposed Einstein aggregation operators and improved VIKOR approach to address uncertain information in decision support problems. Finally, a hepatitis diagnosis case study is examined to show how the suggested approach works in practice. Additionally, a comparison is provided to demonstrate the superiority and efficacy of the suggested decision technique.



    In this paper, we focus on the following hydrodynamically consistent Cahn-Hilliard phase field model for incompressible two-phase flows with variable density

    ρt+(ρu)=0, (1.1a)
    ϕt+uϕγΔw=0, (1.1b)
    w=Δϕ+1ε2(ϕ3ϕ), (1.1c)
    ρ(ut+uu)ηΔu+pλwϕ=0, (1.1d)
    u=0, (1.1e)

    where ϕ, w, u and p are the phase field function, the chemical potential, the velocity of flow and pressure respectively. Moreover, ρ=12(ρ1+ρ2)+ϕ2(ρ1ρ2) is the density, where ρ1 and ρ2 are the densities of the two fluids. The parameters γ, ε, η and λ are the mobility parameter related to the relaxation time scale, the interface thickness [1], the viscosity of the field, and the mixing energy density respectively. The system (1.1) is supplemented with appropriate boundary and initial conditions which are given as

    ϕ(x,0)=ϕ0(x), ρ(x,0)=ρ0(x), u(x,0)=u0(x),ϕ|Ω=0,         u|Ω=0.

    This model is also called the Cahn-Hilliard-Navier-Stokes model when the density is constant, which has many practical applications in physical and engineering, such as wetting, coating, and painting. By considering the influence of variable density, this model has broad applicability, which include highly stratified flows, interfaces between fluids of different densities and some problems of inertial confinement.

    The phase field method which has a wide range of applications is one of the main methods to deal with the fluid interface in two-phase flow modeling; see [2,3,4,5] and the references therein. It was initially developed to simulate solid-liquid phase transitions, where the interface is treated as a thin, smooth transition layer[6,7,8] to remove the singularities between two phases. The basic framework is to use phase field variables to represent the volume fraction of fluid components and then adopt a variational form to derive the model. Recently, it has increasingly attracted researchers' interest, mainly because the phase field method is superior to other available methods in some aspects of two-phase flow [9]. Various material properties or complex interface behaviors can be simulated directly by introducing suitable energy functions. Numerical solutions play a crucial role in their study and applications because the analytic solutions are usually not available.

    A few works have been devoted to the design, analysis and implementation of numerical schemes for the phase field model, although this model is very classical and canonical. We briefly review the available methods used in the phase field model. It is worth noting that there are projection/gauge/penalty methods [10,11,12], scalar auxiliary variables [13,14], linear stability [15,16], convex splitting [17,18,19], invariable energy quadratization(IEQ) [20,21], nonlinear quadratic [22], exponential time differencing [23], etc. In practical applications, we often couple the flow-field equation with the phase field equation. Typically two-phase incompressible flow models are coupled to phase-field models. The Cahn-Hilliard model is taken into account in this paper, because it is effective in the following two aspects: (ⅰ) Cahn-Hilliard models can accurately conserve the volume and dynamics; (ⅱ) the equation is one of the most important models in mathematical physics. Because of these reasons, here we use the Cahn-Hilliard equation developed in [24] to couple the incompressible flows with variable density.

    There also have been a lot of works on numerical approximates for the Cahn-Hilliard phase field model for incompressible two-phase flows with variable density. Hohenberg and Halperin proposed the model in [25] to simulate two incompressible viscous fluids with constant density. In [26], Gurtin et al. obtained the equal model by using the framework of rational continuum mechanics. A fully adaptive energy stabilization scheme is proposed in [27]. An efficient Picard iteration procedure was designed in [28] to further decouple the model. In the last years, many authors have been concerned with designing incompressible two-phase flow models with variable densities. Several efficient and energy-stable time discretization schemes for the coupled nonlinear Cahn-Hilliard phase field system with variable density are constructed; see [29]. Yang and Dong presented an energy-stable scheme in [30] for the numerical approximation of the two-phase governing equations with variable density and viscosity for the two fluids by introducing a scalar-valued variable related to the total of the kinetic energy and the potential free energy. A second-order accurate, coupled, energy-stable schemeis proposed in [31], where the Crank-Nicolson method and the IEQ method were used. In [32], the conservation scheme of the first-order energy law was established, in which the Cahn-Hilliard solver was used to decouple from the two-phase incompressible flows solver through the use of the fractional step method. Ye et al.[33] have designed a fully-decoupled type scheme to solve the Cahn-Hilliard phase field model for a two-phase incompressible fluid flow system with constant density. They only give a detailed practical implementation method and also prove the solvability. None of the various existing schemes have been subjected to error analysis, where the main difficulty lies in the delicate treatment of a several of nonlinear terms. Rigorous error estimates of models with variable densities, using an optimal order error bound, may seem to be a difficult prospect, but is a very interesting direction for future research.

    In this paper, we finally arrive at an unconditionally stable in energy, first-order time-accurate scheme for the incompressible Cahn-Hilliard two phase flows with variable density by coming up with a fractional step method. This method has an advantage over the projection method in that the original boundary conditions of the problem can be implemented in all substeps of the scheme. The popular approach to discretizing the Cahn-Hilliard phase field model (1.1b) and (1.1c) in time is based on the convex-splitting of the free energy functional, i.e., an idea that can be traced back to [32]. In the convex-splitting framework, one treats the contribution from the convex part implicitly and the contribution from the concave part explicitly. This treatment promotes the energy stability of the scheme and this property is unconditional in terms of time steps. We also give a rigorous proof of the convergence results and error estimates in the theoretical analysis. The main contribution of this paper is a rigorous error analysis, particularly under the condition that energy stability is available. To the best of the authors' knowledge, the proof developed in this article is the first to have the description of error estimation. The accuracy and stability are also demonstrated through the simulation of various numerical examples, where the challenge is in creating an efficient and easy to implement numerical scheme that preserves the energy dissipation law.

    The rest of this article is organized as follows. In Section 2, we construct an efficient time discrete scheme for variable density and derive unconditional energy stability. In Section 3, the error analysis of the semi-discrete scheme in time is provided. Some numerical experimentations are given in Section 4. Finally, conclusions are drawn in Section 5.

    For the sake of simplicity, some notations are needed for the following content. We assume that the domain ΩR2 is open, sufficiently smooth and bounded. For any two functions ϕ(x) and ψ(x), their L2 inner product on Ω is denoted by (ϕ,ψ)=Ωϕ(x)ψ(x)dx, and the L2 norm of ϕ(x) is denoted by ϕ=(ϕ,ϕ)12. Let τ>0 be the time step size and set tn=nτ for 0nN with T=Nτ. Moreover, we introduce the following spaces,

    H={uL2(Ω)|divu in Ω, un=0 on Ω},V=H10, V0={uV|divu=0 in Ω},M=L20(Ω)={qL2(Ω)|Ωqdx=0},

    where n is the outward normal vector of Ω. Next, we reformulate the system (1.1) as follows:

    ρn+1ρnτ+ρn+1un=0, (2.1a)
    ϕn+1ϕnτ+˜un+1ϕnγΔwn+1=0, (2.1b)
    wn+1=1ε2((ϕn+1)3ϕn)Δϕn+1, (2.1c)
    ρn˜un+1unτηΔ˜un+1+ρn+1(un)˜un+1λwn+1ϕn+14ρn+1(un)˜un+1=0, (2.1d)
    ρn+1un+1˜un+1τη(Δun+1Δ˜un+1)+pn+1=0, (2.1e)
    un+1=0 (2.1f)

    Remark 1. In (2.1d), the term 14ρn+1(un)˜un+1 is added to obtain the unconditional stability, as it is 0 if un=0.

    Theorem 1. (Stabilityofρ) For any τ>0 and any sequence {un}n=0,,N satisfying the boundary condition unn=0 on Ω, the solution {ρn}n=1,,N to (2.1a) satisfies

    ρN2+N1n=1ρn+1ρN2=ρ02. (2.2)

    Proof. Testing (2.1a) with 2τρn+1 gives

    ρn+12ρn2+ρn+1ρn2+2τΩ(unρn+1+12ρn+1un)ρn+1dx=0. (2.3)

    Owing to the boundary conditions on un, we note that

    Ω(unρn+1+12ρn+1un)ρn+1dx=12Ω(|ρn+1|2un)dx=12Ω|ρn+1|2unndx=0. (2.4)

    Thus, we get

    ρn+12+ρn+1ρn2=ρn2.

    Summing all indices n ranging 0 to N1, the proof is completed.

    Next, the stability of (2.1b)–(2.1e) will be proved in the theorem below. Moreover, since the kinetic energy of the fluid is 12ρnun2, it is more suitable to establish a bound based on ρnun2 than on the velocity itself. For simplicity, let us say that σn=ρn for all 1nN and σ0=ρ0.

    Theorem 2. (Stabilityofenergy) For any τ>0, (2.1b)–(2.1e) satisfy the the conditions of following energy estimates:

    σNuN2+λϕN2λε2ϕN2+λ2ε2ϕN4L4+2τγλN1n=0wn+12+λN1n=0((ϕn+1ϕn)2+12ε2ϕn+1ϕn2)+λ2ε2N1n=0(ϕn+12ϕn)2+N1n=0(σn(˜un+1un)2+σn+1(un+1˜un+1)2)+ητN1n=0(un+12+˜un+12+(un+1˜un+1)2)+λε2N1n=0(ϕn+1(ϕn+1ϕn)2)=σ0u02+λϕ02λε2ϕ02+λ2ε2ϕ04L4.

    Proof. Multiplying (2.1b) by 2τλwn+1 and integrating over Ω, we get

    2λ(ϕn+1ϕn,wn+1)+2τλγwn+12+2τλΩ˜un+1ϕnwn+1dx=0. (2.5)

    Multiplying (2.1c) by 2λ(ϕn+1ϕn) yields

    2λ(ϕn+1ϕn,wn+1)=λ2ε2(ϕn+14L4ϕn4L4)+λ2ε2(ϕn+12ϕn2)2+λε2ϕn+1(ϕn+1ϕn)2+λε2(ϕn2ϕn+12+ϕn+1ϕn2)+λ(ϕn+12ϕn2+(ϕn+1ϕn)2), (2.6)

    where we use the following identity:

    2a(ab)=a2b2+(ab)2,a3(ab)=14(a4b4+(a2b2)2+2a2(ab)2). (2.7)

    Testing (2.1a) with τ|˜un+1|2 leads to

    σn+1˜un+12σn˜un+12+τΩρn+1un|˜un+1|2dx+τ2Ωρn+1(un)|˜un+1|2dx=0. (2.8)

    By taking the inner product of (2.1d) with 2τ˜un+1, we have

    σn˜un+12σnun2+σn(˜un+1un)22τλΩwn+1ϕn˜un+1dx+2τη˜un+12+τΩρn+1(un)|˜un+1|2dx+τ2Ωρn+1(un)|˜un+1|2dx=0. (2.9)

    Testing (2.1e) with 2τun+1 yields

    σn+1un+12σn+1˜un+12+σn+1(un+1˜un+1)2+ητun+12+ητ˜un+12+ητ(un+1˜un+1)2=0. (2.10)

    Taking into account the boundary condition on un and using integration by parts, we have

    Ωρn+1un|˜un+1|2dx+Ωρn+1(un)|˜un+1|2dx+Ωρn+1un|˜un+1|2dx=Ω(ρn+1un|˜un+1|2)dx=Ωρn+1un|˜un+1|2ndx=0. (2.11)

    Summing the above inequality, we arrive at

    2τλγwn+12+λ2ε2(ϕn+14L4ϕn4L4)+λ2ε2(ϕn+12ϕn2)2+λε2ϕn+1(ϕn+1ϕn)2+λε2(ϕn2ϕn+12+ϕn+1ϕn2)+λ(ϕn+12ϕn2+(ϕn+1ϕn)2)+σn+1un+12σnun2+σn(˜un+1un)2+σn+1(un+1˜un+1)2+ητun+12+ητ˜un+12+ητ(un+1˜un+1)2=0. (2.12)

    Adding up the above inequality from n=0 to N1, we obtain Theorem 2.3.

    The bound on the pressure p is proved in the following theorem.

    Theorem 3. (Stabilityofp) For any τ>0, the solution pn+1 to (2.1e) satisfies the following inequality:

    τ2N1n=0pn+12C(σ0u02+λϕ02λε2ϕ02+λ2ε2ϕ04L4)(ρ0+1). (2.13)

    Proof. Under the inf-sup condition, there exists a positive constant β such that

    βpn+1sup (2.14)

    Testing (2.14) with all v \in V leads to

    \begin{equation} \begin{aligned} \left(\nabla \cdot v, p^{n+1}\right) = & \eta\left(\nabla\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right), \nabla v\right)+\tau^{-1}\left(\rho^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right), v\right) \\ \leq & \eta\left(\left\|\nabla \boldsymbol{u}^{n+1}\right\|+\left\|\nabla \tilde{\boldsymbol{u}}^{n+1}\right\|\right)\|\nabla v\| \\ &+\tau^{-1}\left\|\sigma^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right)\right\|_{L^{2}}\left\|\sigma^{n+1}\right\|_{L^{3}}\|v\|_{L^{6}}. \end{aligned} \end{equation} (2.15)

    Given that \left\|\rho^{n}\right\| \leq\left\|\rho_{0}\right\| for all 1 \leq n \leq N , and by using Hölder's inequality, we have

    \begin{equation} \|\sigma^{n+1}\|_{L^{3}} = (\int{\|p^{n+1}\|^{\frac{3}{2}}})^{\frac{1}{3}} \leq C \| \rho^{n+1} \|^{\frac{1}{2}} \leq C \| \rho^{0} \|^{\frac{1}{2}}. \end{equation} (2.16)

    Then, by the Sobolev embedding inequality \|v\|_{L^{6}} \leq C\|\nabla v\|_{L^{2}} for any v \in V , we get

    \begin{equation} \begin{aligned} \left(\nabla \cdot v, p^{n+1}\right) \leq & \eta\left(\left\|\nabla \boldsymbol{u}^{n+1}\right\|+\left\|\nabla \tilde{\boldsymbol{u}}^{n+1}\right\|\right)\|\nabla v\| \\ &+C \tau^{-1}\left\|\sigma^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right)\right\|\left\|\rho_{0}\right\|^{\frac{1}{2}}\|\nabla v\|. \end{aligned} \end{equation} (2.17)

    Substituting the above inequalities into Eq (2.14), we obtain

    \begin{equation} \begin{aligned} \beta\left\|p^{n+1}\right\| \leq & \eta\left(\left\|\nabla \boldsymbol{u}^{n+1}\right\|+\left\|\nabla \tilde{\boldsymbol{u}}^{n+1}\right\|\right) \\ &+C \tau^{-1}\left\|\sigma^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right)\right\|\left\|\rho_{0}\right\|^{\frac{1}{2}}. \end{aligned} \end{equation} (2.18)

    And by using Theorem 2.2, we get the desired result.

    In this section, we will give the time error estimates and show that the scheme has a first-order convergence rate. Although we verified that the scheme (2.1) is unconditionally stable in the previous chapter, we need to make the following assumptions[34,35] when conducting temporal error analysis:

    \begin{equation} \{\rho^n\}_{n = 0, \cdots, N} \ is\ uniformly\ bounded\ in \ L^\infty, \end{equation} (3.1)
    \begin{equation} for \ all\ n = 0, \cdots, N, \ {it\ holds\ that}\ \rho^n\geq\chi \ a.e. \ in \ \Omega , \end{equation} (3.2)

    where \chi is a number in (0, \rho_0^{min}] .

    We assume that the exact solution (\rho, \boldsymbol{u}, \phi, w, p) is sufficiently smooth. To be more precise,

    \begin{align} \rho &\in H^{2}\left(0, T ; L^{2}(\Omega)\right) \cap L^{2}\left(0, T ; W^{1, \infty}(\Omega)\right), \\ \phi &\in L^{\infty}\left(0, T ; H^{3}(\Omega)\right) \cap W^{1, \infty}\left(0, T ; H^{2}(\Omega)\right)\cap W^{2, \infty}\left(0, T ; H^{1}(\Omega)\right) \cap W^{3, \infty}\left(0, T ; L^{2}(\Omega)\right), \\ w& \in L^{\infty}\left(0, T ; H^{2}(\Omega)\right) \cap W^{1, \infty}\left(0, T ; L^{2}(\Omega)\right), \\ \boldsymbol{u} &\in \boldsymbol{H}^{2}\left(0, T ; \boldsymbol{L}^{2}(\Omega)\right) \cap \boldsymbol{L}^{\infty}\left(0, T ; \boldsymbol{V} \cap \boldsymbol{H}^{2}(\Omega)\right), \\ p& \in W^{2, \infty}\left(0, T ; H^{1}(\Omega)\right). \end{align} (3.3)

    We denote

    \begin{equation*} \begin{aligned} &e_{\phi}^{n} = \phi\left(t_{n}\right)-\phi^{n}, e_{w}^{n} = w\left(t_{n}\right)-w^{n}, e_{\rho}^{n} = \rho\left(t_{n}\right)-\rho^{n}, \\ &e_{u}^{n} = \boldsymbol{u} \left(t_{n}\right)-\boldsymbol{u} ^{n}, \tilde{e}_{u}^{n} = \boldsymbol{u} \left(t_{n}\right)-\tilde{\boldsymbol{u} }^{n}, e_{p}^{n} = p\left(t_{n}\right)-p^{n}.\\ \end{aligned} \end{equation*}

    In (1.1), taking t = t_{n+1} and subtracting from (2.1), we get the following error equations:

    \begin{align} &\frac{e_{\rho}^{n+1}-e_{\rho}^{n}}{\tau} +\nabla e_{\rho}^{n+1}\cdot\boldsymbol{u}(t_{n+1}) + \nabla\rho^{n+1}\cdot \left(\boldsymbol{u}(t_{n+1})-\boldsymbol{u}(t_{n})\right) +\nabla\rho^{n+1}\cdot e_{u}^{n} = R_\rho^{n+1}, \end{align} (3.4a)
    \begin{align} &\frac{e^{n+1}_\phi-e^{n}_\phi}{\tau} -\gamma\triangle e^{n+1}_w+\boldsymbol{u}(t_{n+1})\nabla\phi(t_{n+1})-\tilde{\boldsymbol{u}}^{n+1}\nabla\phi^{n} = R_\phi^{n+1}, \end{align} (3.4b)
    \begin{align} &e^{n+1}_w+\triangle e^{n+1}_\phi = \frac{1}{\varepsilon^2}(\phi(t_{n+1})^3-(\phi^{n+1})^3-\phi(t_{n+1})+\phi^n), \end{align} (3.4c)
    \begin{align} &\rho^{n} \frac{\tilde{e}_{u}^{n+1}-e_{u}^{n}}{\tau}-\eta \triangle \tilde{e}_{u}^{n+1}+\nabla p\left(t_{n+1}\right) +\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right) \\ &\qquad \quad -\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1} -\lambda w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)+ \lambda w^{n+1} \nabla \phi^{n} = R_{u}^{n+1}, \end{align} (3.4d)
    \begin{align} &\rho^{n+1} \frac{e_{u}^{n+1}-\tilde{e}_{u}^{n+1}}{\tau}-\eta\triangle \left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)-\nabla p^{n+1} = 0, \end{align} (3.4e)
    \begin{align} &\nabla \cdot e_{u}^{n+1} = 0, \end{align} (3.4f)

    where

    \begin{equation*} \begin{split} &R_{\phi}^{n+1} = \frac{\phi\left(t_{n+1}\right)-\phi\left(t_{n}\right)}{\tau}-\phi_{t}\left(t_{n+1}\right), \\ &R_{\rho}^{n+1} = \frac{\rho\left(t_{n+1}\right)-\rho\left(t_{n}\right)}{\tau}-\rho_{t}\left(t_{n+1}\right), \\ &R_{u}^{n+1} = \rho^{n} \frac{\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)}{\tau}-\rho\left(t_{n+1}\right) \boldsymbol{u}_{t}\left(t_{n+1}\right). \\ \end{split} \end{equation*}

    If the exact solution is sufficiently smooth, it is easy to establish the following estimate of the truncation error.

    Lemma 1. Under the regularity assumptions given by (3.3), the truncation errors satisfy:

    \begin{align*} \left\|R_{\rho}^{n+1}\right\|^{2}& \leq C \tau \int_{t_{n}}^{t_{n+1}}\left\|\rho_{t t}(t)\right\|^{2} d t \leq C \tau^{2}, \\ \left\|R_{\phi}^{n+1}\right\|^{2} &\leq C \tau \int_{t_{n}}^{t_{n+1}}\left\|\phi_{t t}(t)\right\|^{2} d t \leq C \tau^{2}, \\ \left\|R_{u}^{n+1}\right\|^{2} &\leq C \tau \int_{t_{n}}^{t_{n+1}}\left(\left\|\boldsymbol{u}_{t t}(t)\right\|^{2}+\left\|\rho_{t}(t)\right\|^{2}\right) d t+C\left\|e_{\rho}^{n}\right\|^{2} \\ &\leq C \tau^{2}+C\left\|e_{\rho}^{n}\right\|^{2}, \end{align*}

    for all 0 \leq n \leq N-1 .

    Proof. By using the integral residual of the Taylor formula, we have

    \begin{equation} \begin{split} R_\rho^{n+1} = \frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_n)\rho_{tt}(t)dt. \end{split} \end{equation} (3.5)

    By Hölder's inequality, we can derive

    \begin{equation} \begin{split} \Vert R_\rho^{n+1}\Vert^2 = & \int_\Omega(\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_n)\rho_{tt}(t)dt)^2dx\\ \leq&\frac{1}{\tau^2}\int_\Omega\int_{t_{n}}^{t_{n+1}}(t-t_n)^2dt^{\frac{1}{2}\cdot2}\int_{t_{n}}^{t_{n+1}}\rho_{tt}(t)^2dt^{\frac{1}{2}\cdot2}dx\\ \leq&\frac{1}{\tau^2}(\frac{\tau^3}{3}\int_{t_{n}}^{t_{n+1}}\int_\Omega\rho_{tt}(t)^2dxdt)\\ \leq&C\tau\int_{t_{n}}^{t_{n+1}}\Vert \rho_{tt}(t)\Vert^2 dt\\ \leq&C\tau^2.\\ \end{split} \end{equation} (3.6)

    Similarly, we can prove the inequality of R_\phi^{n+1} . For R_u^{n+1} , we can rewrite

    \begin{equation} \begin{split} R_u^{n+1}& = \rho^n(\frac{\boldsymbol{u}(t_{n+1})-\boldsymbol{u}(t_{n})}{\tau}-\boldsymbol{u}_t(t_{n+1})) -(e_\rho^n+\int_{t_{n}}^{t_{n+1}}\rho_{t}(t)dt)\boldsymbol{u}_t(t_{n+1})\\ & = \frac{\rho^n}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_n)\boldsymbol{u}_{tt}(t)dt-(e_\rho^n+\int_{t_{n}}^{t_{n+1}}\rho_{t}(t)dt)\boldsymbol{u}_t(t_{n+1}).\\ \end{split} \end{equation} (3.7)

    Using R_\rho^{n+1} estimation and Hölder's inequality can yield the result for R_u^{n+1} .

    We introduce the following Gronwall's inequality, which will frequently be used in error estimates.

    Lemma 2. Let a_{k} , b_{k} , c_{k} and \gamma_{k} , for integers k\geq0 , be the nonnegative numbers such that

    \begin{equation*} a_{n}+\tau\sum^{n}\limits_{k = 0}b_{k}\leq\tau\sum^{n}\limits_{k = 0}\gamma_{k}a_{k}+\tau\sum^{n}\limits_{k = 0}c_{k}+B\ \ for\ \ \geq0. \end{equation*}

    Suppose that \tau\gamma_{k} < 1 , for all k, and set \sigma_{k} = (1-\tau\gamma_{k})^{-1} . Then,

    \begin{equation*} a_{n}+\tau\sum^{n}\limits_{k = 0}b_{k}\leq \exp\bigg(\tau\sum^{n}\limits_{k = 0}\gamma_{k}\sigma_{k}\bigg)\bigg(\tau\sum^{n}\limits_{k = 0}c_{k}+B\bigg)\ \ for\ \ \geq0. \end{equation*}

    We verify Lemma 5 by the following lemma.

    Lemma 3. Define

    \begin{equation} G_{c}^{n+1} = \left( \phi(t_{n+1})\right)^3-(\phi^{n+1})^3 = 3\left( \phi(t_{n+1})\right)^2 e_{\phi}^{n+1}-3\phi(t_{n+1})\left( e_{\phi}^{n+1}\right)^2+\left( e_{\phi}^{n+1}\right)^3. \end{equation} (3.8)

    Then, for n < \frac{T}{\tau}-1 , we have

    \begin{equation} \begin{split} &\left\|G_{c}^{n+1}\right\| \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{1}}, \\ &\left\|e_{\phi}^{n+1}\right\|_{H^{2}}\leq C\left(\tau+\left\| e_w^{n+1} \right\|+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}+\left\| e_\phi^{n} \right\| \right).\\ \end{split} \end{equation} (3.9)

    Proof. For \left\|G_{c}^{n+1}\right\| , we use (3.8) to conclude that

    \begin{equation} \begin{split} \left\|G_{c}^{n+1}\right\| & \leq C\left(\left\|e_{\phi}^{n+1}\right\|\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}^{2}+\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{2}\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}+\left\|e_{\phi}^{n+1}\right\|_{L^{6}}^{3}\right) \\ & \leq C\left(\left\|e_{\phi}^{n+1}\right\|\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}^{2}+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}^{2}\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}^{3}\right) \\ & \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{1}}, \end{split} \end{equation} (3.10)

    where we have used the a priori bound \left\|\phi^{n}\right\|_{H^{1}} \leq C implied by the stability result given by Theorem 2. Using the H^{2} regularity results for elliptic equations, we conclude that

    \left\|e_{\phi}^{n+1}\right\|_{H^{2}} \leq c\left(\left\|e_{\phi}^{n+1}\right\|_{L^{2}}+\left\|\triangle e_{\phi}^{n+1}\right\|_{L^{2}}\right){;}

    from (3.4c), we know that

    \begin{equation*} \begin{split} \triangle e_{\phi}^{n+1} = -e_w^{n+1}+\frac{1}{\varepsilon^2}G_c^{n+1}-\frac{1}{\varepsilon^2}e_{\phi}^{n} -\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt{;} \end{split} \end{equation*}

    thus,

    \begin{equation*} \begin{split} \left\|e_{\phi}^{n+1}\right\|_{H^{2}}\leq& C\left( \left\|e_{\phi}^{n+1}\right\|+\left\| e_w^{n+1} \right\| +\left\|G_c^{n+1} \right\|+\left\| e_\phi^{n} \right\|+\tau\right)\\ \leq& C\left(\tau+\left\| e_w^{n+1} \right\|+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}+\left\| e_\phi^{n} \right\| \right).\\ \end{split} \end{equation*}

    The error estimate for the discrete density \rho^{n+1} is derived in the following lemma.

    Lemma 4. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) hold. Then, we have

    \begin{equation} \left\|e_{\rho}^{n+1}\right\|^{2}+2 \sum\limits_{m = 0}^{n}\left\|e_{\rho}^{m+1}-e_{\rho}^{m}\right\|^{2} \leq C\left(\tau^{2}+\tau \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}\right) \end{equation} (3.11)

    for all 0 \leq n \leq N-1 .

    Proof. Multiplying (3.4a) by 2 \tau e_{\rho}^{n+1} and integrating over \Omega , we have

    \begin{equation} \begin{split} &\left\|e_{\rho}^{n+1}\right\|^{2}-\left\|e_{\rho}^{n}\right\|^{2}+\left\|e_{\rho}^{n+1}-e_{\rho}^{n}\right\|^{2}\\ = &-2 \tau\left(\nabla e_{\rho}^{n+1} \cdot \boldsymbol{u}\left(t_{n+1}\right), e_{\rho}^{n+1}\right)+2 \tau\left(R_{\rho}^{n+1}, e_{\rho}^{n+1}\right)\\ &-2 \tau\left(\nabla \rho^{n+1} \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right)-2 \tau\left(\nabla \rho^{n+1} \cdot e_{u}^{n}, e_{\rho}^{n+1}\right)\\ = &\sum\limits_{i = 1}^{4} K_{i}.\\ \end{split} \end{equation} (3.12)

    Using \nabla \cdot \boldsymbol{u}\left(t_{n+1}\right) = 0 in \Omega and \boldsymbol{u}\left(t_{n+1}\right) = 0 on \partial \Omega , we have

    \begin{equation} K_{1} = \frac{1}{2} \int_{\partial \Omega}\vert e_{\rho}^{n+1}\vert^{2} \boldsymbol{u}\left(t_{n+1}\right) \cdot \boldsymbol{n} d s = 0 . \end{equation} (3.13)

    By using Young's inequality, the Cauchy-Schwarz inequality and Lemma 1, we have

    K_{2} \leq 2 \tau\left\|R_{\rho}^{n+1}\right\|\left\|e_{\rho}^{n+1}\right\| \leq C \tau\left\|R_{\rho}^{n+1}\right\|^{2}+\varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2} \leq C \tau^{3}+\varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2},

    where \varepsilon > 0 is a sufficiently small constant.

    Then, by the Sobolev inequality and Young inequality, we get

    \begin{equation} \begin{aligned} K_{3} & = -2 \tau\left(\nabla \rho\left(t_{n+1}\right) \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right)+2 \tau\left(\nabla e_{\rho}^{n+1} \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right) \\ & = -2 \tau\left(\nabla \rho\left(t_{n+1}\right) \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right) \\ & \leq 2 \tau\left\|\nabla \rho\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}} \boldsymbol{u}_{t}(t) d t\right\|\left\|e_{\rho}^{n+1}\right\| \\ & \leq \varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2}+C \tau^{2} \int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{t}(t)\right\|^{2} d t \\ & \leq C \tau^{3}+\varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2}. \end{aligned} \end{equation} (3.14)

    Similarly, the last term can be estimated as follows:

    \begin{equation} \begin{aligned} K_{4} & \leq 2 \tau\left\|\nabla \rho\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\frac{1}{\sigma^{n}}\right\|_{L^{\infty}}\|\| \sigma^{n} e_{u}^{n}\|\| e_{\rho}^{n+1} \| \\ & \leq \varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}. \end{aligned} \end{equation} (3.15)

    If \varepsilon is sufficiently small such that \varepsilon \tau \leq \frac{1}{6} , substituting the estimates of K_{i}\ (1 \leq i \leq 4) into (3.12), we have

    \begin{equation} \left\|e_{\rho}^{n+1}\right\|^{2}-\left\|e_{\rho}^{n}\right\|^{2}+2\left\|e_{\rho}^{n+1}-e_{\rho}^{n}\right\|^{2} \leq C \tau^{3}+C \tau\left\|e_{\rho}^{n}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}. \end{equation} (3.16)

    Using the discrete Gronwall inequality, we obtain the desired result.

    Lemma 5. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) are valid. For sufficiently small \tau , there are the following error estimates:

    \begin{equation} \begin{split} &\tau\gamma\sum\limits_{n = 0}^{N-1}\left(\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+\left\| e_{w}^{n+1} \right\|^2 \right)+\left\| e_{\phi}^{N}\right\|^2+\lambda\left\| \nabla e_{\phi}^{N}\right\|^2+\frac{\lambda}{2 \varepsilon^{2}}\left\|e_{\phi}^{N}\right\|_{L^{4}}^{4}\\ &+\left\|\sigma^{N} e_{u}^{N}\right\|^{2}+ \sum\limits_{n = 0}^{N-1}\left(\frac{1}{2}\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\right)\\ &+\tau\eta\sum\limits_{n = 0}^{N-1}\left(\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \frac{1}{2}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \right)\\ \leq& C\tau.\\ \end{split} \end{equation} (3.17)

    Proof. Let us multiply (3.4a) by \tau\vert\tilde{e}_{u}^{n+1}\vert^{2} , (3.4b) by 2 \tau e_{\phi}^{n+1} and 2 \tau \lambda e_{w}^{n+1} , (3.4c) by 2\tau\gamma e_{w}^{n+1} and -2 \lambda(e_{\phi}^{n+1}-e_{\phi}^{n}) , (3.4d) by 2 \tau \tilde{e}_{u}^{n+1} and (3.4e) by 2 \tau e_{u}^{n+1} . Summing up all of the above equations, we have

    \begin{align} &2\tau\gamma\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+ \left\| e_{\phi}^{n+1}\right\|^2- \left\| e_{\phi}^{n}\right\|^2+ \left\| e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &+2\tau\gamma\left\| e_{w}^{n+1} \right\|^2+\lambda\left\| \nabla e_{\phi}^{n+1}\right\|^2- \lambda\left\| \nabla e_{\phi}^{n}\right\|^2+ \lambda\left\| \nabla e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\\ &+ \tau\eta\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \tau\eta\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \tau\eta\left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \\ = &-2 \tau\left(\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) +2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1}, \tilde{e}_{u}^{n+1}\right)\\ &-\tau\left(\nabla \rho^{n+1} \cdot {\bf u}^{n}, \vert e_{u}^{n+1}\vert^{2}\right)-2 \tau\left(\nabla p\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right)+2 \tau\left(R_{u}^{n+1}, \tilde{e}_{u}^{n+1}\right) \\ &+2\tau\lambda \left( w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- w^{n+1} \nabla \phi^{n}, \tilde{e}_{u}^{n+1}\right) \\ &-2\tau\lambda \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)+ 2 \tau\lambda\left(R_{\phi}^{n+1}, e_{w}^{n+1}\right)\\ &-2\tau \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right) + 2 \tau\left(R_{\phi}^{n+1}, e_{\phi}^{n+1}\right)\\ &+\frac{2\tau\gamma}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)^3-\left(\phi^{n+1}\right)^3, e_{w}^{n+1}\right) -\frac{2\tau\gamma}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)-\phi^n, e_{w}^{n+1}\right)\\ &-\frac{2\lambda}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)^3-\left(\phi^{n+1}\right)^3, e_{\phi}^{n+1}-e_{\phi}^{n}\right) +\frac{2\lambda}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)-\phi^n, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ = &\sum\limits_{i = 1}^{i = 14}A_i.\\ \end{align} (3.18)

    Thanks to \nabla \cdot {\bf u}^{n} = 0 in \Omega , we have

    \begin{equation} \begin{aligned} A_{2}+A_{3} & = 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1}, \tilde{e}_{u}^{n+1}\right)-\tau\left(\nabla \rho^{n+1} \cdot \boldsymbol{u}^{n}, \vert\tilde{e}_{u}^{n+1}\vert^{2}\right) \\ & = 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \widetilde{e}_{u}^{n+1}\right)-\tau\left(\nabla \cdot\left(\boldsymbol{u}^{n} \rho^{n+1}\vert\tilde{e}_{u}^{n+1}\vert^{2}\right), 1\right) \\ & = 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right). \end{aligned} \end{equation} (3.19)

    Therefore, we get

    \begin{align} &A_{1}+A_{2}+A_{3} \end{align} (3.20)
    \begin{align} = &2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) \\ = &-2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right)-2 \tau\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(\rho\left(t_{n+1}\right)\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) \\ \leq& C \tau\left\|\rho^{n+1}\right\|_{L^{\infty}}\left\|\frac{1}{\sigma^{n}}\right\|_{L^{\infty}}\left\|\sigma^{n} e_{u}^{n}\right\|\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|\tilde{e}_{u}^{n+1}\right\|_{L^{6}} \end{align} (3.21)
    \begin{align} &+C \tau\left\|e_{\rho}^{n+1}\right\|\left\|\boldsymbol{u}\left(t_{n}\right)\right\|_{L^{\infty}} \left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|\tilde{e}_{u}^{n+1}\right\|_{L^{6}}\\ &+C \tau\left\|\rho\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}} \boldsymbol{u}_{t}(t) d t\right\|\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|\tilde{e}_{u}^{n+1}\right\|_{L^{6}} \\ \leq& C \tau\left(\left\|\sigma^{n} e_{u}^{n}\right\|+\left\|e_{\rho}^{n+1}\right\|+\tau \int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{t}(t)\right\| d t\right)\left\|\nabla \tilde{e}_{u}^{n+1}\right\| \\ \leq& \frac{\eta \tau}{12}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+C \tau\left(\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|e_{\rho}^{n+1}\right\|^{2}+\tau^{2}\right) \\ \leq& C \tau^{3}+\frac{\eta \tau}{12}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2} , \\ \end{align} (3.22)

    where the Sobolev embedding inequality \|v\|_{L^{6}} \leq C\|\nabla v\|_{L^{2}} for any v \in V is used.

    For A_4 , we have

    \begin{equation} \begin{split} A_4&\leq2\tau \left\|\nabla p\left(t_{n+1}\right)\right\|\left\|\frac{1}{\sigma^n}\right\|_{L^\infty} \left\|\sigma^n\left(\tilde e_u^{n+1}-e_u^{n} \right)\right\|\\ &\leq C\tau^2+\frac{1}{2} \left\| \sigma^n\left(\tilde e_u^{n+1}-e_u^{n} \right) \right\|^2 .\\ \end{split} \end{equation} (3.23)

    Using the Poincar\acute{e} inequality and Cauchy-Schwarz inequality, we can deal with A_6 , A_7 and A_9 :

    \begin{align} A_6 = &2\tau\lambda \left( w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- w^{n+1} \nabla \phi^{n}, \tilde{e}_{u}^{n+1}\right)\\ = &2\tau\lambda \left( w\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right), \tilde{e}_{u}^{n+1}\right) +2\tau\lambda \left( w\left(t_{n+1}\right)\nabla e_\phi^n, \tilde{e}_{u}^{n+1}\right)\\ &+2\tau\lambda \left( \nabla\phi^n e_w^{n+1}, \tilde{e}_{u}^{n+1}\right)\\ \leq& 2\tau\lambda \left\|w\left(t_{n+1}\right)\right\|_{L^\infty} \left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_tdt\right\|\left\|\tilde{e}_{u}^{n+1}\right\| +2\tau\lambda \left\|w\left(t_{n+1}\right)\right\|_{L^\infty} \left\| \nabla e_\phi^n \right\| \left\|\tilde{e}_{u}^{n+1} \right\|\\ &+2\tau\lambda \left( \nabla\phi^n e_w^{n+1}, \tilde{e}_{u}^{n+1}\right)\\ \leq& C\tau^3+\frac{\eta\tau}{12}\left\|\nabla\tilde{e}_{u}^{n+1}\right\|^2+C\tau\left\| \nabla e_\phi^n \right\|^2+2\tau\lambda \left( \nabla\phi^n e_w^{n+1}, \tilde{e}_{u}^{n+1}\right), \\ A_7 = &-2\tau\lambda \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)\\ = &-2\tau\lambda \left( \boldsymbol{u}\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right) +\boldsymbol{u}\left(t_{n+1}\right)\nabla e_\phi^n+\nabla\phi^n\tilde{e}_{u}^{n+1}, e_{w}^{n+1}\right)\\ \leq& 2\tau\lambda\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_tdt\right\| \left\|e_{w}^{n+1}\right\| +2\tau\lambda \left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty} \left\| \nabla e_\phi^n \right\| \left\|e_{w}^{n+1} \right\|\\ & -2\tau\lambda \left(\nabla\phi^n\tilde{e}_{u}^{n+1}, e_{w}^{n+1}\right)\\ \leq& C\tau^3+\frac{\tau\gamma}{3} \left\|e_w^{n+1}\right\|^2+C\tau \left\|\nabla e_\phi^{n}\right\|^2 -2\tau\lambda \left(\nabla\phi^n\tilde{e}_{u}^{n+1}, e_{w}^{n+1}\right), \end{align} (3.24)

    and

    \begin{equation} \begin{split} A_9 = &-2\tau \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right)\\ \leq& 2\tau \left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_tdt\right\| \left\|e_{\phi}^{n+1}\right\| +2\tau \left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty} \left\| \nabla e_\phi^n \right\| \left\|e_{\phi}^{n+1} \right\|\\ & +2\tau\left\|\nabla\phi^n \right\|\left\| \tilde{e}_{u}^{n+1} \right\|_{L_6} \left\| e_{\phi}^{n+1} \right\|_{L_6}\\ \leq& C\tau^3+C\tau \left\|\nabla e_\phi^{n+1}\right\|^2+C\tau \left\|\nabla e_\phi^{n}\right\|^2 +\frac{\eta\tau}{12} \left\|\nabla\tilde{e}_{u}^{n+1}\right\|^2.\\ \end{split} \end{equation} (3.25)

    From Lemma 1, we have

    \begin{equation} \begin{split} A_5+A_8+A_{10} \leq& 2\tau \left\|R_u^{n+1} \right\| \left\| \tilde{e}_{u}^{n+1} \right\| +2\tau\lambda\left\|R_\phi^{n+1} \right\| \left\| e_w^{n+1} \right\| +2\tau \left\| R_\phi^{n+1} \right\| \left\| e_\phi^{n+1} \right\|\\ \leq& C\tau^3+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2} +\frac{\eta\tau}{12} \left\|\nabla \tilde{e}_{u}^{n+1} \right\|^2+\frac{\tau\gamma}{3} \left\| e_w^{n+1} \right\|^2 +\frac{1}{2} \left\| e_\phi^{n+1} \right\|^2.\\ \end{split} \end{equation} (3.26)

    From Lemma 3, we obtain

    \begin{equation} \begin{aligned} A_{11}+A_{12} & = \frac{2\tau\gamma}{\varepsilon^2}\left(G_{c}^{n+1}-e_{\phi}^n-\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{w}^{n+1}\right) \\ & \leq \frac{2\tau\gamma}{\varepsilon^2}\left(\left\|e_{\phi}^n\right\|+\left\|G_{c}^{n+1}\right\|+ \left\|\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt\right\|\right)\left\|e_{w}^{n+1}\right\| \\ & \leq C\tau^3+C \tau\left(\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}\right) +\frac{\tau\gamma}{3}\left\|e_{w}^{n+1}\right\|^{2}, \\ \end{aligned} \end{equation} (3.27)

    and

    \begin{align} A_{13}+A_{14} = &- \frac{2\lambda}{\varepsilon^2}\left(G_{c}^{n+1}-e_{\phi}^n-\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{\phi}^{n+1}-e_{\phi}^{n}\right) \\ = &-\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|\left(e_{\phi}^{n+1}\right)^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e_{\phi}^{n}\right)\right\|^{2}\right)-\frac{2\lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2} \\ &-\frac{2\lambda}{\varepsilon^{2}}\left(3\left(\phi\left(t_{n+1}\right)\right)^{2} e_{\phi}^{n+1}-3 \phi\left(t_{n+1}\right)\left(e_{\phi}^{n+1}\right)^{2}-e_{\phi}^{n+1}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ &+\frac{2\lambda}{\varepsilon^{2}} \left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{\phi}^{n+1}-e_{\phi}^{n}\right), \end{align} (3.28)

    where we use this identity

    \begin{equation} \begin{split} a^{3}(a-b) = \frac{1}{4}\left(a^{4}-b^{4}+\left(a^{2}-b^{2}\right)^{2}+2 a^{2}(a-b)^{2}\right). \end{split} \end{equation} (3.29)

    We denote

    \begin{equation} \tilde{G}^{n+1} = -\frac{2\lambda}{\varepsilon^{2}}\left(3\left(\phi\left(t_{n+1}\right)\right)^{2} e_{\phi}^{n+1}-3 \phi\left(t_{n+1}\right)\left(e_{\phi}^{n+1}\right)^{2}-e_{\phi}^{n+1}\right) \text {. } \end{equation} (3.30)

    Similar to the method estimated from Lemma 3, we can obtain

    \begin{equation} \left\|\tilde{G}^{n+1}\right\| \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{1}} . \end{equation} (3.31)

    Taking the gradient of \tilde{G}^{n+1} , we get

    \begin{equation} \begin{aligned} \nabla \tilde{G}^{n+1} = &-\frac{2\lambda}{\varepsilon^{2}}\left[\left(3\left(\phi\left(t_{n+1}\right)\right)^{2}-1\right) \nabla e_{\phi}^{n+1}+6 \phi\left(t_{n+1}\right) e_{\phi}^{n+1} \nabla \phi\left(t_{n+1}\right)\right.\\ &\left.-3\left(e_{\phi}^{n+1}\right)^{2} \nabla \phi\left(t_{n+1}\right)-6 \phi\left(t_{n+1}\right) e_{\phi}^{n+1} \nabla e_{\phi}^{n+1}\right]. \end{aligned} \end{equation} (3.32)

    Since H^{2}(\Omega) \subset L^{\infty}(\Omega) and by utilizing the bound of \left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}} implied by Theorem 2, we conclude that

    \begin{equation} \left\|e_{\phi}^{n+1} \nabla e_{\phi}^{n+1}\right\|_{L^{2}} \leq\left\|e_{\phi}^{n+1}\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}} \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{2}} . \end{equation} (3.33)

    In view of (3.10) and the bound \left\|\phi^{n}\right\|_{H^{1}} < C , we have

    \begin{equation} \begin{array}{rl} \left\|\nabla \tilde{G}^{n+1}\right\| \leq & C\left[\left(\left\|\phi\left(t_{n+1}\right)\right\|_{L^{\infty}}^{2}+1\right)\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}}\right.\\ &+\left\|\phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{3}}\left\|e_{\phi}^{n+1}\right\|_{L^{6}} \\ & \left.+\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{6}}\left\|e_{\phi}^{n+1}\right\|_{L^{6}}^{2}+\left\|\phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|e_{\phi}^{n+1} \nabla e_{\phi}^{n+1}\right\|_{L^{2}}\right] \\ \leq &C\left(\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}}+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}\right. \left.+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}^{2}+\left\| e_{\phi}^{n+1} \nabla e_{\phi}^{n+1} \right\|_{L^{2}}\right) \\ \leq &C\left(\tau+\left\|e_{w}^{n+1}\right\|+\left\|e_{\phi}^{n+1}\right\|+\left\|\nabla e_{\phi}^{n+1}\right\| +\left\|e_{\phi}^{n}\right\| \right) . \end{array} \end{equation} (3.34)

    Dealing with the penultimate term of (3.28) and in view of (3.31) and (3.33), we get

    \begin{equation} \begin{split} &\left(\tilde{G}^{n+1}, e_{\phi}^{n+1}-e_{\phi}^{n}\right) \\ = &\tau\left(\tilde{G}^{n+1}, \gamma \Delta e_{w}^{n+1}-\boldsymbol{u}\left(t_{n+1}\right) \nabla e_{\phi}^{n}-\tilde{e}_{u}^{n+1} \nabla \phi^{n}\right)\\ &+\tau\left(\tilde{G}^{n+1}, -\boldsymbol{u}\left(t_{n+1}\right)\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt+R_{\phi}^{n+1}\right)\\ \leq& \gamma\tau\left\|\nabla e_{w}^{n+1}\right\|\left\|\nabla \tilde{G}^{n+1}\right\| +\tau\left\| \nabla\phi^n\right\|\left\| \tilde{e}_u^{n+1}\right\|_{H_1}\left\| \tilde{G}^{n+1} \right\|_{H_1}\\ &+\tau\left\|\tilde{G}^{n+1}\right\|\left(\left\|R_{\phi}^{n+1}\right\|+\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\|+\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt\right\| \right) \\ \leq& C\tau^3+\frac{\gamma\lambda\tau}{2}\left\|\nabla e_{w}^{n+1}\right\|^2+\frac{\eta\tau}{12}\left\| \nabla\tilde{e}_u^{n+1}\right\|^2\\ &+C\tau\left(\left\|e_{w}^{n+1}\right\|^2+\left\|e_{\phi}^{n+1}\right\|^2+\left\|\nabla e_{\phi}^{n+1}\right\|^2 +\left\|e_{\phi}^{n}\right\| ^2+\left\|\nabla e_{\phi}^{n}\right\| ^2 \right).\\ \end{split} \end{equation} (3.35)

    Then, we estimate the last term of (3.28),

    \begin{align} & \frac{2\lambda}{\varepsilon^2}\left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{\phi}^{n+1}-e_{\phi}^{n}\right) \\ = &\frac{2\lambda\tau}{\varepsilon^2}\left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, \gamma \Delta e_{w}^{n+1}-\boldsymbol{u}\left(t_{n+1}\right) \nabla e_{\phi}^{n}-\tilde{e}_{u}^{n+1} \nabla \phi^{n}\right)\\ &+\frac{2\lambda\tau}{\varepsilon^2}\left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, -\boldsymbol{u}\left(t_{n+1}\right)\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt+R_{\phi}^{n+1}\right)\\ \leq& \frac{2\lambda\gamma\tau}{\varepsilon^2}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt\right\|\left\|\nabla e_w^{n+1}\right\| +\frac{2\lambda\tau}{\varepsilon^2}\left\| \nabla\phi^n\right\|\left\| \tilde{e}_u^{n+1}\right\|_{H_1}\left\| \int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt \right\|_{H_1}\\ &+\frac{2\lambda\tau}{\varepsilon^2}\left\|\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt\right\|\left(\left\|R_{\phi}^{n+1}\right\|+\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\| \right) \\ &+\frac{2\lambda\tau}{\varepsilon^2}\left\|\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt\right\|\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt\right\| \\ \leq& C\tau^3+\frac{\gamma\lambda\tau}{2}\left\|\nabla e_{w}^{n+1}\right\|^2+\frac{\eta\tau}{12}\left\|\nabla \tilde{e}_{u}^{n+1}\right\| ^2 +C\tau\left\|\nabla e_{\phi}^{n}\right\| ^2 .\\ \end{align} (3.36)

    Combining (3.35) and (3.36), we can obtain

    \begin{equation} \begin{split} A_{13}+A_{14} \leq&-\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|\left(e_{\phi}^{n+1}\right)^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e_{\phi}^{n}\right)\right\|^{2}\right)-\frac{2\lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2} \\ &+ C\tau^3+\gamma\lambda\tau\left\|\nabla e_{w}^{n+1}\right\|^2+\frac{\eta\tau}{6}\left\|\nabla \tilde{e}_{u}^{n+1}\right\| ^2\\ &+C\tau \left( \left\|\nabla e_{\phi}^{n}\right\| ^2+\left\|e_{w}^{n+1}\right\|^2+\left\|e_{\phi}^{n+1}\right\|^2+\left\|\nabla e_{\phi}^{n+1}\right\|^2 +\left\|e_{\phi}^{n}\right\| ^2+\left\|\nabla e_{\phi}^{n}\right\| ^2 \right) .\\ \end{split} \end{equation} (3.37)

    Substituting the above estimates into (3.18), we have

    \begin{equation} \begin{split} &\tau\gamma\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+ \frac{1}{2}\left\| e_{\phi}^{n+1}\right\|^2- \left\| e_{\phi}^{n}\right\|^2+ (1+\frac{2\lambda}{\varepsilon^{2}})\left\| e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &+\tau\gamma\left\| e_{w}^{n+1} \right\|^2+\lambda\left\| \nabla e_{\phi}^{n+1}\right\|^2- \lambda\left\| \nabla e_{\phi}^{n}\right\|^2+ \lambda\left\| \nabla e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &+\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\frac{1}{2}\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\\ &+ \tau\eta\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \frac{\tau\eta}{2}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \tau\eta\left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \\ &+\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|\left(e_{\phi}^{n+1}\right)^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e_{\phi}^{n}\right)\right\|^{2}\right)\\ \leq& C \tau^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}\\ &+C\tau \left( \left\|\nabla e_{\phi}^{n}\right\| ^2+\left\|e_{w}^{n+1}\right\|^2+\left\|e_{\phi}^{n+1}\right\|^2+\left\|\nabla e_{\phi}^{n+1}\right\|^2 +\left\|e_{\phi}^{n}\right\| ^2+\left\|\nabla e_{\phi}^{n}\right\| ^2 \right).\\ \end{split} \end{equation} (3.38)

    Adding up from 0 to N-1 , and applying Gronwall's inequality, we infer that

    \begin{align*} \label{eq12} &\tau\gamma\sum\limits_{n = 0}^{N-1}\left(\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+\left\| e_{w}^{n+1} \right\|^2 \right)+\left\| e_{\phi}^{N}\right\|^2+\lambda\left\| \nabla e_{\phi}^{N}\right\|^2+\frac{\lambda}{2 \varepsilon^{2}}\left\|e_{\phi}^{N}\right\|_{L^{4}}^{4}\\ &+\left\|\sigma^{N} e_{u}^{N}\right\|^{2}+ \sum\limits_{n = 0}^{N-1}\left(\frac{1}{2}\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\right)\\ &+\tau\eta\sum\limits_{n = 0}^{N-1}\left(\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \frac{1}{2}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \right)\\ \leq& C\tau. \end{align*}

    Lemma 5 shows that the fractional scheme converges at a rate of O(\tau^{\frac{1}{2}}) . However, since we use the first-order backward Euler method of time discretization, it is not optimal from the perspective of theoretical analysis. Next, we will improve the convergence speed to the first order.

    Lemma 6. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) are valid. For sufficiently small \tau , there are the following error estimates:

    \begin{equation} \begin{split} &\tau\gamma\sum\limits_{n = 0}^{N-1}\left(\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+\left\| e_{w}^{n+1} \right\|^2 \right)+\left\| e_{\phi}^{N}\right\|^2+\lambda\left\| \nabla e_{\phi}^{N}\right\|^2\\ &+\frac{\lambda}{2 \varepsilon^{2}}\left\|e_{\phi}^{N}\right\|_{L^{4}}^{4} +\left\|\sigma^{N} e_{u}^{N}\right\|^{2}+\tau\eta\sum\limits_{n = 0}^{N-1}\left\|\nabla e_{u}^{n+1}\right\|^{2}\\ \leq& C\tau^2.\\ \end{split} \end{equation} (3.39)

    Proof. Taking the sum of (3.4d) and (3.4e), we get

    \begin{equation} \begin{split} &\rho^{n} \frac{e_u^{n+1}-e_u^ n}{\tau}+\frac{\rho^{n+1}-\rho^{n}}{\tau}\left(e_{u}^{n+1}-\tilde{e}_{n}^{n+1}\right)-\eta \triangle e_{u}^{n+1}+\rho\left(t_{n+1}\right)\left(u\left(t_{n+1}\right) \cdot \nabla\right) u(t_{n+1})\\ &-\rho^{n+1}\left(u^{n} \cdot \nabla\right) \tilde{u}^{n+1}+\Delta e_{p}^{n+1}-\lambda w\left(t _{n+1}\right) \nabla \phi(t_{n+1})+\lambda w^{n+1}\nabla\phi^{n} = R_{u}^{n+1}.\\ \end{split} \end{equation} (3.40)

    Let us multiply (3.4a) by \tau\tilde{e}_{u}^{n+1} , (3.4b) by 2 \tau e_{\phi}^{n+1} and 2 \tau \lambda e_{w}^{n+1} , (3.4c) by 2\tau\gamma e_{w}^{n+1} and -2 \lambda(e_{\phi}^{n+1}-e_{\phi}^{n}) and (3.40) by 2 \tau e_{u}^{n+1} . Summing up all of the above equations, we obtain

    \begin{align} &\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|\sigma^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+ 2 \eta\tau\left\|\nabla e_{u}^{n+1}\right\|^{2}+2 \tau \gamma \lambda\left\|\nabla e_{w}^{n+1}\right\|^{2}\\ &+\left\|e_{\phi}^{n+1}\right\|^{2}-\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2}+2 \tau \gamma\left\|e_{w}^{n+1}\right\|^{2}+\lambda\left\|\nabla e_{\phi}^{n+1}\right\|^{2}-\lambda\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+\lambda\left\|\nabla e _{\phi}^{n+1}-\nabla e_{\phi}^{n} \right\|^{2}\\ = &-2\tau\left(\nabla\rho^{n+1}\cdot u^n, \tilde{e}_{u}^{n+1}\cdot e_{u}^{n+1} \right)+\tau\left( \nabla\rho^{n+1}\cdot u^n, \vert e_{u}^{n+1}\vert^2 \right)\\ &-2 \tau\left(\rho\left(t_{n+1}\right)\left(u\left(t_{n+1}\right)\cdot\nabla\right) u\left(t_{n+1}\right)-\rho^{n+1}\left(u^{n}\cdot\nabla\right) \tilde{u}^{n+1}, e_u^{n+1}\right)\\ &+2 \tau \lambda\left(w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-w^{n+1} \nabla \phi^{n}, e_{u}^{n+1}\right)+2 \tau\left(R_{u}^{n+1}, e_u^{n+1}\right)\\ &-2 \tau \lambda\left(u\left(t_{n+1}\right) \nabla \phi\left(t _{n+1}\right)-\tilde{u}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)+2 \tau \lambda \left(R_{\phi}^{n+1}, e_w^{n+1}\right)\\ &-2 \tau\left(u\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-\tilde{u}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right)+2 \tau\left(R_{\phi}^{n+1}, e_{\phi}^{n+1}\right)\\ &+\frac{2 \tau r}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)^{3}-\left(\phi^{n+1}\right)^{3}, e_w^{n+1}\right)-\frac{2 \tau r}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)-\phi^{n}, e_{w}^{n+1}\right)\\ &-\frac{2 \lambda}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)^{3}-\left(\phi^{n+1}\right)^{3}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)+\frac{2 \lambda}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)-\phi^{n}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ = &\sum\limits_{i = 1}^{13} L_{i}. \end{align} (3.41)

    Using integration by parts, we have

    \begin{equation} \begin{split} L_{1}+L_{2} = & 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{e}_{u}^{n+1}, e_{ u}^{n+1}\right)+2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}\left(t_{n}\right) \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right){;} \end{split} \end{equation} (3.42)

    we rewrite the L_{3} as

    \begin{equation} \begin{split} L_{3} = &-2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde e_{u}^{n+1}, e_{u}^{n+1}\right)-2 \tau\left(\rho^{n+1}\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right)-2 \tau\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right){;} \end{split} \end{equation} (3.43)

    then,

    \begin{equation} \begin{split} L_{1}+L_{2}+L_{3} = & 2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}\left(t_{n}\right) \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ = & \sum\limits_{i = 1}^{5} J_{i}. \\ \end{split} \end{equation} (3.44)

    Lemma 5 shows that

    \begin{equation} \begin{split} \left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\| \leq C.\quad\quad \end{split} \end{equation} (3.45)

    According to Young's inequality and the Cauchy-Schwarz inequality, we can deduce that

    \begin{align} J_{1} & \leq \frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{u}^{n}\right\|^{2}\left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\| \\ & \leq \frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau^{\frac{3}{2}}\left\|\nabla e_{u}^{n}\right\|^{2}, \\ J_{2} & \leq 2 \tau\left\|e_{\rho}^{n+1}\right\|\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|e_{u}^{n+1}\right\|_{L^{6}} \\ & \leq \frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|e_{\rho}^{n+1}\right\|^{2} \\ & \leq C \tau^{3}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}, \\ J_{3} & \leq C \tau\left\|\nabla e_{u}^{n+1}\right\|\left\|\boldsymbol{u}\left(t_{n}\right)\right\|_{{L}^{\infty}}\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\| \\ & \leq C \tau^{3}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2, \\ J_{4} & \leq 2 \tau\left\|\rho^{n+1}\right\|_{{L}^{\infty}}\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}} \boldsymbol{u}_{t} d t\right\|\left\|e_{u}^{n+1}\right\| \\ & \leq C \tau^{3}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}, \\ J_{5} & \leq 2 \tau\left\|\rho^{n+1}\right\|_{L^{\infty}}\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\frac{1}{\sigma^{n}}\right\|_{L^{\infty}}\left\|\sigma^{n} e_{u}^{n}\right\|\left\|e_{u}^{n+1}\right\| \\ & \leq C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}. \end{align} (3.46)

    Therefore, we have

    \begin{equation} \begin{split} M_{1}+M_{2}+M_{3} \leq& C \tau^{3}+C \tau^{\frac{3}{2}}\left\|\nabla e_{u}^{n}\right\|^{2}+\frac{\eta \tau}{2}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}\\ &+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2. \end{split} \end{equation} (3.47)

    Using the embedding inequality and Cauchy-Schwarz inequality, we can deal with M_4 M_9

    \begin{equation} \begin{split} M_{4} = &2 \tau \lambda\left(w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-w^{n+1} \nabla \phi^{n}, e_{u}^{n+1}\right)\\ = &2 \tau \lambda\left(w\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right)+w\left(t_{n+1}\right) \nabla e_{\phi}^{n}+\nabla \phi^{n} e_{w}^{n+1}, e_{u}^{n+1}\right)\\ \leq& 2 \tau \lambda \left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\| \int_{t_{n}}^{t_{n+1}} \nabla \phi_{t}(t) d t\right\|\| e_{u}^{n+1}\|+2 \tau \lambda\left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\|\left\|e_{u}^{n+1}\right\|\\ &+2 \tau \lambda\left\|\nabla \phi^{n}\right\|_{L^2}\left\|e_w^{n+1}\right\|_{L^{3}}\left\|e_u^{n+1}\right\|_{L^{6}}\\ \leq& C \tau^{3}+\tau \eta\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}+C \tau\left\|\nabla e_{w}^{n+1}\right\|^{2}, \\ \end{split} \end{equation} (3.48)
    \begin{equation} \begin{split} M_{5}+M_{7}+M_{9}& = 2 \tau\left(R_{u}^{n+1}, e_{u}^{n+1}\right)+2\tau \lambda\left(R_{\phi}^{n+1}, e_{w}^{n+1}\right)+2\tau\left(R_{\phi}^{n+1}, e_{\phi}^{n+1}\right)\\ &\leq2\tau\left\|R_{u}^{n+1}\right\|\left\|e_{u}^{n+1}\right\|+2 \tau \lambda\left\|R_{\phi}^{n+1}\right\|\left\|e_{w}^{n+1}\right\|+2\tau\left\|R_{\phi}^{n+1}\right\|\left\|e_{\phi}^{n+1}\right\|\\ &\leq C \tau^{3}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}+\frac{\tau\lambda\gamma}{4}\left\|\nabla e_w^{n+1}\right\|^{2}+\frac{1}{4}\| e_\phi^{n+1} \|^{2}, \end{split} \end{equation} (3.49)
    \begin{equation} \begin{split} M_{6} = &-2 \tau \lambda\left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-\tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)\\ = &-2 \tau \lambda\left(\left(\tilde{e}_{u}^{n+1}-e_{u}^{n+1}\right) \nabla \phi\left(t_{n+1}\right), e_w^{n+1}\right)-2 \tau\lambda\left(e_{u}^{n+1} \nabla \phi\left(t_{n+1}\right), e_{w}^{n+1}\right)\\ &-2 \tau \lambda\left(\tilde{\boldsymbol{u}}^{n+1} \int_{t_n}^{t_{n+1}} \nabla \phi_t(t) d t, e_{w}^{n+1}\right)-2 \tau \lambda \left(\tilde{\boldsymbol{u}}^{n+1} \nabla e_{\phi}^{n}, e_{w}^{n+1}\right)\\ \leq& 2 \tau \lambda\left\|\frac{1}{\sigma^{n+1}}\right\|_{L^{\infty}}\left\|\sigma^{n+1}\left(\tilde{e}_{u}^{{n+1}}-e_{u}^{n+1}\right)\right\|\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|e_{w}^{n+1}\right\|\\ &+2\tau \lambda\left\|\frac{1}{\sigma^{n+1}}\right\|_{L^{\infty}}\left\|\sigma^{n+1} e_{u}^{n+1}\right\|\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|e_w^{n+1}\right\|\\ &+2 \tau \lambda\left\|\int_{t_{n}}^{t_{n+1}} \nabla \phi_{t}(t) d t\right\|_{L^{2}}\left\|\tilde{\boldsymbol{u}}^{n+1}\right\|_{L^{3}}\left\|e_{w}^{n+1}\right\|_{L^{6}}\\ &+2 \tau \lambda\left\|\nabla e_{\phi}^{n}\right\|\left\|\tilde{\boldsymbol{u}}^{n+1}\right\| _{L^{3}}\left\|e_{w}^{n+1}\right\|_{L^{6}}\\ \leq& C \tau^{3}+\frac{\tau\lambda}{2}\left\|e_w^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}+\frac{\tau\lambda\gamma}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2, \\ \end{split} \end{equation} (3.50)

    and

    \begin{equation} \begin{split} M_8 & = -2 \tau\left(u\left(t_{n+1}\right) \nabla \phi(t_{n+1} )-\tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right) \\ & \leq C\tau^{3}+\frac{1}{4}\left\|e_{\phi}^{n+1}\right\|^{2}+C\tau\left\|\sigma^{n+1} e_u^{n+1}\right\|^{2}+\frac{\lambda}{2}\left\|\nabla e_{\phi}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}.\\ \end{split} \end{equation} (3.51)

    Next, we estimate M_{10} + M_{11} and M_{12} + M_{13} as follows:

    \begin{equation} \begin{split} M_{10}+M_{11} \leq C\tau^{3}+C \tau\left(\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}\right)+\frac{\tau\gamma}{2} \left\|e_{w}^{n+1}\right\|^{2}. \end{split} \end{equation} (3.52)
    \begin{align} &\left(\tilde{G}^{n+1}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ = &\tau\left(\tilde{G}^{n+1}, \gamma \triangle e_w^{n+1}-u(t_{n+1} ) \nabla \phi(t _{n+1})+\tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}+R_{\phi}^{n+1}\right)\\ = &\tau\left(\tilde{G}^{n+1}, \gamma \triangle e_w^{n+1}-\left(\tilde{e}_{u}^{n+1}-e_{u}^{n+1}\right) \nabla \phi\left(t_{n+1}\right)-e_{u}^{n+1} \nabla \phi(t _{n+1})-\tilde{\boldsymbol{u}}^{n+1} \int_{t _n}^{t _{n+1}} \nabla \phi_t(t) d t\right)\\ &-\tau\left(\tilde{G}^{n+1}, \tilde{\boldsymbol{u}} \nabla e_{\phi}^{n}-R_{\phi}^{n+1}\right) \\ \leq &C \tau^{3}+C\tau\left\|\tilde{G}^{n+1}\right\|^{2}+\frac{\tau\gamma \lambda}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+C\tau\left\|\nabla \tilde{G}^{n+1}\right\|^{2}+\frac{1}{4}\left\|\sigma^{n+1} e_u^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2\\ \leq& C\tau^{3}+\frac{\tau \gamma \lambda}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+\frac{1}{4}\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2\\ &+C\tau\left(\left\|e_{w}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{\phi}^{n}\right\|^{2}\right), \end{align} (3.53)

    and

    \begin{align} \left(\int_{t_{n}}^{t_{n+1}} \phi_{t}(t) d t, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\leq C \tau^{3}+\frac{\tau \gamma \lambda}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+\frac{1}{4}\left\|\sigma^{n+1} e^{n+1}\right\|^{2}+C \tau\|\nabla e _{\phi}^{n}\|^{2}. \end{align}

    From (3.29), combining the two inequalities above, we deduce that

    \begin{equation} \begin{split} M_{12}+M_{13} \leq&-\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|(e_{\phi}^{n+1})^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e _{\phi}^{n}\right)\right\|^{2}\right)-\frac{2 \lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e _{\phi}^{n}\right\|^{2}\\ &+C \tau^{3}+\frac{\tau\gamma \lambda}{2}\left\|\nabla e_{w}^{n+1}\right\|^{2}+\frac{1}{2}\left\|\sigma^{n+1} e_u^{n+1}\right\|^{2} \\ &+C\tau\left(\left\|e_{w}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{\phi}^{n}\right\|^{2}\right).\\ \end{split} \end{equation} (3.54)

    Plugging the above inequality into (3.41) gives

    \begin{align} &\frac{1}{2}\left(\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|\sigma^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}\right)+ \eta\tau\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \tau \gamma \lambda\left\|\nabla e_{w}^{n+1}\right\|^{2} \\ &+\frac{1}{2}\left(\left\|e_{\phi}^{n+1}\right\|^{2}-\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2}\right)+ \tau \gamma\left\|e_{w}^{n+1}\right\|^{2}+\frac{\lambda}{2}\left\|\nabla e_{\phi}^{n+1}\right\|^{2}-\lambda\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+\lambda\left\|\nabla e _{\phi}^{n+1}-\nabla e_{\phi}^{n} \right\|^{2}+\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|(e_{\phi}^{n+1})^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e _{\phi}^{n}\right)\right\|^{2}\right)+\frac{2 \lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e _{\phi}^{n}\right\|^{2}\\ \leq& C \tau^{3}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2\\ &+C\tau\left(\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{w}^{n+1}\right\|^{2} +\left\|e_{w}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}\right){;} \end{align} (3.55)

    for sufficiently small \tau , taking the sum of (3.55) from 0 to N -1 and using the discrete Gronwall inequality, we can obtain Lemma 6.

    Theorem 4. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) are valid. For sufficiently small \tau , there are the following error estimates:

    \begin{equation} \begin{aligned} &\left\|\sigma\left(t_{n}\right)\boldsymbol{u}\left(t_{n}\right)-\sigma^{n}\boldsymbol{u}^{n}\right\|^{2}+\left\|\rho\left(t_{n}\right)-\rho^{n}\right\|^{2}+\left\|\nabla\left(w\left(t_{n}\right)-w^{n}\right)\right\|^{2} \\ &+\left\|\nabla\left(\phi\left(t_{n}\right)-\phi^{n}\right)\right\|^{2}+\eta \tau \sum\limits_{m = 1}^{n}\left\|\nabla\left(\boldsymbol{u}\left(t_{m}\right)-\boldsymbol{u}^{m}\right)\right\|^{2} \leq C \tau^{2} \end{aligned} \end{equation} (3.56)

    for all 1 \leq n \leq N .

    Proof. From Lemmas 5 and 6, we obtain

    \begin{equation} \left\|\rho\left(t_{n}\right)-\rho^{n}\right\|^{2}+\left\|\nabla\left(w\left(t_{n}\right)-w^{n}\right)\right\|^{2}+\left\|\nabla\left(\phi\left(t_{n}\right)-\phi^{n}\right)\right\|^{2} +\eta \tau \sum\limits_{m = 1}^{n}\left\|\nabla\left(\boldsymbol{u}\left(t_{m}\right)-\boldsymbol{u}^{m}\right)\right\|^{2} \leq C \tau^{2}. \end{equation} (3.57)

    Thus, we only prove that

    \begin{equation} \left\|\sigma\left(t_{n}\right)\boldsymbol{u}\left(t_{n}\right)-\sigma^{n}\boldsymbol{u}^{n}\right\|^{2} \leq C \tau^{2} . \end{equation} (3.58)

    In fact, we have

    \begin{equation} \begin{aligned} \left\|\sigma\left(t_{n}\right)\boldsymbol{u}\left(t_{n}\right)-\sigma^{n}\boldsymbol{u}^{n}\right\|^{2} & \leq C\left\|\left(\sigma\left(t_{n}\right)-\sigma^{n}\right)\boldsymbol{u}\left(t_{n}\right)\right\|^{2}+\left\|\sigma^{n} e_{u}^{n}\right\|^{2} \\ & \leq C\left(\left\|e_{\rho}^{n}\right\|^{2}+\left\|\sigma^{n} e_{u}^{n}\right\|^{2}\right) \leq C \tau^{2}. \end{aligned} \end{equation} (3.59)

    Theorem 3.6 states that both \sigma^{n}\boldsymbol{u}^{n}, \; \rho^{n}, \; \boldsymbol{u}^{n} , \phi^{n} and w^{n} are order 1 approximations to \sigma\boldsymbol{u}, \; \rho, \; \boldsymbol{u} , \phi and w in l^{\infty}\left(L^{2}(\Omega)\right) , l^{\infty}\left(L^{2}(\Omega)\right) , l^{2}\left(H_{0}^{1}(\Omega)\right) , l^{\infty}\left(H_0^{1}(\Omega)\right) and l^{\infty}\left(H_0^{1}(\Omega)\right) , respectively. Finally, we can obtain order \frac{1}{2} error estimates for p approximation in l^{\infty}\left(L^{2}(\Omega)\right) .

    Theorem 5. Under the assumptions in Theorem 4, the following holds true:

    \begin{equation} \tau\sum\limits_{m = 1}^{N}\left\|p\left(t_{m}\right)-p^{m}\right\| \leq C \tau. \end{equation} (3.60)

    Proof. Let us rewrite (3.40) as

    \begin{equation} \begin{split} -\nabla e_{p}^{n+1} = &\rho^{n} \frac{e_{u}^{n+1}-e_{u}^{n}}{\tau}+\frac{\rho^{n+1}-\rho^{n}}{\tau}\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)-\eta \triangle e_{u}^{n+1}-R_{u}^{n+1}\\ &+\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1}\\ &- \lambda w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)+ \lambda w^{n+1} \nabla \phi^{n}.\\ \end{split} \end{equation} (3.61)

    To prove the theorem, we introduce the discrete inf-sup condition, i.e.,

    \begin{equation} \beta\left\|e_{p}^{n+1}\right\| \leq \frac{\left(\nabla e_{p}^{n+1}, v\right)}{\|\nabla v\|}. \end{equation} (3.62)

    Then, we can constrain the products of the right-hand side of (3.61) with an arbitrary v \in V as follows:

    \begin{equation*} \begin{aligned} \frac{1}{\tau}\left(\rho^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right), v\right) &\leq C \tau^{-1}\left\|\sigma^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right)\right\| \| \nabla v \|, \\ -\eta\left(\Delta e_{u}^{n+1}, v\right)-\left(R_{u}^{n+1}, v\right) &\leq C\left(\tau^{2}+\left\|\nabla e_{u}^{n+1}\right\|\right) \| \nabla v \|, \end{aligned} \end{equation*}

    and

    \begin{equation} \begin{split} &- \lambda\left(w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-w^{n+1} \nabla \phi^{n}, v\right)\\ = & -\lambda\left(w\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right)+w\left(t_{n+1}\right) \nabla e_{\phi}^{n}+\nabla \phi^{n} e_{w}^{n+1}, v\right)\\ \leq& \lambda \left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\| \int_{t_{n}}^{t_{n+1}} \nabla \phi_{t}(t) d t\right\|\| v\|+ \lambda\left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\|\left\|v\right\|\\ &+ \lambda\left\|\nabla \phi^{n}\right\|_{L^2}\left\|e_w^{n+1}\right\|_{L^{3}}\left\|v\right\|_{L^{6}}\\ \leq& C\left( \tau+\left\|\nabla e_{\phi}^{n}\right\|+\left\|\nabla e_{w}^{n+1}\right\|\right)\left\| \nabla v \right\|.\\ \end{split} \end{equation} (3.63)

    For the second term on the right-hand side, using \|v\|_{L^{3}} \leq\|v\|^{\frac{1}{2}}\|\nabla v\|^{\frac{1}{2}} , we have

    \begin{equation} \begin{aligned} &\left(\frac{\rho^{n+1}-\rho^{n}}{\tau}\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right), v\right) \\ &\leq C \tau^{-1}\left\|\rho^{n+1}-\rho^{n}\right\|\left\|e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right\|_{L^{3}}\|v\|_{L^{6}} \\ &\leq C \tau^{-1}\left\|\rho^{n+1}-\rho^{n}\right\|\left\|e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right\|^{\frac{1}{2}}\left\|\nabla\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)\right\|^{\frac{1}{2}}\|\nabla v\| \\ &\leq C\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)\right\|^{\frac{1}{2}}\left\|\nabla\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)\right\|^{\frac{1}{2}}\|\nabla v\|. \end{aligned} \end{equation} (3.64)

    By the split method, we arrive at

    \begin{equation} \begin{aligned} &\left(\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \widetilde{\boldsymbol{u}}^{n+1}, v\right) \\ & = \left(\rho^{n+1}\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), v\right)+\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{e}_{u}^{n+1}, v\right) \\ &\quad+\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), v\right)+\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n+1}\right) \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), v\right). \end{aligned} \end{equation} (3.65)

    Utilizing H \ddot{o} lder's inequality and Young's inequality, we derive

    \begin{equation} \begin{split} &\left(\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \widetilde{\boldsymbol{u}}^{n+1}, v\right) \\ \leq& C\left(\tau+\left\| e_{\rho}^{n+1}\right\|+\left\|\nabla e_{u}^{n}\right\|+\left\|\nabla \tilde{e}_{u}^{n+1}\right\|\right) \| \nabla v\|.\\ \end{split} \end{equation} (3.66)

    By adding up the inequalities and incorporating (3.62), we get the desired result.

    In this section, we will present some numerical experiments to prove the validity and accuracy of our method. In the following simulation, for phase field \phi , chemical potential w and pressure p , we take the P1 finite element space (continuous piecewise linear), and for fluid velocity \boldsymbol{u} , we take the P2 finite element space. All experiments were conducted in Freefem++. We fixed \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1, T = 0.1, \rho_1 = 1 and \rho_2 = 3. The computational domain and the initial conditions were taken as

    \begin{equation*} \begin{split} &\Omega = \{(x, y)\in R^2:x^2+y^2 < 1\}, \\ &\phi_0 = cos(\pi x)cos(\pi y), \\ &\boldsymbol{u}_0 = ( \pi cos(\pi x)sin(\pi y), -\pi sin(\pi x)cos(\pi y)).\\ \end{split} \end{equation*}

    Tables 1 and 2 verify that (3.4) is the first-order convergence rate O(\tau) of (\phi, \sigma\boldsymbol{u}, \boldsymbol{u}, \rho) , which is consistent with the conclusion obtained from theoretical analysis. It is only proved in Theorem 3.2 that the pressure is the half-order convergence rate O(\tau^{\frac{1}{2}}) because of technical reasons. However, the numerical results on p in Table 1 still reach the first-order optimal convergence rate O(\tau) . Tables 3 and 4 show the convergence rate with another set of parameters.

    Table 1.  The order of temporal convergence with \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1 .
    \tau \Vert \phi\Vert_{H^1} Rate \Vert p\Vert_{L^2} Rate
    0.007812 0.102393 0.469438
    0.003906 0.0653769 0.647262 0.282333 0.733534
    0.001953 0.0367805 0.829842 0.143976 0.971574
    0.000976 0.0196276 0.906053 0.0725853 0.988076

     | Show Table
    DownLoad: CSV
    Table 2.  The order of temporal convergence with \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1 .
    \tau \Vert \sigma \boldsymbol{u}\Vert_{L^2} Rate \Vert \boldsymbol{u}\Vert_{H^1} Rate \Vert \rho\Vert_{L^2} Rate
    0.007812 0.114448 0.554732 0.0779992
    0.003906 0.0610184 0.90737 0.2956 0.90814 0.0463154 0.751968
    0.001953 0.02533 1.2684 0.132451 1.15819 0.022386 1.04889
    0.000976 0.0127883 0.98602 0.071371 0.89204 0.0107557 1.05749

     | Show Table
    DownLoad: CSV
    Table 3.  The order of temporal convergence with \eta = 0.1, \lambda = 0.2, \gamma = 0.0003, \varepsilon = 0.08 .
    \tau \Vert \phi\Vert_{H^1} Rate \Vert p\Vert_{L^2} Rate
    0.007812 0.041787 0.0794722
    0.003906 0.0215011 0.958642 0.0435045 0.869285
    0.001953 0.0106192 1.01774 0.0210772 1.04548
    0.000976 0.0054229 0.96953 0.0106356 0.986787

     | Show Table
    DownLoad: CSV
    Table 4.  The order of temporal convergence with \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1 .
    \tau \Vert \sigma \boldsymbol{u}\Vert_{L^2} Rate \Vert \boldsymbol{u}\Vert_{H^1} Rate \Vert \rho\Vert_{L^2} Rate
    0.007812 0.0308723 0.387768 0.0585448
    0.003906 0.0161585 0.934018 0.201806 0.942225 0.0300705 0.961192
    0.001953 0.00587585 1.45943 0.0765909 1.39772 0.00932566 1.68907
    0.000976 0.00302534 0.957701 0.0398452 0.942768 0.00468607 0.992826

     | Show Table
    DownLoad: CSV

    Figure 1 shows the evolution of the total energy at \tau = 0.02 . The downward trend of the energy curve confirms that our scheme is unconditionally energy stable. We also see a downward trend in energy when using different parameters. The energy curves for different time steps are shown in Figure 2 as a result of keeping the other parameters unchanged. It can been found that the curves are very similar, which means that the scheme is robust against different time steps.

    Figure 1.  Energy evolution with different parameters.
    Figure 2.  Energy evolution with different time step sizes.

    To solve the Cahn-Hilliard phase field model for two-phase incompressible flows with variable density, we have designed a novel time marching scheme, which can significantly improve the calculation efficiency. The method is efficient because we decoupled the pressure from the velocity and phase field. We have also proved the unconditional energy stability, presented the error analysis and provided various numerical examples to demonstrate the stability and accuracy of the scheme. In addition, the decoupling method developed in this paper is universally applicable, and this method is always applicable for the generation of an effective fully decoupling scheme.

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

    This work was supported by the Research Project Supported by Shanxi Scholarship Council of China (No. 2021-029), Shanxi Provincial International Cooperation Base and Platform Project (202104041101019), Shanxi Province Natural Science Research (202203021211129), Shanxi Province Natural Science Research (No. 202203021212249) and Special/Youth Foundation of Taiyuan University of Technology (No. 2022QN101).

    All authors declare that they have no competing interests in this paper.



    [1] L. A. Zadeh, Fuzzy sets, in: Fuzzy sets, fuzzy logic, and fuzzy systems: Selected papers by Lotfi A Zadeh, 1996,394–432. https://doi.org/10.1142/9789814261302 0021
    [2] K. T. Atanassov, Intuitionistic fuzzy sets, Fuzzy Set. Syst., 20 (1986), 87–96. https://doi.org/10.1016/S0165-0114(86)80034-3 doi: 10.1016/S0165-0114(86)80034-3
    [3] R. R. Yager, Pythagorean membership grades in multi-criteria decision making, IEEE T. Fuzzy Syst., 2014,958–965. https://doi.org/10.1109/TFUZZ.2013.2278989 doi: 10.1109/TFUZZ.2013.2278989
    [4] R. R. Yager, Generalized orthopair fuzzy sets, IEEE T. Fuzzy Syst., 25 (2017), 1222–1230. https://doi.org/10.1109/TFUZZ.2016.2604005 doi: 10.1109/TFUZZ.2016.2604005
    [5] M. Riaz, M. R. Hashmi, Linear Diophantine fuzzy set and its applications towards multi-attribute decision-making problems, J. Intell. Fuzzy Syst., 37 (2019), 5417–5439. https://doi.org/10.3233/JIFS-190550 doi: 10.3233/JIFS-190550
    [6] A. O. Almagrabi, S. Abdullah, M. Shams, Y. D. Al-Otaibi, S. Ashraf, A new approach to q-linear Diophantine fuzzy emergency decision support system for COVID19, J. Amb. Intel. Hum. Comp., 13 (2022), 1687–1713. https://doi.org/10.1007/s12652-021-03130-y doi: 10.1007/s12652-021-03130-y
    [7] S. Ashraf, S. Abdullah, Spherical aggregation operators and their application in multi-attribute group decision‐making, Int. J. Intell. Syst., 34 (2019), 493–523. https://doi.org/10.1002/int.22062 doi: 10.1002/int.22062
    [8] S. Ashraf, S. Abdullah, T. Mahmood, F. Ghani, T. Mahmood, Spherical fuzzy sets and their applications in multi-attribute decision making problems, J. Intell. Fuzzy Syst., 36 (2019), 2829–2844. https://doi.org/10.3233/JIFS-172009 doi: 10.3233/JIFS-172009
    [9] M. Riaz, M. R. Hashmi, D. Pamucar, Y. M. Chu, Spherical linear Diophantine fuzzy sets with modeling uncertainties in MCDM, Comput. Model. Eng. Sci., 126 (2021), 1125–1164. https://doi.org/10.32604/cmes.2021.013699 doi: 10.32604/cmes.2021.013699
    [10] J. Dombi, A general class of fuzzy operators, the demorgan class of fuzzy operators and fuzziness measures induced by fuzzy operators, Fuzzy Set. Syst., 8 (1982), 149–163. https://doi.org/10.1016/0165-0114(82)90005-7 doi: 10.1016/0165-0114(82)90005-7
    [11] C. Jana, G. Muhiuddin, M. Pal, D. Al-Kadi, Intuitionistic fuzzy dombi hybrid decision-making method and their applications to enterprise financial performance evaluation, Math. Probl. Eng., 2021.
    [12] C. Jane, M. Pal, G. Wei, Multiple attribute decision making method based on intuitionistic Dombi operators and its application in mutual fund evaluation, Arch. Control Sci., 2020,437–470.
    [13] Z. Li, H. Gao, G. Wei, Methods for multiple attribute group decision making based on intuitionistic fuzzy dombi hamy mean operators, Symmetry, 10 (2018), 574. https://doi.org/10.3390/sym10110574 doi: 10.3390/sym10110574
    [14] T. Senapati, V. Simic, A. Saha, M. Dobrodolac, Y. Rong, E. B. Tirkolaee, Intuitionistic fuzzy power Aczel-Alsina model for prioritization of sustainable transportation sharing practices, Eng. Appl. Artif. Intell., 119 (2023), 105716. https://doi.org/10.1016/j.engappai.2022.105716 doi: 10.1016/j.engappai.2022.105716
    [15] C. Jana, T. Senapati, M. Pal, Pythagorean fuzzy Dombi aggregation operators and its applications in multiple attribute decision-making, Int. J. Intell. Syst., 34 (2019), 2019–2038. https://doi.org/10.1002/int.22125 doi: 10.1002/int.22125
    [16] M. Akram, W. A. Dudek, J. M. Dar, Pythagorean Dombi fuzzy aggregation operators with application in multicriteria decision-making, Int. J. Intell. Syst., 34 (2019), 3000–3019. https://doi.org/10.1002/int.22183 doi: 10.1002/int.22183
    [17] P. Liu, Q. Khan, T. Mahmood, R. A. Khan, H. U. Khan, Some improved pythagorean fuzzy Dombi power aggregation operators with application in multiple-attribute decision making, J Intel. Fuzzy Syst., 40 (2021), 9237–9257. https://doi.org/10.3233/JIFS-201723 doi: 10.3233/JIFS-201723
    [18] M. Akram, A. Khan, A. Luqman, T. Senapati, D. Pamucar, An extended MARCOS method for MCGDM under 2-tuple linguistic q-rung picture fuzzy environment, Eng. Appl. Artif. Intel., 120 (2023), 105892.
    [19] T. Senapati, Approaches to multi-attribute decision-making based on picture fuzzy Aczel-Alsina average aggregation operators, Comput. Appl. Math., 41 (2022), 40.
    [20] C. Jana, G. Muhiuddin, M. Pal, Some Dombi aggregation of Q-rung orthopair fuzzy numbers in multiple-attribute decision making, Int. J. Intell. Syst., 34 (2019), 3220–3240. https://doi.org/10.1002/int.22191 doi: 10.1002/int.22191
    [21] S. Ashraf, S. Abdullah, T. Mahmood, Spherical fuzzy Dombi aggregation operators and their application in group decision making problems, J. Amb. Intel. Hum. Comp., 11 (2020), 2731–2749. https://doi.org/10.1007/s12652-019-01333-y doi: 10.1007/s12652-019-01333-y
    [22] J. Aldring, S. Santhoshkumar, D. Ajay, A decision making approach using linear diophantine fuzzy sets with Dombi operations, In: International Conference on Intelligent and Fuzzy Systems, 2022,684–692.
    [23] S. Ashraf, S. Abdullah, M. Aslam, M. Qiyas, M. A. Kutbi, Spherical fuzzy sets and its representation of spherical fuzzy t-norms and t-conorms, J. Intell. Fuzzy Syst., 36 (2019), 6089–6102. https://doi.org/10.3233/JIFS-181941 doi: 10.3233/JIFS-181941
    [24] P. Liu, B. Zhu, P. Wang, M. Shen, An approach based on linguistic spherical fuzzy sets for public evaluation of shared bicycles in China, Eng. Appl. Artif. Intell., 87 (2020), 103295. https://doi.org/10.1016/j.engappai.2019.103295 doi: 10.1016/j.engappai.2019.103295
    [25] M. Rafiq, S. Ashraf, S. Abdullah, T. Mahmood, S. Muhammad, The cosine similarity measures of spherical fuzzy sets and their applications in decision making, J. Intell. Fuzzy Syst., 36 (2019), 6059–6073. https://doi.org/10.3233/JIFS-181922 doi: 10.3233/JIFS-181922
    [26] S. Ashraf, S. Abdullah, T. Mahmood, Spherical fuzzy Dombi aggregation operators and their application in group decision making problems, J. Amb. Intel. Hum. Comp., 11 (2020), 2731–2749. https://doi.org/10.1007/s12652-019-01333-y doi: 10.1007/s12652-019-01333-y
    [27] Q. Khan, T. Mahmood, K. Ullah, Applications of improved spherical fuzzy Dombi aggregation operators in decision support system, Soft Comput., 25 (2021), 9097–9119. https://doi.org/10.1007/s00500-021-05829-8 doi: 10.1007/s00500-021-05829-8
    [28] G. Elif, H. Aygün, Generalized spherical fuzzy Einstein aggregation operators: Application to multi-criteria group decision-making problems, In: Conference Proceedings of Science and Technology, 3 (2020), 227–235.
    [29] S. Ashraf, S. Abdullah, Decision aid modeling based on sine trigonometric spherical fuzzy aggregation information, Soft Comput., 25 (2021), 8549–8572. https://doi.org/10.1007/s00500-021-05712-6 doi: 10.1007/s00500-021-05712-6
    [30] M. Qiyas, S. Abdullah, S. Khan, M. Naeem, Multi-attribute group decision making based on sine trigonometric spherical fuzzy aggregation operators, Granular Comput., 7 (2022), 141–162. https://doi.org/10.1007/s41066-021-00256-4 doi: 10.1007/s41066-021-00256-4
    [31] R. Chinram, S. Ashraf, S. Abdullah, P. Petchkaew, Decision support technique based on spherical fuzzy Yager aggregation operators and their application in wind power plant locations: A case study of Jhimpir, Pakistan, J. Math., 2020. https://doi.org/10.1155/2020/8824032 doi: 10.1155/2020/8824032
    [32] M. Riaz, H. M. Athar Farid, D. Pamucar, S. Tanveer, Spherical fuzzy information aggregation based on Aczel-Alsina operations and data analysis for supply chain, Math. Probl. Eng., 2022.
    [33] M. Qiyas, S. Abdullah, M. Naeem, Spherical uncertain linguistic Hamacher aggregation operators and their application on achieving consistent opinion fusion in group decision making, Int. J. Intell. Comput., 2021.
    [34] G. Wei, H. Zhang, X. Chen, Spherical fuzzy Hamacher power aggregation operators based on entropy for multiple attribute group decision making.
    [35] F. Kutlu-Gündoğdu, C. Kahraman, A novel VIKOR method using spherical fuzzy sets and its application to warehouse site selection, J. Intell. Fuzzy Syst., 37 (2019), 1197–1211. https://doi.org/10.3233/JIFS-182651 doi: 10.3233/JIFS-182651
    [36] F. Kutlu-Gündoğdu, C. Kahraman, A. Karaşan, Spherical fuzzy VIKOR method and its application to waste management, In: International Conference on Intelligent and Fuzzy Systems, Springer Cham., 2019,997–1005.
    [37] M. Akram, C. Kahraman, K. Zahid, Group decision-making based on complex spherical fuzzy VIKOR approach, Knowl.-Based Syst., 216 (2021), 106793. https://doi.org/10.1016/j.knosys.2021.106793 doi: 10.1016/j.knosys.2021.106793
    [38] M. Akram, M. Shabir, A. Adeel, A. N. Al-Kenani, A multiattribute decision-making framework: VIKOR method with complex spherical fuzzy-soft sets, Math. Probl. Eng., 2021.
    [39] I. M. Sharaf, Spherical fuzzy VIKOR with SWAM and SWGM operators for MCDM, In: Decision Making with Spherical Fuzzy Sets, Springer, Cham. 2021,217–240.
    [40] E. Ayyildiz, A. Taskin, A novel spherical fuzzy AHP-VIKOR methodology to determine serving petrol station selection during COVID-19 lockdown: A pilot study for İstanbul, Socio-Econ. Plan. Sci., 2022, 101345.
    [41] T. Y. Chen, An evolved VIKOR method for multiple-criteria compromise ranking modeling under T-spherical fuzzy uncertainty, Adv. Eng. Inform., 54 (2022), 101802. https://doi.org/10.1016/j.aei.2022.101802 doi: 10.1016/j.aei.2022.101802
    [42] M. Qiyas, S. Abdullah, S. Ashraf, M. Aslam, Utilizing linguistic picture fuzzy aggregation operators for multiple-attribute decision-making problems, Int J Fuzzy Syst., 22 (2020), 310–320. https://doi.org/10.1007/s40815-019-00726-7 doi: 10.1007/s40815-019-00726-7
    [43] V. T. Nguyen, R. Chaysiri, Spherical fuzzy AHP-VIKOR model application in solar energy location selection problem: A case study in Vietnam, In : {International Joint Symposium on Artificial Intelligence and Natural Language Processing} (iSAI-NLP), IEEE, 2022, 1–5.
    [44] I. M. Sharaf, A new approach for spherical fuzzy TOPSIS and spherical fuzzy VIKOR applied to the evaluation of hydrogen storage systems, Soft Comput., 27 (2023), 1–21. https://doi.org/10.1007/s00500-022-07749-7 doi: 10.1007/s00500-022-07749-7
    [45] I. Riali, M. Fareh, M. C. Ibnaissa, M. Bellil, A semantic-based approach for hepatitis C virus prediction and diagnosis using a fuzzy ontology and a fuzzy Bayesian network, J. Intell. Fuzzy Syst., 44 (2022), 1–15. https://doi.org/10.3233/JIFS-213563 doi: 10.3233/JIFS-213563
    [46] K. Naeem, M. Riaz, F. Karaaslan, A mathematical approach to medical diagnosis via Pythagorean fuzzy soft TOPSIS, VIKOR and generalized aggregation operators, Complex Intell. Syst., 7 (2021), 2783–2795. https://doi.org/10.1007/s40747-021-00458-y doi: 10.1007/s40747-021-00458-y
    [47] B. Khaoula, B. Imene, G. Wahiba, G. Abdelkader, L. Slimane, Intelligent analysis of some factors accompanying hepatitis B, Mol. Sci. Appl., 2 (2022), 61–71. https://doi.org/10.37394/232023.2022.2.7 doi: 10.37394/232023.2022.2.7
    [48] W. Liu, X. Liu, M. Peng, G. Q. Chen, P. H. Liu, X. W. Cui, et al., Artificial intelligence for hepatitis evaluation, World J. Gastroentero., 27 (2021), 5715.
    [49] M. Riaz, S. T. Tehrim, A robust extension of VIKOR method for bipolar fuzzy sets using connection numbers of SPA theory based metric spaces, Artif. Intell. Rev., 54 (2021), 561–591. https://doi.org/10.1007/s10462-020-09859-w doi: 10.1007/s10462-020-09859-w
    [50] A. Vaillant, Transaminase elevations during treatment of chronic hepatitis B infection: Safety considerations and role in achieving functional cure, Viruses, 13 (2021), 745. https://doi.org/10.3390/v13050745 doi: 10.3390/v13050745
    [51] J. Wei, X. Lin, The multiple attribute decision-making VIKOR method and its application, In: 2008 4th international conference on wireless communications, networking and mobile computing, IEEE, 2008, 1–4.
  • This article has been cited by:

    1. Jaeyong Choi, Seokjun Ham, Soobin Kwak, Youngjin Hwang, Junseok Kim, Stability analysis of an explicit numerical scheme for the Allen-Cahn equation with high-order polynomial potentials, 2024, 9, 2473-6988, 19332, 10.3934/math.2024941
    2. Soobin Kwak, Seokjun Ham, Jian Wang, Hyundong Kim, Junseok Kim, Effective perpendicular boundary conditions in phase-field models using Dirichlet boundary conditions, 2025, 0177-0667, 10.1007/s00366-025-02109-z
  • 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(1368) PDF downloads(48) Cited by(9)

Figures and Tables

Figures(10)  /  Tables(6)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog