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

On (p,q)-fractional linear Diophantine fuzzy sets and their applications via MADM approach

  • The integration of internationally sustainable practices into supply chain management methodologies is known as "green supply chain management". Reducing the supply chain's overall environmental impact is the main objective in order to improve corporate connections and the social, ecological, and economic ties with other nations. To accomplish appropriate and accurate measures to address the issue of emergency decision-making, the paper is divided into three major sections. First, the (p,q)-fractional linear Diophantine fuzzy set represents a new generalization of several fuzzy set theories, including the Pythagorean fuzzy set, q-rung orthopair fuzzy set, linear Diophantine fuzzy set, and q-rung linear Diophantine fuzzy set, with its key features thoroughly discussed. Additionally, aggregation operators are crucial for handling uncertainty in decision-making scenarios. Consequently, algebraic norms for (p,q)-fractional linear Diophantine fuzzy sets were established based on operational principles. In the second part of the study, we introduced a range of geometric aggregation operators and a series of averaging operators under the (p,q)-fractional linear Diophantine fuzzy set, all grounded in established operational rules. We also explained some flexible aspects for the invented operators. Furthermore, using the newly developed operators for (p,q)-fractional linear Diophantine fuzzy information, we constructed the multi-attribute decision-making (MADM) technique to assess the green supply chain management challenge. Last, we compared the ranking results of the produced approaches with the obtained ranking results of the techniques using several numerical instances to demonstrate the validity and superiority of the developed techniques. Finally, a few comparisons between the findings were made.

    Citation: Hanan Alohali, Muhammad Bilal Khan, Jorge E. Macías-Díaz, Fahad Sikander. On (p,q)-fractional linear Diophantine fuzzy sets and their applications via MADM approach[J]. AIMS Mathematics, 2024, 9(12): 35503-35532. doi: 10.3934/math.20241685

    Related Papers:

    [1] Hongying Jiao, Shuhai Zhu, Jinguo Zhang . Existence of infinitely many solutions for critical sub-elliptic systems via genus theory. Communications in Analysis and Mechanics, 2024, 16(2): 237-261. doi: 10.3934/cam.2024011
    [2] Jizheng Huang, Shuangshuang Ying . Hardy-Sobolev spaces of higher order associated to Hermite operator. Communications in Analysis and Mechanics, 2024, 16(4): 858-871. doi: 10.3934/cam.2024037
    [3] Xiulan Wu, Yaxin Zhao, Xiaoxin Yang . On a singular parabolic p-Laplacian equation with logarithmic nonlinearity. Communications in Analysis and Mechanics, 2024, 16(3): 528-553. doi: 10.3934/cam.2024025
    [4] Erlend Grong, Irina Markina . Harmonic maps into sub-Riemannian Lie groups. Communications in Analysis and Mechanics, 2023, 15(3): 515-532. doi: 10.3934/cam.2023025
    [5] Leandro Tavares . Solutions for a class of problems driven by an anisotropic (p,q)-Laplacian type operator. Communications in Analysis and Mechanics, 2023, 15(3): 533-550. doi: 10.3934/cam.2023026
    [6] Velimir Jurdjevic . Time optimal problems on Lie groups and applications to quantum control. Communications in Analysis and Mechanics, 2024, 16(2): 345-387. doi: 10.3934/cam.2024017
    [7] Zhiyong Wang, Kai Zhao, Pengtao Li, Yu Liu . Boundedness of square functions related with fractional Schrödinger semigroups on stratified Lie groups. Communications in Analysis and Mechanics, 2023, 15(3): 410-435. doi: 10.3934/cam.2023020
    [8] Ming Liu, Binhua Feng . Grand weighted variable Herz-Morrey spaces estimate for some operators. Communications in Analysis and Mechanics, 2025, 17(1): 290-316. doi: 10.3934/cam.2025012
    [9] Xiaotian Hao, Lingzhong Zeng . Eigenvalues of the bi-Xin-Laplacian on complete Riemannian manifolds. Communications in Analysis and Mechanics, 2023, 15(2): 162-176. doi: 10.3934/cam.2023009
    [10] Eleonora Amoroso, Angela Sciammetta, Patrick Winkert . Anisotropic (p,q)-Laplacian problems with superlinear nonlinearities. Communications in Analysis and Mechanics, 2024, 16(1): 1-23. doi: 10.3934/cam.2024001
  • The integration of internationally sustainable practices into supply chain management methodologies is known as "green supply chain management". Reducing the supply chain's overall environmental impact is the main objective in order to improve corporate connections and the social, ecological, and economic ties with other nations. To accomplish appropriate and accurate measures to address the issue of emergency decision-making, the paper is divided into three major sections. First, the (p,q)-fractional linear Diophantine fuzzy set represents a new generalization of several fuzzy set theories, including the Pythagorean fuzzy set, q-rung orthopair fuzzy set, linear Diophantine fuzzy set, and q-rung linear Diophantine fuzzy set, with its key features thoroughly discussed. Additionally, aggregation operators are crucial for handling uncertainty in decision-making scenarios. Consequently, algebraic norms for (p,q)-fractional linear Diophantine fuzzy sets were established based on operational principles. In the second part of the study, we introduced a range of geometric aggregation operators and a series of averaging operators under the (p,q)-fractional linear Diophantine fuzzy set, all grounded in established operational rules. We also explained some flexible aspects for the invented operators. Furthermore, using the newly developed operators for (p,q)-fractional linear Diophantine fuzzy information, we constructed the multi-attribute decision-making (MADM) technique to assess the green supply chain management challenge. Last, we compared the ranking results of the produced approaches with the obtained ranking results of the techniques using several numerical instances to demonstrate the validity and superiority of the developed techniques. Finally, a few comparisons between the findings were made.



    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, Inf. Control, 8 (1965), 338–353. https://doi.org/10.1016/S00199958(65)90241-X doi: 10.1016/S00199958(65)90241-X
    [2] G. Q. Huang, L. M. Xiao, W. Pedrycz, G. B. Zhang, L. Martinez, Failure mode and effect analysis using T-spherical fuzzy maximizing deviation and combined comparison solution methods, IEEE T. Reliab., 2022, 1–22. https://doi.org/10.1109/TR.2022.3194057
    [3] A. Sotoudeh-Anvari, A critical review on theoretical drawbacks and mathematical incorrect assumptions in fuzzy OR methods: Review from 2010 to 2020, Appl. Soft Comput., 93 (2020), 106354. https://doi.org/10.1016/j.asoc.2020.106354 doi: 10.1016/j.asoc.2020.106354
    [4] L. M. Xiao, G. Q. Huang, G. B. Zhang, An integrated risk assessment method using Z-fuzzy clouds and generalized TODIM, Qual. Reliab. Eng., 38 (2022), 1909–1943. https://doi.org/10.1002/qre.3062 doi: 10.1002/qre.3062
    [5] 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
    [6] K. T. Atanassov, G. Gargov, Interval valued intuitionistic fuzzy sets, Fuzzy Set. Syst., 31 (1989), 343–349. https://doi.org/10.1016/0165-0114(89)90205-4 doi: 10.1016/0165-0114(89)90205-4
    [7] Z. S. Xu, R. R. Yager, Some geometric aggregation operators based on intuitionistic fuzzy sets, Int. J. Gen. Syst., 35 (2006), 417–433. https://doi.org/10.1080/03081070600574353 doi: 10.1080/03081070600574353
    [8] M. B. Khan, A. M. Deaconu, J. Tayyebi, D. E. Spridon, Diamond intuitionistic fuzzy sets and their applications, IEEE Access, 2024. https://doi.org/10.1109/ACCESS.2024.3502202
    [9] H. Garg, Some series of intuitionistic fuzzy interactive averaging aggregation operators, Springer Plus, 5 (2016), 999. https://doi.org/10.1186/s40064-016-2591-9 doi: 10.1186/s40064-016-2591-9
    [10] Y. G. Xue, Y. Deng, Decision making under measure-based granular uncertainty with intuitionistic fuzzy sets, Appl. Intell., 51 (2021), 6224–6233. https://doi.org/10.1007/s10489-021-02216-6 doi: 10.1007/s10489-021-02216-6
    [11] S. M. Khalil, M. A. H. Hasab, Decision making using new distances of intuitionistic fuzzy sets and study their application in the universities, In: Intelligent and fuzzy techniques: Smart and innovative solutions, Cham: Springer, 2020. https://doi.org/10.1007/978-3-030-51156-246
    [12] R. R. Yager, Pythagorean fuzzy subsets, In: 2013 joint IFSA world congress and NAFIPS annual meeting (IFSA/NAFIPS), 2013, 57–61. https://doi.org/10.1109/IFSA-NAFIPS.2013.6608375
    [13] Y. Zhang, X. Wang, X. Wang, H. A. Mang, Virtual displacement based discontinuity layout optimization, Int. J. Numer. Meth. Eng., 123 (2022), 5682–5694. https://doi.org/10.1002/nme.7084 doi: 10.1002/nme.7084
    [14] Y. Zhang, Z. Gao, X. Wang, Q. Liu, Image representations of numerical simulations for training neural networks, Comput. Model. Eng. Sci., 134 (2023), 821–833. https://doi.org/10.32604/cmes.2022.022088 doi: 10.32604/cmes.2022.022088
    [15] B. Farhadinia, Similarity-based multi-criteria decision making technique of pythagorean fuzzy set, Artif. Intell. Rev., 55 (2022), 2103–2148. https://doi.org/10.1007/s10462-021-10054-8 doi: 10.1007/s10462-021-10054-8
    [16] R. R. Yager, Pythagorean membership grades in multi-criteria decision making, IEEE T. Fuzzy Syst., 22 (2014), 958–965. https://doi.org/10.1109/TFUZZ.2013.2278989 doi: 10.1109/TFUZZ.2013.2278989
    [17] H. Garg, New logarithmic operational laws and their aggregation operators for Pythagorean fuzzy set and their applications, Int. J. Intell. Syst., 34 (2019), 82–106. https://doi.org/10.1002/int.22043 doi: 10.1002/int.22043
    [18] M. J. Khan, P. Kumam, P. D. Liu, W. Kumam, S. Ashraf, A novel approach to generalized intuitionistic fuzzy soft sets and its application in decision support system, Mathematics, 7 (2019), 742. https://doi.org/10.3390/math7080742 doi: 10.3390/math7080742
    [19] H. Garg, A new generalized Pythagorean fuzzy information aggregation using Einstein operations and its application to decision making, Int. J. Intell. Syst., 31 (2016), 886–920. https://doi.org/10.1002/int.21809 doi: 10.1002/int.21809
    [20] H. Garg, Generalized Pythagorean fuzzy geometric aggregation operators using Einstein t-norm and t-conorm for multicriteria decision-making process, Int. J. Intell. Syst., 32 (2017), 597–630. https://doi.org/10.1002/int.21860 doi: 10.1002/int.21860
    [21] Z. M. Ma, Z. S. Xu, Symmetric Pythagorean fuzzy weighted geometric/averaging operators and their application in multicriteria decision-making problems, Int. J. Intell. Syst., 31 (2016), 1198-1219. https://doi.org/10.1002/int.21823 doi: 10.1002/int.21823
    [22] S. Z. Zeng, Pythagorean fuzzy multiattribute group decision making with probabilistic information and OWAapproach, Int. J. Intell. Syst., 32 (2017), 1136–1150. https://doi.org/10.1002/int.21886 doi: 10.1002/int.21886
    [23] D. Q. Li, W. Y. Zeng, Distance measure of Pythagorean fuzzy sets, Int. J. Intell. Syst., 33 (2018), 348–361. https://doi.org/10.1002/int.21934 doi: 10.1002/int.21934
    [24] M. A. Firozja, B. Agheli, E. B. Jamkhaneh, A new similarity measure for Pythagorean fuzzy sets, Complex Intell. Syst., 6 (2020), 67–74. https://doi.org/10.1007/s40747-019-0114-3 doi: 10.1007/s40747-019-0114-3
    [25] S. Zhang, Y. Hou, S. Zhang, M. Zhang, Fuzzy control model and simulation for nonlinear supply chain system with lead times, Complexity, 2017 (2017), 2017634. https://doi.org/10.1155/2017/2017634 doi: 10.1155/2017/2017634
    [26] S. Zhang, S. Li, S. Zhang, M. Zhang, Decision of lead‐time compression and stable operation of supply chain, Complexity, 2017 (2017), 7436764. https://doi.org/10.1155/2017/7436764 doi: 10.1155/2017/7436764
    [27] S. Zhang, C. Zhang, S. Zhang, M. Zhang, Discrete switched model and fuzzy robust control of dynamic supply chain network, Complexity, 2018 (2018), 3495096. https://doi.org/10.1155/2018/3495096 doi: 10.1155/2018/3495096
    [28] K. Y. Bai, X. M. Zhu, J. Wang, R. T. Zhang, Some partitioned Maclaurin symmetric mean based on q-rung orthopair fuzzy information for dealing with multi-attribute group decision making, Symmetry, 10 (2018), 383. https://doi.org/10.3390/sym10090383 doi: 10.3390/sym10090383
    [29] O. Barukab, S. Abdullah, S. Ashraf, M. Arif, S. A. Khan, A new approach to fuzzy TOPSIS method based on entropy measure under spherical fuzzy information, Entropy, 21 (2019), 1231. https://doi.org/10.3390/e21121231 doi: 10.3390/e21121231
    [30] E. Alsuwat, S. Alzahrani, H. Alsuwat, Detecting COVID-19 utilizing probabilistic graphical models, Int. J. Adv. Comput. Sci. Appl., 12 (2021), 786–793. https://doi.org/10.14569/IJACSA.2021.0120692 doi: 10.14569/IJACSA.2021.0120692
    [31] 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
    [32] Z. Zhao, H. Zhang, J. Shiau, W. Du, L. Ke, F. Wu, et al., Failure envelopes of rigid tripod pile foundation under combined vertical-horizontal-moment loadings in clay, Appl. Ocean Res., 150 (2024), 104131. https://doi.org/10.1016/j.apor.2024.104131 doi: 10.1016/j.apor.2024.104131
    [33] Y. Gao, Q. Liu, Y. Yang, K. Wang, Latent representation discretization for unsupervised text style generation, Inf. Process. Manag., 61 (2024), 103643. https://doi.org/10.1016/j.ipm.2024.103643 doi: 10.1016/j.ipm.2024.103643
    [34] Y. Luo, H. Zhang, Z. Chen, Q. Li, S. Ye, Q. Liu, Novel multidimensional composite development for aging resistance of SBS-modified asphalt by attaching zinc oxide on expanded vermiculite, Energ. Fuel., 38 (2024), 16772–16781. https://doi.org/10.1021/acs.energyfuels.4c02685 doi: 10.1021/acs.energyfuels.4c02685
    [35] A. Iampan, G. S. Garc, M. Riaz, H. M. A. Farid, R. Chinram, Linear Diophantine fuzzy Einstein aggregation operators for multi-criteria decision making problems, J. Math., 2021 (2021), 5548033. https://doi.org/10.1155/2021/5548033 doi: 10.1155/2021/5548033
    [36] S. Ayub, M. Shabir, M. Riaz, M. Aslam, R. Chinram, Linear Diophantine fuzzy relations and their algebraic properties with decision making, Symmetry, 13 (2021), 945. https://doi.org/10.3390/sym13060945 doi: 10.3390/sym13060945
    [37] T. Mahmood, Z. Ali, M. Aslam, R. Chinram, Generalized Hamacher aggregation operators based on linear Diophantine uncertain linguistic setting and their applications in decision-making problems, IEEE Access, 9 (2021), 126748–126764. https://doi.org/10.1109/ACCESS.2021.3110273 doi: 10.1109/ACCESS.2021.3110273
    [38] M. B. Khan, M. A. Noor, T. Abdeljawad, B. Abdalla, A. Althobaiti, Some fuzzy-interval integral inequalities for harmonically convex fuzzy-interval-valued functions, AIMS Math., 7 (2022), 349–370. https://doi.org/10.3934/math.2022024 doi: 10.3934/math.2022024
    [39] M. B. Khan, M. A. Noor, M. M. Al-Shomrani, L. Abdullah, Some novel inequalities for LR-h-convex interval-valued functions by means of pseudo order relation, Math. Meth. Appl. Sci., 45 (2022), 1310–1340. https://doi.org/10.1002/mma.7855 doi: 10.1002/mma.7855
    [40] M. B. Khan, M. A. Noor, K. I. Noor, K. S. Nisar, K. A. Ismail, A. Elfasakhany, Some inequalities for LR-(h1, h2)-convex interval-valued functions by means of pseudo order relation, Int. J. Comput. Intell. Syst., 14 (2021), 1–15. https://doi.org/10.1007/s44196-021-00032-x doi: 10.1007/s44196-021-00032-x
    [41] M. B. Khan, M. A. Noor, L. Abdullah, Y. M. Chu, Some new classes of preinvex fuzzy-interval-valued functions and inequalities, Int. J. Comput. Intell. Syst., 14 (2021), 1403–1418. https://doi.org/10.2991/ijcis.d.210409.001 doi: 10.2991/ijcis.d.210409.001
    [42] P. Liu, M. B. Khan, M. A. Noor, K. I. Noor, New Hermite-Hadamard and Jensen inequalities for log-s-convex fuzzy-interval-valued functions in the second sense, Complex Intell. Syst., 8 (2021), 1–15. https://doi.org/10.1007/s40747-021-00379-w doi: 10.1007/s40747-021-00379-w
    [43] Z. Y. Xiao, Y. J. Li, W. Zhang, Y. J. Han, D. Li, Q. Chen, et al., Enhancement of torque efficiency and spin Hall angle driven collaboratively by orbital torque and spin-orbit torque, Appl. Phys. Lett., 121 (2022), 075801. https://doi.org/10.1063/5.0086125 doi: 10.1063/5.0086125
    [44] Y. Li, J. Li, C. Feng, M. Wen, Y. Zhang, An interface constitutive model of plastic tensile-compressive damage under impact loading based on continuous-discontinuous framework, Comput. Geotech., 173 (2024), 106502. https://doi.org/10.1016/j.compgeo.2024.106502 doi: 10.1016/j.compgeo.2024.106502
    [45] C. Xu, M. Zhu, Q. Wang, J. Cui, Y. Huang, X. Huang, et al., TROP2-directed nanobody-drug conjugate elicited potent antitumor effect in pancreatic cancer, J. Nanobiotechnol., 21 (2023), 410. https://doi.org/10.1186/s12951-023-02183-9 doi: 10.1186/s12951-023-02183-9
    [46] 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. Intell. Human. Comp., 13 (2022), 1687–1713. https://doi.org/10.1007/s12652-021-03130-y doi: 10.1007/s12652-021-03130-y
    [47] 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
    [48] M. Qiyas, M. Naeem, S. Abdullah, N. Khan, A. Ali, Similarity measures based on q-rung linear Diophantine fuzzy sets and their application in logistics and supply chain management, J. Math., 2022 (2022), 4912964. https://doi.org/10.1155/2022/4912964 doi: 10.1155/2022/4912964
    [49] M. Gulistan, W. Pedrycz, Introduction to q-fractional fuzzy set, Int. J. Fuzzy Syst., 2024, 1–18. https://doi.org/10.1007/s40815-023-01633-8
    [50] T. Mahmood, M. Asif, U. ur Rehman, J. Ahmmad, T-bipolar soft semigroups and related results, Spectr. Mech. Eng. Oper. Res., 1 (2024), 258–271. https://doi.org/10.31181/smeor11202421 doi: 10.31181/smeor11202421
    [51] A. R. Mishra, P. Rani, F. Cavallaro, A. F. Alrasheedi, Assessment of sustainable wastewater treatment technologies using interval-valued intuitionistic fuzzy distance measure-based MAIRCA method, Facta Univ. -Ser. Mech., 21 (2023), 359–386. https://doi.org/10.22190/FUME230901034M doi: 10.22190/FUME230901034M
    [52] A. Hussain, K. Ullah, An intelligent decision support system for spherical fuzzy sugeno-weber aggregation operators and real-life applications, Spectr. Mech. Eng. Oper. Res., 1 (2024), 177–188. https://doi.org/10.31181/smeor11202415 doi: 10.31181/smeor11202415
    [53] J. Kannan, V. Jayakumar, M. Pethaperumal, Advanced fuzzy-based decision-making: the linear diophantine fuzzy CODAS method for logistic specialist selection, Spectr. Oper. Res., 2 (2025) 41–60. https://doi.org/10.31181/sor2120259 doi: 10.31181/sor2120259
    [54] P. Wang, B. Zhu, Y. Yu, Z. Ali, B. Almohsen, Complex intuitionistic fuzzy DOMBI prioritized aggregation operators and their application for resilient green supplier selection, Facta Univ. -Ser. Mech., 21 (2023) 339–357. https://doi.org/10.22190/FUME230805029W doi: 10.22190/FUME230805029W
    [55] H. M. A. Farid, M. Riaz, M. J. Khan, P. Kumam, K. Sitthithakerngkiet, Sustainable thermal power equipment supplier selection by Einstein prioritized linear Diophantine fuzzy aggregation operators, AIMS Math., 7 (2022), 11201–11242. https://doi.org/10.3934/math.2022627 doi: 10.3934/math.2022627
    [56] 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
    [57] L. M. Xiao, G. Q. Huang, W. Pedrycz, D. Pamucar, L. Martınez, G. B. Zhang, A q-rung orthopair fuzzy decision-making model with new score function and best-worst method for manufacturer selection, Inf. Sci., 608 (2022), 153–177. https://doi.org/10.1016/j.ins.2022.06.061 doi: 10.1016/j.ins.2022.06.061
    [58] S. Zhang, P. Zhang, M. Zhang, Fuzzy emergency model and robust emergency strategy of supply chain system under random supply disruptions, Complexity, 2019 (2019), 3092514. https://doi.org/10.1155/2019/3092514 doi: 10.1155/2019/3092514
  • This article has been cited by:

    1. Yuxia Guo, Yichen Hu, Infinitely many solutions for Hamiltonian system with critical growth, 2024, 13, 2191-950X, 10.1515/anona-2023-0134
    2. Xingyue He, Chenghua Gao, Jingjing Wang, k-convex solutions for multiparameter Dirichlet systems with k-Hessian operator and Lane-Emden type nonlinearities, 2024, 13, 2191-950X, 10.1515/anona-2023-0136
    3. Hongying Jiao, Shuhai Zhu, Jinguo Zhang, Existence of infinitely many solutions for critical sub-elliptic systems via genus theory, 2024, 16, 2836-3310, 237, 10.3934/cam.2024011
    4. Jinli Yang, Jiajing Miao, Algebraic Schouten solitons of Lorentzian Lie groups with Yano connections, 2023, 15, 2836-3310, 763, 10.3934/cam.2023037
  • Reader Comments
  • © 2024 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(806) PDF downloads(49) Cited by(0)

Figures and Tables

Figures(9)  /  Tables(12)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog