Keywords. FENE-P model, finite extensibility, polymer molecular diffusion, conformation tensor, trace barrier, entropy-compatible reconstruction, asymptotic-preserving schemes

MSC 2020. 65M12, 65M60, 76A10, 76M20

1 Introduction↩︎

The Oldroyd–B conformation tensor is constrained only by positive definiteness. In finitely extensible nonlinear elastic models the admissible set is smaller. If \(\boldsymbol{C}\) denotes the conformation tensor and \(L^2\) is the squared maximum extension parameter, a FENE-type closure requires \[\boldsymbol{C}\in \mathcal{D}_L:=\{\boldsymbol{C}\in\mathbb{S}_{++}^d:\;\operatorname{tr}\boldsymbol{C}<L^2\}.\] Equivalently, in index notation the model requires \(C_{kk}<L^2\). The coefficient multiplying the elastic stress becomes singular as \(\operatorname{tr}\boldsymbol{C}\uparrow L^2\). This changes the numerical problem qualitatively. A log-conformation update can preserve \(\boldsymbol{C}\succ0\) while still allowing \(\operatorname{tr}\boldsymbol{C}\) to approach or exceed \(L^2\) after reconstruction, interpolation, or nonlinear iteration. The finite-extensibility barrier must therefore be treated as part of the structure, not as an afterthought.

This paper develops a structure-preserving numerical analysis for this constraint. We focus on the FENE-P type conformation closure \[f_L(\boldsymbol{C})=\frac{L^2-d}{L^2-\operatorname{tr}\boldsymbol{C}},\qquad \boldsymbol{T}(\boldsymbol{C})=f_L(\boldsymbol{C})\boldsymbol{C}-\boldsymbol{I},\] where \(L^2>d\) and \(\boldsymbol{T}(\boldsymbol{I})=0\). The singular denominator is precisely the difficulty. The associated entropy variable is \[\boldsymbol{H}_L(\boldsymbol{C})=f_L(\boldsymbol{C})\boldsymbol{I}-\boldsymbol{C}^{-1}.\] It is the gradient of a convex barrier energy on \(\mathcal{D}_L\), and it yields the same polymeric-work cancellation as the Oldroyd–B entropy variable, provided the stress tensor used in the momentum equation is the same accepted tensor used in the entropy calculation.

The central idea is to treat \(\mathcal{D}_L\) as the computational state space. We combine five devices. First, a trace-barrier entropy prevents accepted states from reaching the finite-extensibility boundary. Second, a barrier logarithmic parametrization maps unconstrained symmetric matrices bijectively into \(\mathcal{D}_L\), so positivity and \(\operatorname{tr}\boldsymbol{C}<L^2\) can be enforced during nonlinear solves. Third, high-order reconstructed barrier-log states are accepted only if they satisfy a FENE entropy budget. Fourth, polymer molecular diffusion is discretized in a way that dissipates the same barrier entropy. Fifth, the momentum equation is assembled in the scaled FENE stress \[\boldsymbol{S}=\frac{\boldsymbol{T}(\boldsymbol{C})}{\mathrm{Wi}},\] which remains regular in the Newtonian relaxation limit. The resulting formulation is simultaneously barrier preserving, energy stable, and asymptotic preserving.

The contribution is not the introduction of the FENE-P model itself. Rather, it is a compatibility analysis for finite-extensibility constraints at the fully discrete level. The main results are:

1. a convex trace-barrier entropy whose gradient produces the FENE stress cancellation;

2. a barrier-log parametrization that enforces \(\boldsymbol{C}\succ0\) and \(\operatorname{tr}\boldsymbol{C}<L^2\) by construction;

3. an entropy-compatible barrier-log reconstruction selected by the largest admissible parameter on a logarithmic path;

4. a fully discrete free-energy inequality with relaxation and polymer molecular-diffusion dissipation;

5. an AP stress-closure estimate showing \(\boldsymbol{S}_h=\boldsymbol{D}(\boldsymbol{u}_h)+\mathcal{O}(\mathrm{Wi})\) at fixed \(h\) and \(\Delta t\);

6. a conditional relative-entropy estimate on compact subsets of the finite-extensibility domain;

7. numerical diagnostics that test the barrier, entropy correction, AP limit, coupled pressure–stress feedback, and high-Weissenberg behavior.

Related work on FENE closures, log-conformation variables, energy-stable viscoelastic schemes, and asymptotic-preserving relaxation discretizations is extensive; see, for example, [1]–[9]. The present paper is closest in spirit to structure-preserving conformation-tensor methods, but the finite-extensibility constraint introduces a trace barrier that is not present in Oldroyd–B. This extra constraint is the focus of the analysis.

The paper retains the entropy identity, finite-extensibility preservation, entropy-compatible barrier-log reconstruction, the discrete energy estimate, the AP closure, and the core benchmarks. This organization keeps the finite-extensibility mechanism visible while connecting each diagnostic to a structural claim.

Position relative to logarithmic and factor reconstructions↩︎

Classical log-conformation variables [2]–[4] and square-root or symmetric-factor variables [10] are designed primarily to preserve positive definiteness and to improve robustness in strongly stretched flows. For FENE-type models this is only the first admissibility condition. A positive tensor with \(\operatorname{tr}\boldsymbol{C}\ge L^2\) is outside the model domain, and a positive tensor with \(\operatorname{tr}\boldsymbol{C}<L^2\) may still be unacceptable if a reconstruction step injects a mesh-scale burst of FENE entropy. Table 2 summarizes the distinction.

Richter, Iaccarino, and Shaqfeh [11] used a different and important finite-extensibility safeguard in DNS of FENE-P flow past a cylinder: the trace equation is advanced first and the positive root of a scalar algebraic relation is selected so that the FENE denominator remains positive. That idea is a trace-first admissibility update. The present construction has a different purpose. It treats the whole tensor as the constrained state variable, enforces \(\boldsymbol{C}\succ0\) and \(\operatorname{tr}\boldsymbol{C}<L^2\) simultaneously by a barrier-log map, and then accepts high-order reconstructed states only through a FENE entropy budget. Thus the guarantee is not only denominator positivity in a component update, but compatibility with the fully discrete free-energy inequality.

3pt

Thus the proposed reconstruction has two layers. The barrier-log map enforces the geometry of the FENE state space, while the entropy acceptance rule enforces compatibility with the free-energy estimate. Both layers are needed for the fully discrete theorem below.

2 FENE-P model with polymer molecular diffusion↩︎

Let \(b=L^2>d\) and define \[f_b(\boldsymbol{C})=\frac{b-d}{b-\operatorname{tr}\boldsymbol{C}},\qquad \mathcal{D}_b=\{\boldsymbol{C}\in\mathbb{S}_{++}^d:\operatorname{tr}\boldsymbol{C}<b\}.\] The molecularly diffusive FENE-P type system is written as \[\begin{align} \mathrm{Re}(\partial_t\boldsymbol{u}+\boldsymbol{u}\cdot\nabla\boldsymbol{u})+\nabla p-\beta\Delta\boldsymbol{u} &=(1-\beta)\nabla\!\cdot\boldsymbol{S}+\boldsymbol{f}, \tag{1}\\ \nabla\!\cdot\boldsymbol{u}&=0, \tag{2}\\ \partial_t\boldsymbol{C}+\boldsymbol{u}\cdot\nabla\boldsymbol{C}-(\nabla\boldsymbol{u})\boldsymbol{C}-\boldsymbol{C}(\nabla\boldsymbol{u})^T &=-\frac{1}{\mathrm{Wi}}\boldsymbol{T}(\boldsymbol{C})+\varepsilon\Delta\boldsymbol{C}, \tag{3}\\ \boldsymbol{T}(\boldsymbol{C})&=f_b(\boldsymbol{C})\boldsymbol{C}-\boldsymbol{I},\qquad \boldsymbol{S}=\frac{\boldsymbol{T}(\boldsymbol{C})}{\mathrm{Wi}}. \tag{4} \end{align}\] The normalization \(f_b(\boldsymbol{I})=1\) makes \(\boldsymbol{C}=\boldsymbol{I}\) the equilibrium conformation. The parameter \(\varepsilon=\mathrm{Pe}_p^{-1}\) represents polymer center-of-mass molecular diffusion in the constant-polymer-density FENE-P closure. Periodic boundary conditions or the no-flux condition \[\nabla\boldsymbol{C}\,\boldsymbol{n}=0\] are assumed for this diffusive flux. If the polymer number density is not constant, the closure must be augmented by a concentration equation; that extension is outside the present analysis. The point here is that molecular diffusion must be paired with the FENE barrier entropy rather than treated as an arbitrary componentwise smoothing operator.

The elastic entropy is \[\Phi_b(\boldsymbol{C}) = -\log\det\boldsymbol{C}-(b-d)\log\left(\frac{b-\operatorname{tr}\boldsymbol{C}}{b-d}\right). \label{eq:fene-entropy}\tag{5}\] It satisfies \(\Phi_b(\boldsymbol{I})=0\) and \[D\Phi_b(\boldsymbol{C})[\boldsymbol{B}] = \left(f_b(\boldsymbol{C})\boldsymbol{I}-\boldsymbol{C}^{-1}\right):\boldsymbol{B} =\boldsymbol{H}_b(\boldsymbol{C}):\boldsymbol{B}.\] The Hessian is positive: \[D^2\Phi_b(\boldsymbol{C})[\boldsymbol{B},\boldsymbol{B}] = \operatorname{tr}(\boldsymbol{C}^{-1}\boldsymbol{B}\boldsymbol{C}^{-1}\boldsymbol{B}) + \frac{b-d}{(b-\operatorname{tr}\boldsymbol{C})^2}\bigl(\operatorname{tr}\boldsymbol{B}\bigr)^2 . \label{eq:hessian}\tag{6}\] Thus \(\Phi_b\) is convex on \(\mathcal{D}_b\) and blows up at both parts of the boundary: \(\lambda_{\min}(\boldsymbol{C})\downarrow0\) and \(\operatorname{tr}\boldsymbol{C}\uparrow b\).

  Let \(\boldsymbol{C}\in\mathcal{D}_b\) and \(\Phi_b(\boldsymbol{C})\le M\). Then there is a positive function \(\delta_b(M)\) such that \[b-\operatorname{tr}\boldsymbol{C}\ge \delta_b(M)>0 .\] One possible choice is \[\delta_b(M) = (b-d)\exp\!\left[-\frac{M+d\log(b/d)}{b-d}\right].\]

  Since \(\det\boldsymbol{C}\le(\operatorname{tr}\boldsymbol{C}/d)^d\le(b/d)^d\), one has \(-\log\det\boldsymbol{C}\ge-d\log(b/d)\). From the entropy bound, \[-(b-d)\log\left(\frac{b-\operatorname{tr}\boldsymbol{C}}{b-d}\right) \le M+d\log(b/d),\] which gives the stated lower bound on \(b-\operatorname{tr}\boldsymbol{C}\). \(\square\)

  For every \(\boldsymbol{C}\in\mathcal{D}_b\), \[\boldsymbol{T}(\boldsymbol{C}):\boldsymbol{H}_b(\boldsymbol{C}) = \sum_{i=1}^d \lambda_i\left(f_b(\boldsymbol{C})-\lambda_i^{-1}\right)^2 \ge0,\] where \(\lambda_i\) are the eigenvalues of \(\boldsymbol{C}\).

  Because \(f_b(\boldsymbol{C})\) is a scalar function of \(\operatorname{tr}\boldsymbol{C}\), it commutes with \(\boldsymbol{C}\). Diagonalizing \(\boldsymbol{C}\) gives \[\bigl(f_b\lambda_i-1\bigr)\bigl(f_b-\lambda_i^{-1}\bigr) = \lambda_i\bigl(f_b-\lambda_i^{-1}\bigr)^2 .\] Summing over \(i\) proves the claim. \(\square\)

  Assume periodic boundary conditions or \(\nabla\boldsymbol{C}\,\boldsymbol{n}=0\) on \(\partial\Omega\). For every smooth admissible conformation field, \[\int_\Omega \Delta\boldsymbol{C}:\boldsymbol{H}_b(\boldsymbol{C})\,\,\mathrm dx = -\int_\Omega \mathcal{I}_b(\boldsymbol{C})\,\,\mathrm dx ,\] where \[\mathcal{I}_b(\boldsymbol{C}) = \sum_{j=1}^d \left[ \operatorname{tr}\!\left(\boldsymbol{C}^{-1}\partial_j\boldsymbol{C}\,\boldsymbol{C}^{-1}\partial_j\boldsymbol{C}\right) + \frac{b-d}{(b-\operatorname{tr}\boldsymbol{C})^2}\bigl(\partial_j\operatorname{tr}\boldsymbol{C}\bigr)^2 \right]\ge0 .\]

  Integrating by parts gives \[\int_\Omega \Delta\boldsymbol{C}:\boldsymbol{H}_b(\boldsymbol{C})\,\,\mathrm dx = -\sum_{j=1}^d\int_\Omega \partial_j\boldsymbol{C}:D\boldsymbol{H}_b(\boldsymbol{C})[\partial_j\boldsymbol{C}]\,\,\mathrm dx .\] Since \(D\boldsymbol{H}_b=D^2\Phi_b\), the Hessian formula 6 gives the stated expression. \(\square\)

Trace-singular FENE-family closures↩︎

Although the formulas below are written for the FENE-P coefficient \(f_b\), the same algebra applies to a larger FENE-family class. Let \[\boldsymbol{T}_g(\boldsymbol{C})=g(\operatorname{tr}\boldsymbol{C})\boldsymbol{C}-\boldsymbol{I},\qquad \boldsymbol{H}_g(\boldsymbol{C})=g(\operatorname{tr}\boldsymbol{C})\boldsymbol{I}-\boldsymbol{C}^{-1},\] where \(g\in C^1([d,b))\), \(g(d)=1\), \(g'(s)\ge0\), and \(g(s)\to\infty\) as \(s\uparrow b\). Define \[G_g(s)=\int_d^s g(r)\,\,\mathrm dr,\qquad \Phi_g(\boldsymbol{C})=-\log\det\boldsymbol{C}+G_g(\operatorname{tr}\boldsymbol{C}).\] The FENE-P entropy corresponds to \(g=f_b\) and \(G_g(s)=-(b-d)\log((b-s)/(b-d))\).

  For every \(\boldsymbol{C}\in\mathcal{D}_b\) and every symmetric \(\boldsymbol{B}\), \[\begin{align} D\Phi_g(\boldsymbol{C})[\boldsymbol{B}]&=\boldsymbol{H}_g(\boldsymbol{C}):\boldsymbol{B},\\ D^2\Phi_g(\boldsymbol{C})[\boldsymbol{B},\boldsymbol{B}] &= \operatorname{tr}(\boldsymbol{C}^{-1}\boldsymbol{B}\boldsymbol{C}^{-1}\boldsymbol{B})+g'(\operatorname{tr}\boldsymbol{C})(\operatorname{tr}\boldsymbol{B})^2 . \end{align}\] Moreover, \[\boldsymbol{T}_g(\boldsymbol{C}):\boldsymbol{H}_g(\boldsymbol{C}) = \sum_{i=1}^d \lambda_i\left(g(\operatorname{tr}\boldsymbol{C})-\lambda_i^{-1}\right)^2\ge0,\] and, for \(\nabla\!\cdot\boldsymbol{u}=0\), \[-\big((\nabla\boldsymbol{u})\boldsymbol{C}+\boldsymbol{C}(\nabla\boldsymbol{u})^T\big):\boldsymbol{H}_g(\boldsymbol{C}) = -2\boldsymbol{T}_g(\boldsymbol{C}):\nabla\boldsymbol{u}.\]

  The derivative and Hessian follow from the chain rule and the derivative of \(-\log\det\boldsymbol{C}\). Since \(g(\operatorname{tr}\boldsymbol{C})\) is scalar, \(\boldsymbol{T}_g\) and \(\boldsymbol{H}_g\) are diagonal in the eigenbasis of \(\boldsymbol{C}\), giving the relaxation identity. The stretching identity follows from \(\boldsymbol{C}:\nabla\boldsymbol{u}=\boldsymbol{C}:\boldsymbol{D}(\boldsymbol{u})\) and \(\boldsymbol{I}:\nabla\boldsymbol{u}=\nabla\!\cdot\boldsymbol{u}=0\). \(\square\)

This proposition is the structural reason for calling the method FENE-type rather than only FENE-P. Once a closure has the trace-singular entropy variable \(\boldsymbol{H}_g\), the barrier-log state space, entropy-compatible reconstruction, and discrete coupling argument below carry over with \(f_b\) replaced by \(g(\operatorname{tr}\boldsymbol{C})\). The numerical section uses FENE-P because it is the standard Peterlin closure and makes the finite-extensibility denominator explicit.

3 Continuous energy law↩︎

The FENE trace barrier is compatible with the usual polymeric-work cancellation. Testing 1 by \(\boldsymbol{u}\) gives the polymeric work term \[- (1-\beta)\int_\Omega \boldsymbol{S}:\nabla\boldsymbol{u}\,\,\mathrm dx = -\frac{1-\beta}{\mathrm{Wi}}\int_\Omega \boldsymbol{T}(\boldsymbol{C}):\nabla\boldsymbol{u}\,\,\mathrm dx .\] Testing 3 by \(\boldsymbol{H}_b(\boldsymbol{C})\) yields the stretching contribution \[-\big((\nabla\boldsymbol{u})\boldsymbol{C}+\boldsymbol{C}(\nabla\boldsymbol{u})^T\big):\boldsymbol{H}_b(\boldsymbol{C}).\] Since \(\boldsymbol{H}_b(\boldsymbol{C})=f_b\boldsymbol{I}-\boldsymbol{C}^{-1}\) and \(\nabla\!\cdot\boldsymbol{u}=0\), \[\begin{align} &-\big((\nabla\boldsymbol{u})\boldsymbol{C}+\boldsymbol{C}(\nabla\boldsymbol{u})^T\big):(f_b\boldsymbol{I}-\boldsymbol{C}^{-1})\\ &\qquad = -2 f_b\boldsymbol{C}:\nabla\boldsymbol{u}+2\boldsymbol{I}:\nabla\boldsymbol{u} = -2\boldsymbol{T}(\boldsymbol{C}):\nabla\boldsymbol{u}. \end{align}\] Multiplication by \((1-\beta)/(2\mathrm{Wi})\) therefore cancels the polymeric work exactly.

  Assume a smooth solution with \(\boldsymbol{C}(t,x)\in\mathcal{D}_b\) and periodic or entropy-compatible boundary conditions. Then \[\begin{align} \frac{\,\mathrm d}{\,\mathrm dt}\left[ \frac{\mathrm{Re}}{2}\left\|\boldsymbol{u}\right\|^2_{L^2} +\frac{1-\beta}{2\mathrm{Wi}}\int_\Omega\Phi_b(\boldsymbol{C})\,\,\mathrm dx \right] &+\beta\left\|\nabla\boldsymbol{u}\right\|_{L^2}^2\\ &+\frac{1-\beta}{2\mathrm{Wi}^2} \int_\Omega\boldsymbol{T}(\boldsymbol{C}):\boldsymbol{H}_b(\boldsymbol{C})\,\,\mathrm dx\\ &+\frac{(1-\beta)\varepsilon}{2\mathrm{Wi}} \int_\Omega \mathcal{I}_b(\boldsymbol{C})\,\,\mathrm dx = (\boldsymbol{f},\boldsymbol{u}), \end{align} \label{eq:continuous-energy}\qquad{(1)}\] where \(\mathcal{I}_b\) is the molecular-diffusion entropy density in Lemma 2.3.

  The kinetic-energy identity follows from the incompressibility cancellation of the convection term. The entropy derivative follows from \(D\Phi_b=\boldsymbol{H}_b\). The stretching contribution cancels the polymeric work as shown above. Relaxation gives the nonnegative term in Lemma 2.2, and polymer molecular diffusion gives the nonnegative barrier dissipation in Lemma 2.3. Combining the terms gives ?? . \(\square\)

The trace-barrier term in \(\mathcal{I}_b\) is the new contribution relative to Oldroyd–B. It penalizes gradients of \(\operatorname{tr}\boldsymbol{C}\) more strongly near the finite-extensibility boundary. This term is also the reason that a discrete diffusion operator should be paired with the same barrier entropy rather than treated as a componentwise Laplacian detached from the entropy calculation.

4 Barrier logarithmic parametrization↩︎

A standard log-conformation parametrization \(\boldsymbol{C}=\operatorname{Exp}\Psi\) enforces \(\boldsymbol{C}\succ0\) but does not enforce \(\operatorname{tr}\boldsymbol{C}<b\). We use instead the map \[\mathcal{B}_b(\Psi) = \frac{b\,\operatorname{Exp}\Psi}{1+\operatorname{tr}(\operatorname{Exp}\Psi)}, \qquad \Psi=\Psi^T . \label{eq:barrier-log-map}\tag{7}\] Then \(\mathcal{B}_b(\Psi)\in\mathcal{D}_b\) for every symmetric \(\Psi\), because \[\operatorname{tr}\mathcal{B}_b(\Psi) = \frac{b\,\operatorname{tr}(\operatorname{Exp}\Psi)}{1+\operatorname{tr}(\operatorname{Exp}\Psi)}<b .\] The map is onto \(\mathcal{D}_b\). If \(\boldsymbol{C}\in\mathcal{D}_b\), set \[\operatorname{Exp}\Psi=\frac{\boldsymbol{C}}{b-\operatorname{tr}\boldsymbol{C}}.\] Then \[\mathcal{B}_b(\Psi) = \frac{b\,\boldsymbol{C}/(b-\operatorname{tr}\boldsymbol{C})}{1+\operatorname{tr}\boldsymbol{C}/(b-\operatorname{tr}\boldsymbol{C})} = \boldsymbol{C}.\] Thus 7 is a natural finite-extensibility analogue of the log-conformation map.

  If the nonlinear algebraic solve is performed in a symmetric variable \(\Psi_h\) and the accepted conformation tensor is defined by \(\boldsymbol{C}_h=\mathcal{B}_b(\Psi_h)\) at entropy quadrature points, then \[\boldsymbol{C}_h(x_q)\succ0,\qquad \operatorname{tr}\boldsymbol{C}_h(x_q)<b\] at every quadrature point \(x_q\), independent of the Newton iterate size before acceptance.

  The matrix exponential is positive definite, and the scalar denominator in 7 is strictly positive. The trace calculation above gives the strict upper bound. \(\square\)

This parametrization is not required for the energy proof; one may also use a line-search Newton method that rejects candidates outside \(\mathcal{D}_b\). Its advantage is practical: the finite-extensibility constraint is built into the unknown. The proof below is written in the accepted conformation tensor, so it applies to either realization.

5 Entropy-compatible barrier-log reconstruction↩︎

The barrier-log map solves the geometric part of admissibility, but it does not by itself guarantee compatibility with the discrete FENE free energy. A high-order reconstruction can remain inside \(\mathcal{D}_b\) and still add artificial elastic entropy. We therefore add an entropy-compatible reconstruction layer.

Let \(\widehat\Psi_q\) be a physical predictor at quadrature point \(x_q\), and let \(\widetilde{\Psi}_q\) be a raw high-order barrier-log reconstruction. Define \[\Psi_q(\theta)=\widehat\Psi_q+\theta(\widetilde{\Psi}_q-\widehat\Psi_q), \qquad \boldsymbol{C}_q(\theta)=\mathcal{B}_b(\Psi_q(\theta)),\qquad 0\le\theta\le1.\] For a nonnegative budget \(\tau_h\), the accepted parameter is \[\theta_\star = \max\left\{\theta\in[0,1]: \sum_qw_q\Phi_b(\boldsymbol{C}_q(\theta)) \le \sum_qw_q\Phi_b(\boldsymbol{C}_q(0))+\tau_h \right\}. \label{eq:theta-star}\tag{8}\] The accepted tensor is \(\boldsymbol{C}_q^\star=\boldsymbol{C}_q(\theta_\star)\). If the raw reconstruction is already entropy-compatible, then \(\theta_\star=1\). Otherwise the method damps only the entropy-incompatible part along the barrier-log segment.

  Let \(\boldsymbol{C}(\theta)=\mathcal{B}_b(\Psi_0+\theta E)\). Then \(J(\theta)=\Phi_b(\boldsymbol{C}(\theta))\) is convex on \([0,1]\). More explicitly, \[\Phi_b(\mathcal{B}_b(\Psi)) = b\log\bigl(1+\operatorname{tr}(e^\Psi)\bigr)-\operatorname{tr}\Psi+\gamma_b, \label{eq:barrier-log-entropy}\qquad{(2)}\] where \(\gamma_b\) is independent of \(\Psi\).

  Set \(Z(\Psi)=1+\operatorname{tr}(e^\Psi)\). Since \(\mathcal{B}_b(\Psi)=b e^\Psi/Z\), one has \[\det\mathcal{B}_b(\Psi)=b^d e^{\operatorname{tr}\Psi}Z^{-d}, \qquad b-\operatorname{tr}\mathcal{B}_b(\Psi)=b/Z .\] Substitution into 5 gives ?? . The map \(\Psi\mapsto\log(1+\operatorname{tr}e^\Psi)\) is the spectral log-sum-exp function with an additional zero mode and is convex on symmetric matrices. The term \(-\operatorname{tr}\Psi\) is affine. Hence \(J\) is convex along every affine path. \(\square\)

  The set in 8 is a closed interval containing \(0\). Hence \(\theta_\star\) exists, is the largest admissible parameter, and can be computed by bisection. The accepted tensor satisfies \[\boldsymbol{C}_q^\star\in\mathcal{D}_b,\qquad \sum_qw_q\Phi_b(\boldsymbol{C}_q^\star) \le \sum_qw_q\Phi_b(\boldsymbol{C}_q(0))+\tau_h .\]

  Admissibility in \(\mathcal{D}_b\) follows from Proposition 4.1 for every \(\theta\in[0,1]\). By Lemma 5.1, the entropy profile is continuous and convex; therefore its sublevel set is a closed interval. Since \(\theta=0\) satisfies the inequality, the interval is nonempty and has a largest element. Bisection applies because admissibility is monotone along this interval. \(\square\)

  Assume the raw barrier-log reconstruction defect \(E_q=\widetilde{\Psi}_q-\widehat\Psi_q\) is \(\mathcal{O}(h^{k+1})\) and satisfies \[\sum_qw_q\,D\{\Phi_b\circ\mathcal{B}_b\}(\widehat\Psi_q)[E_q] = \mathcal{O}(h^{2k+2}).\] If \(\tau_h=c_\tau h^{2k+2}\) with \(c_\tau\) large enough for the leading consistency constant, then \(\theta_\star=1\) for sufficiently small \(h\). In the marginal active case, \(1-\theta_\star=\mathcal{O}(h^{k+1})\) under a nondegenerate endpoint derivative.

  Taylor expansion of the convex entropy profile gives \[\sum_qw_q\{\Phi_b(\boldsymbol{C}_q(1))-\Phi_b(\boldsymbol{C}_q(0))\} = \mathcal{O}(h^{2k+2})\] under the stated first-variation condition and compactness of the path. The mesh-scaled budget therefore accepts the raw endpoint for sufficiently small \(h\). If the endpoint is marginally active, convexity and the nondegenerate derivative give the stated bound on \(1-\theta_\star\). \(\square\)

When \(\boldsymbol{C}^\star\) is used consistently in the stress force, stretching term, relaxation term, and entropy quadrature, the fully discrete energy estimate below is unchanged except for the explicit budget contribution \((1-\beta)\tau_h/(2\mathrm{Wi})\). Thus the reconstruction is not merely trace-preserving; it is compatible with the FENE free-energy balance.

6 Fully discrete scheme↩︎

Let \(V_h\times Q_h\) be an inf-sup stable velocity–pressure pair and \(M_h\) a symmetric tensor space. Let \((\cdot,\cdot)_Q\) denote a positive quadrature rule used consistently in the entropy terms and the polymeric work. At time level \(n+1\), find \[(\boldsymbol{u}_h^{n+1},p_h^{n+1},\boldsymbol{C}_h^{n+1})\in V_h\times Q_h\times M_h\] such that \(\boldsymbol{C}_h^{n+1}(x_q)\in\mathcal{D}_b\) and \[\boldsymbol{T}_h^{n+1}=f_b(\boldsymbol{C}_h^{n+1})\boldsymbol{C}_h^{n+1}-\boldsymbol{I},\qquad \boldsymbol{S}_h^{n+1}=\frac{\boldsymbol{T}_h^{n+1}}{\mathrm{Wi}}\] at quadrature points. The momentum step is \[\begin{align} \mathrm{Re}\left(\frac{\boldsymbol{u}_h^{n+1}-\boldsymbol{u}_h^n}{\Delta t},\boldsymbol{v}_h\right) &+\mathrm{Re}\,c_h(\boldsymbol{u}_h^{n+1};\boldsymbol{u}_h^{n+1},\boldsymbol{v}_h) +\beta(\nabla\boldsymbol{u}_h^{n+1},\nabla\boldsymbol{v}_h)\\ &-(p_h^{n+1},\nabla\!\cdot\boldsymbol{v}_h) -(1-\beta)(\boldsymbol{S}_h^{n+1},\nabla\boldsymbol{v}_h)_Q = (\boldsymbol{f}^{n+1},\boldsymbol{v}_h), \end{align} \label{eq:scheme-mom}\tag{9}\] with \[(q_h,\nabla\!\cdot\boldsymbol{u}_h^{n+1})=0. \label{eq:scheme-div}\tag{10}\] The conformation update is \[\begin{align} \left(\frac{\boldsymbol{C}_h^{n+1}-\boldsymbol{C}_h^n}{\Delta t},\boldsymbol{B}_h\right)_Q &+a_h(\boldsymbol{u}_h^{n+1};\boldsymbol{C}_h^{n+1},\boldsymbol{B}_h)\\ &-\left((\nabla\boldsymbol{u}_h^{n+1})\boldsymbol{C}_h^{n+1} +\boldsymbol{C}_h^{n+1}(\nabla\boldsymbol{u}_h^{n+1})^T,\boldsymbol{B}_h\right)_Q\\ &+\frac{1}{\mathrm{Wi}}(\boldsymbol{T}_h^{n+1},\boldsymbol{B}_h)_Q +\varepsilon d_h(\boldsymbol{C}_h^{n+1},\boldsymbol{B}_h)=0 . \end{align} \label{eq:scheme-conf}\tag{11}\] The forms are chosen so that the velocity convection is skew-symmetric, conformation transport has zero entropy contribution for discretely incompressible velocity, and the stretching term is evaluated at the same quadrature points as the stress force. The form \(d_h\) is the discrete polymer molecular-diffusion operator and is paired with the barrier entropy variable: \[d_h(\boldsymbol{C}_h,\boldsymbol{H}_b(\boldsymbol{C}_h))=\mathcal{I}_{b,h}(\boldsymbol{C}_h)\ge0\] for conforming periodic or no-flux discretizations, or at least \(d_h(\boldsymbol{C}_h,\boldsymbol{H}_b(\boldsymbol{C}_h))\ge \mathcal{I}_{b,h}(\boldsymbol{C}_h)\) for stabilized variants. This is the discrete counterpart of Lemma 2.3. If a high-order reconstruction is used before the conformation state enters the coupled step, \(\boldsymbol{C}_h^{n+1}\) denotes the entropy-compatible accepted tensor from Section 5.

Define the discrete energy \[\mathcal{E}_h^n = \frac{\mathrm{Re}}{2}\left\|\boldsymbol{u}_h^n\right\|_{L^2}^2 + \frac{1-\beta}{2\mathrm{Wi}} \sum_q w_q\Phi_b(\boldsymbol{C}_h^n(x_q)).\]

  Assume the accepted conformation tensor satisfies \(\boldsymbol{C}_h^{n+1}(x_q)\in\mathcal{D}_b\) at all entropy quadrature points and the discrete transport and stretching forms have the compatibility properties stated above. Then \[\begin{align} \mathcal{E}_h^{n+1}-\mathcal{E}_h^n &+\frac{\mathrm{Re}}{2}\left\|\boldsymbol{u}_h^{n+1}-\boldsymbol{u}_h^n\right\|_{L^2}^2 +\Delta t\,\beta\left\|\nabla\boldsymbol{u}_h^{n+1}\right\|_{L^2}^2\\ &+\Delta t\,\frac{1-\beta}{2\mathrm{Wi}^2} \sum_qw_q\,\boldsymbol{T}_h^{n+1}(x_q):\boldsymbol{H}_b(\boldsymbol{C}_h^{n+1}(x_q))\\ &+\Delta t\,\frac{(1-\beta)\varepsilon}{2\mathrm{Wi}} \mathcal{I}_{b,h}(\boldsymbol{C}_h^{n+1}) \le \Delta t\,(\boldsymbol{f}^{n+1},\boldsymbol{u}_h^{n+1}) +\frac{1-\beta}{2\mathrm{Wi}}\tau_h . \end{align} \label{eq:discrete-energy}\qquad{(3)}\] In particular, if the right-hand side is controlled by Young’s inequality, the discrete entropy remains bounded, and Lemma 2.1 gives a positive trace buffer at quadrature points.

  Set \(\boldsymbol{v}_h=\boldsymbol{u}_h^{n+1}\) in 9 . The skew-symmetric convection and pressure terms vanish, and the standard identity \[(a-b,a)=\frac{1}{2}\bigl(\left\|a\right\|^2-\left\|b\right\|^2+\left\|a-b\right\|^2\bigr)\] gives the kinetic increment. In 11 use the entropy test \[\boldsymbol{B}_h=\boldsymbol{H}_b(\boldsymbol{C}_h^{n+1}) = f_b(\boldsymbol{C}_h^{n+1})\boldsymbol{I}-(\boldsymbol{C}_h^{n+1})^{-1}\] at quadrature points. Convexity of \(\Phi_b\) gives \[\boldsymbol{H}_b(\boldsymbol{C}_h^{n+1}):(\boldsymbol{C}_h^{n+1}-\boldsymbol{C}_h^n) \ge \Phi_b(\boldsymbol{C}_h^{n+1})-\Phi_b(\boldsymbol{C}_h^n).\] If the conformation tensor is obtained by the entropy-compatible reconstruction of Section 5, the accepted endpoint contributes at most the budget \(\tau_h\) to the quadrature entropy; otherwise \(\tau_h=0\). The stretching term cancels the polymeric work in the momentum equation by the same pointwise identity used in the continuous proof. Relaxation is nonnegative by Lemma 2.2, and diffusion contributes the discrete Hessian dissipation \(\mathcal{I}_{b,h}\). Summing the momentum and entropy relations gives ?? . \(\square\)

7 AP stress closure and Newtonian limit↩︎

The finite-extensibility stress is nonlinear in \(\boldsymbol{C}\), but the AP variable is simple: \[\boldsymbol{S}=\frac{f_b(\boldsymbol{C})\boldsymbol{C}-\boldsymbol{I}}{\mathrm{Wi}}.\] The conformation equation can be rewritten as a relaxation equation for \(\boldsymbol{S}\). At fixed \(h\) and \(\Delta t\), insert \(\boldsymbol{T}_h^{n+1}=\mathrm{Wi}\boldsymbol{S}_h^{n+1}\) into 11 . The relaxation term becomes \((\boldsymbol{S}_h^{n+1},\boldsymbol{B}_h)_Q\), while the leading stretching term gives the rate-of-strain tensor. All time-difference, transport, diffusion, and nonlinear stretching remainders are multiplied by \(\mathrm{Wi}\) after the stress substitution.

  Assume the discrete states remain uniformly bounded in the norms entering the conformation residual and stay in a compact subset of \(\mathcal{D}_b\). Then, for fixed \(h\) and \(\Delta t\), \[\left\|\boldsymbol{S}_h^{n+1}-\boldsymbol{D}(\boldsymbol{u}_h^{n+1})\right\|_{M_h'} \le C\mathrm{Wi}, \label{eq:ap-closure}\qquad{(4)}\] where \(C\) is independent of \(\mathrm{Wi}\) in the tested small-Weissenberg range.

  After substituting \(\boldsymbol{T}_h^{n+1}=\mathrm{Wi}\boldsymbol{S}_h^{n+1}\) in 11 , the leading balance is \[(\boldsymbol{S}_h^{n+1},\boldsymbol{B}_h)_Q = (\boldsymbol{D}(\boldsymbol{u}_h^{n+1}),\boldsymbol{B}_h)_Q + \mathrm{Wi}\,\mathcal{R}_h^{n+1}(\boldsymbol{B}_h),\] with the normalization of \(\boldsymbol{D}\) matching the stretching convention. The residual \(\mathcal{R}_h^{n+1}\) contains the time difference, transport, diffusion, and higher-order stretching terms. The assumed uniform bound gives \[|\mathcal{R}_h^{n+1}(\boldsymbol{B}_h)|\le C\left\|\boldsymbol{B}_h\right\|_{M_h},\] which proves ?? . \(\square\)

  Let \(\mathrm{Wi}_j\to0\) and assume the corresponding admissible discrete solutions are uniformly bounded and remain in a compact subset of \(\mathcal{D}_b\). Then any convergent subsequence has a limit satisfying the discrete incompressible Navier–Stokes scheme with total viscosity equal to the solvent viscosity plus the polymeric contribution from the closure \(\boldsymbol{S}_h=\boldsymbol{D}(\boldsymbol{u}_h)\).

  The finite-dimensional compactness gives a convergent subsequence. The AP closure identifies the stress limit as \(\boldsymbol{D}(\boldsymbol{u}_h)\). Passing to the limit in the momentum equation gives the discrete Newtonian system with the polymeric stress absorbed into the viscous operator. \(\square\)

8 Relative-entropy error estimate↩︎

For convergence rates one needs a compact spectral and trace set \[\mathcal{K}_{\lambda,\delta} = \{\boldsymbol{C}\in\mathcal{D}_b:\lambda\boldsymbol{I}\preceq\boldsymbol{C},\;b-\operatorname{tr}\boldsymbol{C}\ge\delta\}, \qquad \lambda,\delta>0 .\] On this set the entropy Hessian is bounded above and below, so the relative entropy \[\Phi_b(\boldsymbol{C}|\widehat\boldsymbol{C}) = \Phi_b(\boldsymbol{C})-\Phi_b(\widehat\boldsymbol{C}) -\boldsymbol{H}_b(\widehat\boldsymbol{C}):(\boldsymbol{C}-\widehat\boldsymbol{C})\] is equivalent to \(\left\|\boldsymbol{C}-\widehat\boldsymbol{C}\right\|_F^2\). The constants deteriorate as \(\lambda\downarrow0\) or \(\delta\downarrow0\), which is unavoidable because the entropy becomes singular at both boundaries.

  Let \((\boldsymbol{u},p,\boldsymbol{C})\) be a smooth solution on \([0,T]\) with \(\boldsymbol{C}(t,x)\in\mathcal{K}_{\lambda,\delta}\). Assume stable projections, a quadrature-consistent discretization of the coupling terms, and nonlinear solver residuals of order \(h^{k+1}+\Delta t\). If the discrete solution remains in a slightly larger compact subset of \(\mathcal{D}_b\), then \[\max_{0\le n\le N} \left[ \left\|\boldsymbol{u}(t_n)-\boldsymbol{u}_h^n\right\|_{L^2}^2 + \sum_qw_q \Phi_b(\boldsymbol{C}(t_n,x_q)|\boldsymbol{C}_h^n(x_q)) \right] \le C_T(h^{2k+2}+\Delta t^2).\]

  The proof follows the relative-entropy stability argument. Subtract the projected exact equations from the discrete equations, test the velocity error by itself, and test the conformation error in the entropy variable. The compatible quadrature cancels the leading polymeric-work error. The Hessian bounds on \(\mathcal{K}_{\lambda,\delta}\) control the nonlinear FENE coefficients and convert relative entropy into squared tensor error. Consistency and solver residuals contribute \(C(h^{2k+2}+\Delta t^2)\), while lower-order terms are bounded by the current error. A discrete Gronwall inequality gives the estimate. \(\square\)

The theorem is intentionally conditional: it gives a classical rate in a resolved regime away from the finite-extensibility boundary. The energy theorem and trace-barrier buffer are the mechanisms that help keep the computation in such a regime.

9 Numerical verification↩︎

The numerical diagnostics are designed to test the mechanisms used in the analysis: trace-barrier preservation, entropy-compatible reconstruction, energy decay, AP closure, coupled velocity–pressure–stress feedback, and high-Weissenberg behavior near the finite-extensibility boundary. The reported benchmarks are kept compact and tied directly to the structural claims.

9.1 Barrier preservation, entropy correction, and energy decay↩︎

The first diagnostic compares a plain log-conformation map with the barrier-log map 7 . The background dimension is \(d=2\) and \(L^2=12\). The raw log map preserves positive definiteness but can produce \(\operatorname{tr}\boldsymbol{C}>L^2\) for large logarithmic stretches. The barrier-log map keeps the same unconstrained variable while mapping it into \(\mathcal{D}_b\).

The next diagnostic isolates the reconstruction layer and compares it with trace-only denominator repair. We start from a near-barrier tensor with \(\operatorname{tr}\boldsymbol{C}/b=0.88\) and apply the same symmetric high-order perturbation in different variables. The standard log state can cross the trace barrier immediately. A trace repair rescales the state back below \(b\), which mimics the role of trace-first safeguards in denominator control, but it leaves a very large FENE entropy jump because the repaired state is placed extremely close to the singular barrier. The barrier-log endpoint remains admissible by construction, yet it may still add substantial FENE entropy. The entropy-compatible barrier-log reconstruction selects the largest admissible \(\theta\) satisfying 8 , so the accepted state controls both the trace barrier and the entropy increment.

4pt

Next we evolve a homogeneous relaxation-diffusion diagnostic with no forcing. Figure 1 shows the normalized free energy and the maximum trace ratio. The free energy decreases, while the trace remains below the finite-extensibility threshold.

9.2 AP closure↩︎

The AP test fixes \(\Delta t=0.05\) and varies \(\mathrm{Wi}\). The error is the dual norm of \(\boldsymbol{S}_h-\boldsymbol{D}(\boldsymbol{u}_h)\) at the final time. The observed slope is essentially first order in \(\mathrm{Wi}\), as predicted by the AP closure estimate ?? .

9.3 Fully coupled manufactured solution↩︎

The manufactured solution uses a periodic velocity, mean-zero pressure, and conformation tensor defined through the barrier map \[\boldsymbol{C}_{\rm ex}(t,x)=\mathcal{B}_b(\Psi_{\rm ex}(t,x)).\] This construction guarantees \(\boldsymbol{C}_{\rm ex}\in\mathcal{D}_b\) and allows the manufactured forcing to exercise the nonlinear FENE stress, pressure projection, tensor transport, stretching, relaxation, and diffusion terms simultaneously. Table 8 reports the observed errors at final time \(T=0.05\).

9.4 High-Weissenberg robustness near the trace barrier↩︎

The final benchmark stresses the finite-extensibility boundary. We prescribe an extensional velocity field and initialize \(\boldsymbol{C}\) so that \(\operatorname{tr}\boldsymbol{C}/L^2\) is already large. The test is run at increasing \(\mathrm{Wi}\) with \(L^2=20\), \(\varepsilon=2\times10^{-3}\), and \(T=0.12\). Larger \(\mathrm{Wi}\) weakens relaxation and increases stretch, but the barrier update keeps the conformation tensor inside \(\mathcal{D}_b\).

The last row shows one line-search rejection when \(\mathrm{Wi}=200\). The accepted state remains strictly within the barrier. This diagnostic is useful because a positivity-only method would not detect the finite-extensibility failure until the FENE coefficient becomes singular or changes sign numerically.

10 Discussion↩︎

The FENE trace constraint changes what a structure-preserving method must prove. Positivity is still necessary, but it is no longer the complete admissibility condition. The stress coefficient \(f_b(\boldsymbol{C})\) depends on the distance to the trace boundary, and the polymer molecular-diffusion term is entropy dissipative only when it is paired with the FENE barrier variable. A stable method must therefore preserve the accepted tensor in the set \(\mathcal{D}_b\), control reconstruction-level FENE entropy increments, and use that same tensor consistently in the stress force, stretching term, relaxation term, molecular diffusion, and entropy quadrature.

The barrier-log map 7 is one practical way to enforce the finite-extensibility geometry. It is a direct analogue of the log-conformation map, but its range is the FENE admissible set rather than the entire positive definite cone. The entropy-compatible reconstruction 8 adds the missing balance condition: among admissible barrier-log states, it chooses the least-damping state whose FENE entropy increment is within the prescribed budget. The energy proof itself is independent of the nonlinear parametrization; it requires only that the accepted tensor be admissible, entropy-compatible, and used consistently in the discrete coupling identities.

The AP part of the analysis also differs from Oldroyd–B in detail. The conformation perturbation \(\boldsymbol{C}-\boldsymbol{I}\) is not the best stress variable because the FENE factor changes with \(\operatorname{tr}\boldsymbol{C}\). The regular stress variable is instead \(\boldsymbol{T}(\boldsymbol{C})/\mathrm{Wi}\). This choice avoids singular scaling in the momentum equation and gives the correct Newtonian limit at fixed discretization parameters.

11 Conclusion↩︎

We have proposed and analyzed a barrier-preserving, entropy-compatible, asymptotic-preserving discretization for molecularly diffusive FENE-P type viscoelastic flow. The key finite-extensibility feature is the admissible set \(\mathcal{D}_b=\{\boldsymbol{C}\succ0,\operatorname{tr}\boldsymbol{C}<b\}\), which is enforced by a trace-barrier entropy and can be built into the nonlinear solve through a barrier logarithmic parametrization. High-order barrier-log reconstructions are filtered by a least-damping FENE entropy constraint, so admissibility and entropy compatibility are enforced simultaneously. Compatible quadrature makes the cancellation between FENE polymeric work and conformation stretching exact at the fully discrete level, while polymer molecular diffusion contributes the barrier Hessian dissipation. The resulting scheme satisfies a free-energy inequality with relaxation and molecular-diffusion barrier dissipation, preserves a positive trace buffer under an entropy bound, and has a regular small-Weissenberg limit in the scaled FENE stress variable. The numerical diagnostics verify the main mechanisms in both moderate and near-barrier regimes.

Acknowledgments↩︎

The author acknowledges financial support from the National Natural Science Foundation of China (NSFC, Grant No. 12501602), the Education Department of Hunan Province (Grant No. 24C0055), the Science and Technology Department of Hunan Province (Grant No. 2025JJ60052), and the Scientific Research Start-up Fund of Xiangtan University (Grant No. KZ0810769).

References↩︎