2025-11-12T15:04:10.606403

A decoupled Crank-Nicolson leap-frog scheme for the unsteady bioconvection flows problem with concentration dependent viscosity

Li
A fully discrete Crank--Nicolson Leap--Frog (CNLF) scheme is proposed and analyzed for the unsteady bioconvection flow problem with concentration-dependent viscosity. Spatial discretization is handled via the Galerkin finite element method (FEM), while temporal discretization employs the CNLF method for the linear terms and a semi-implicit approach for the nonlinear terms. The scheme is proven to be unconditionally stable, i.e., the time step is not subject to a restrictive upper bound. Using the energy method, $L^2$-optimal error estimates are derived for the velocity and concentration . Finally, numerical experiments are presented to validate the theoretical results.
academic

A decoupled Crank-Nicolson leap-frog scheme for the unsteady bioconvection flows problem with concentration dependent viscosity

Basic Information

  • Paper ID: 2510.14034
  • Title: A decoupled Crank-Nicolson leap-frog scheme for the unsteady bioconvection flows problem with concentration dependent viscosity
  • Author: Chenyang Li (School of Mathematical Sciences, East China Normal University)
  • Classification: math.NA cs.NA
  • Publication Date: October 15, 2025 (arXiv preprint)
  • Paper Link: https://arxiv.org/abs/2510.14034

Abstract

This paper proposes and analyzes a fully discrete Crank-Nicolson leap-frog (CNLF) scheme for the unsteady bioconvection flow problem with concentration-dependent viscosity. The spatial discretization employs the Galerkin finite element method, while the temporal discretization applies the CNLF method to linear terms and a semi-implicit method to nonlinear terms. The scheme is proven to be unconditionally stable, meaning the time step is not subject to restrictive upper bounds. Using energy methods, optimal L2L^2 error estimates for velocity and concentration are derived. Finally, numerical experiments validate the theoretical results.

Research Background and Motivation

Problem Background

Bioconvection refers to the fluid convection phenomenon induced by microorganism motion, which has significant importance in biology, environmental science, and engineering applications. This phenomenon is described by coupled Navier-Stokes type equations and convection-diffusion equations:

  1. Fluid dynamics equations: Describing the flow of incompressible viscous culture fluid
  2. Microorganism transport equations: Describing the transport process of microorganisms

Core Challenges

  1. Concentration-dependent viscosity: Unlike classical Newtonian fluids, the viscosity of real suspensions depends on microorganism concentration
  2. Nonlinear coupling: Strong coupling exists between velocity and concentration fields
  3. Numerical stability: Requires designing numerical schemes that are both stable and efficient

Research Motivation

Existing methods have the following limitations when dealing with bioconvection problems with concentration-dependent viscosity:

  • Most studies assume constant viscosity
  • Existing numerical schemes may require strict time step restrictions
  • Lack of optimal error estimates for variable viscosity cases

Core Contributions

  1. Proposed a fully discrete CNLF scheme: First application of the Crank-Nicolson leap-frog method to bioconvection problems with concentration-dependent viscosity
  2. Proved unconditional stability: Time step is not subject to restrictive upper bounds, improving the practicality of the numerical scheme
  3. Established optimal error estimates: Achieved optimal convergence orders for velocity and concentration in the L2L^2 norm
  4. Provided a decoupled algorithm: Semi-implicit treatment requires solving only linear systems at each time step, improving computational efficiency

Method Details

Problem Formulation

Consider the bioconvection model on a bounded domain ΩRd\Omega \subset \mathbb{R}^d (d=2d=2 or 33):

utdiv(ν(c)D(u))+uu+p=g(1+γc)i2+f\frac{\partial u}{\partial t} - \text{div}(\nu(c)D(u)) + u \cdot \nabla u + \nabla p = -g(1+\gamma c)i_2 + f

u=0\nabla \cdot u = 0

ctθΔc+uc+Ucx2=0\frac{\partial c}{\partial t} - \theta\Delta c + u \cdot \nabla c + U\frac{\partial c}{\partial x_2} = 0

where:

  • uu: velocity field, pp: pressure, cc: concentration field
  • ν(c)\nu(c): concentration-dependent viscosity function
  • D(u)=12(u+uT)D(u) = \frac{1}{2}(\nabla u + \nabla u^T): strain rate tensor

Model Architecture

1. Spatial Discretization

Mixed finite element method is employed:

  • Velocity-pressure: Mini element (P1b-P1)
  • Concentration: Piecewise linear element (P1)

Finite element spaces are defined as: Vh={vhC(Ω)2VvhK(P1(K)b(K))2,KTh}V_h = \{v_h \in C(\Omega)^2 \cap V | v_h|_K \in (P_1(K) \oplus b(K))^2, \forall K \in T_h\}Mh={qhC(Ω)H1(Ω)qhKP1(K),KTh,Ωqhdx=0}M_h = \{q_h \in C(\Omega) \cap H^1(\Omega) | q_h|_K \in P_1(K), \forall K \in T_h, \int_\Omega q_h dx = 0\}

2. Temporal Discretization: CNLF Scheme

First step (Backward Euler): uh1uh0τ+ν(ch0+α)uh1+B(uh0,uh1,vh)(vh,ph1)=RHS\frac{u_h^1 - u_h^0}{\tau} + \nu(c_h^0 + \alpha)\nabla u_h^1 + B(u_h^0, u_h^1, v_h) - (\nabla \cdot v_h, p_h^1) = \text{RHS}

Subsequent steps (CNLF): uhn+1uhn12τ+A(chn,uhn+1+uhn12,vh)+B(uhn,uhn+1+uhn12,vh)=RHS\frac{u_h^{n+1} - u_h^{n-1}}{2\tau} + A(c_h^n, \frac{u_h^{n+1} + u_h^{n-1}}{2}, v_h) + B(u_h^n, \frac{u_h^{n+1} + u_h^{n-1}}{2}, v_h) = \text{RHS}

Technical Innovations

  1. Decoupling strategy: Achieves decoupling of velocity and concentration equations through semi-implicit treatment of nonlinear terms
  2. Leap-frog temporal integration: Applies second-order accurate Crank-Nicolson scheme to linear terms
  3. Variable coefficient handling: Specially designed projection operators for handling concentration-dependent viscosity

Experimental Setup

Computational Domain

Domain Ω=[0,1]×[0,1]\Omega = [0,1] \times [0,1] with parameter settings:

  • θ=γ=1\theta = \gamma = 1
  • Final time T=1.0T = 1.0
  • Analytical solutions: u(x,y,t)=(yet(2y1)(y1),xet(2x1)(x1))Tu(x,y,t) = (ye^{-t}(2y-1)(y-1), -xe^{-t}(2x-1)(x-1))^Tp(x,y,t)=et(2x1)(2y1)p(x,y,t) = e^{-t}(2x-1)(2y-1)c(x,y,t)=etsin(πx)sin(πy)c(x,y,t) = e^{-t}\sin(\pi x)\sin(\pi y)

Evaluation Metrics

  • L2L^2 norm error: rrhL2=r(tN)rhNL2\|r - r_h\|_{L^2} = \|r(t_N) - r_h^N\|_{L^2}
  • H1H^1 norm error: rrhH1\|r - r_h\|_{H^1}
  • Convergence rate: Computed through mesh refinement

Comparison Methods

Three different viscosity models are considered:

  1. ν=1\nu = 1 (constant viscosity)
  2. ν=1+0.1c\nu = 1 + 0.1c (linear dependence)
  3. ν=ec\nu = e^c (exponential dependence)

Implementation Details

  • Time step: τ=h\tau = h
  • Mesh refinement: h=1/4,1/8,1/16,1/32,1/64,1/128h = 1/4, 1/8, 1/16, 1/32, 1/64, 1/128
  • Implementation tool: FreeFEM++

Experimental Results

Main Results

Stability Verification

For all three viscosity models, numerical solutions remain stable across different mesh scales, verifying the unconditional stability of the scheme.

Convergence Analysis

L2L^2 norm convergence:

  • Velocity: Achieves second-order convergence for all viscosity models
  • Concentration: Achieves second-order convergence
  • Pressure: Achieves first-order convergence

Specific numerical results (with ν=1\nu = 1 as example):

hhuuhL2\|u-u_h\|_{L^2}RatecchL2\|c-c_h\|_{L^2}RatepphL2\|p-p_h\|_{L^2}Rate
1/40.0087769-0.0182156-0.033836-
1/80.0022631.960.00888621.040.01309761.37
1/160.00062861.850.0023941.890.00711040.88
1/320.00016641.920.0006031.990.00366560.96

Ablation Studies

Through comparison of different viscosity models, the following are verified:

  1. Robustness of the CNLF scheme to different viscosity functions
  2. Concentration-dependent viscosity does not affect scheme convergence
  3. Theoretical predicted convergence rates match numerical results

Experimental Findings

  1. Optimal convergence: Numerical experiments fully verify the theoretical O(τ2+h2)O(\tau^2 + h^2) convergence rate
  2. Robustness: The scheme exhibits good stability and convergence for different types of viscosity functions
  3. Efficiency advantages: The decoupled scheme significantly improves computational efficiency

Theoretical Analysis

Stability Analysis

Theorem 3.1 (Unconditional Stability): uhn+1L22+chn+1L22+κτn=1N(uhn+1+uhn1)L22C\|u_h^{n+1}\|_{L^2}^2 + \|c_h^{n+1}\|_{L^2}^2 + \kappa\tau\sum_{n=1}^N \|\nabla(u_h^{n+1} + u_h^{n-1})\|_{L^2}^2 \leq C

The proof employs energy methods with key steps:

  1. Utilizing skew-symmetric properties to handle nonlinear terms
  2. Applying discrete Gronwall inequality

Error Estimates

Theorem 4.1 (Convergence): Under assumptions A1 and A2, there exists a constant CC such that: max0iN(uiuhiL22+cichiL22)C(τ4+h4)\max_{0 \leq i \leq N}(\|u^i - u_h^i\|_{L^2}^2 + \|c^i - c_h^i\|_{L^2}^2) \leq C(\tau^4 + h^4)

The proof employs mathematical induction combined with:

  1. Error estimates of projection operators
  2. Truncation error analysis of temporal discretization
  3. Refined treatment of nonlinear terms

Numerical Methods for Bioconvection

  1. Constant viscosity case: 23,24 established existence of solutions, 27 provided finite element error estimates
  2. Variable viscosity case: 26 proved existence and uniqueness of weak solutions, 9 proposed BDF2 scheme
  3. Higher-order methods: 19 developed linearized Crank-Nicolson scheme
  1. First CNLF application: Introduces leap-frog method to bioconvection problems
  2. Unconditional stability: Eliminates time step restrictions compared to existing methods
  3. Decoupled design: Improves computational efficiency and facilitates parallel implementation

Conclusions and Discussion

Main Conclusions

  1. Successfully constructed a CNLF scheme for bioconvection problems with concentration-dependent viscosity
  2. Theoretically proved unconditional stability and optimal convergence of the scheme
  3. Numerical experiments validated the correctness of theoretical results

Limitations

  1. Dimensional restriction: Theoretical analysis primarily focuses on two-dimensional cases
  2. Viscosity function assumptions: Requires Lipschitz continuity and boundedness conditions
  3. Boundary conditions: Only considers homogeneous Dirichlet boundary conditions

Future Directions

The authors propose extending the CNLF framework to:

  1. Chemotaxis-Navier-Stokes systems
  2. Patlak-Keller-Segel-Navier-Stokes systems
  3. Chemo-Repulsion-Navier-Stokes systems

In-Depth Evaluation

Strengths

  1. Theoretical rigor: Complete stability and convergence analysis with detailed proofs
  2. Methodological innovation: First application of CNLF method to variable viscosity bioconvection problems
  3. Practical value: Unconditional stability makes the scheme more flexible in practical applications
  4. Comprehensive numerical verification: Testing with multiple viscosity models validates method robustness

Weaknesses

  1. Strong theoretical assumptions: High regularity requirements for solutions may limit practical applicability
  2. Missing three-dimensional extension: Theoretical analysis primarily limited to two dimensions
  3. Insufficient computational complexity analysis: Lacks comparison of computational efficiency with other methods
  4. Limited physical parameter sensitivity analysis: Insufficient discussion of method sensitivity to physical parameter variations

Impact

  1. Academic contribution: Provides new theoretical tools for bioconvection numerical methods
  2. Application prospects: Has potential applications in bioengineering, environmental science, and related fields
  3. Method generalization: CNLF framework may be applicable to other similar coupled systems

Applicable Scenarios

  1. Microorganism suspension modeling: Suitable for biological fluids requiring viscosity variation consideration
  2. Environmental fluid mechanics: Can be used to simulate natural water body flows containing microorganisms
  3. Bioreactor design: Provides numerical tools for bioreactor optimization design

Technical Details Supplement

Key Mathematical Tools

  1. Skew-symmetric trilinear form: B(u,v,w)=12Ω(uv)wdx12Ω(uw)vdxB(u,v,w) = \frac{1}{2}\int_\Omega (u \cdot \nabla v) \cdot w dx - \frac{1}{2}\int_\Omega (u \cdot \nabla w) \cdot v dx
  2. Variable coefficient projection operator: ν(c)((uPhn+1u),vh)+(vh,pρhn+1p)=0\nu(c)(\nabla(u-P_h^{n+1}u), \nabla v_h) + (\nabla \cdot v_h, p-\rho_h^{n+1}p) = 0
  3. Discrete Gronwall inequality: Key tool for stability analysis

Numerical Implementation Key Points

  1. Initial value treatment: First step employs backward Euler method to ensure accuracy
  2. Mass conservation: Preserves total microorganism mass through appropriate function space selection
  3. Linear system solving: Each time step requires solving only linear systems, improving efficiency

This paper makes important contributions in both theoretical and numerical aspects, providing an effective numerical method for handling complex biological fluid problems. Although there are some theoretical assumptions and dimensional limitations, its innovative CNLF scheme and rigorous analysis establish an important foundation for related research in the field.