High-order stabilized matrix-free simulation of rotating mixing devices using the Mortar Element Method

B. Campos; P. Munch; V. O. Ferreira; D. C. Boffito; X. Banquy; B. Blais

Source-linked AI summary

High-order stabilized matrix-free simulation of rotating mixing devices using the Mortar Element Method

B. Campos, P. Munch, V. O. Ferreira, D. C. Boffito, X. Banquy, B. Blais

arXiv:2608.27423v1Aug 27, 2026physics.flu-dynmath.NA

TL;DR

Rotating mixing-device simulations require accurate treatment of moving, nonconforming rotor–stator interfaces and efficient computation at high order. The paper combines MEM with a stabilized, matrix-free finite-element Navier–Stokes framework and reports optimal convergence, strong scalability, and agreement with numerical and experimental mixer data. The reported numerical dissipation is ∼1% at Re = 200 and 10% at Re = 2000.

Problem
Rotating mixing devices involve complex geometries, turbulent flows, mesh deformation, and costly mortar-cell and boundary-integral computations.
Method
The model combines MEM with a high-order matrix-free finite-element Navier–Stokes solver using ALE motion, SUPG–PSPG stabilization, and DG-style interface terms.
Results
Two-dimensional cases recover optimal convergence without significant interface error, while three-dimensional and mixer cases show strong scalability and agreement with numerical and experimental results.
Takeaways & Limitations
The framework supports accurate and efficient simulations of rotating geometries across verification, scalability, and mixer-validation cases.
Takeaways & Limitations
The current model restricts the interface to equal-sized cells, limiting mesh flexibility for more complex geometries.

Abstract

from arXiv · show

We present a finite element framework to simulate rotating mixing devices using the Mortar Element Method as a domain decomposition strategy. The model is implemented within a matrix-free Navier-Stokes framework which uses a high-order Continuous Galerkin method. The discretized domain is subdivided into rotor and stator parts. An Arbitrary Lagrangian-Eulerian approach accounts for the relative rotor-stator motion, and stabilization is ensured through the Streamline-Upwind/Petrov-Galerkin and Pressure-Stabilizing Petrov-Galerkin methods. The rotor-stator domains are connected by an interface composed of mortar cells, and continuity is weakly enforced in a Discontinuous Galerkin fashion by accounting for boundary integrals at the rotor-stator interface. Verifications of the convergence order in two-dimensional steady and transient examples report optimal rates. The geometric non-conformity created at the mortar interface due to rotor rotation does not introduce significant error in the solution. A three-dimensional example is used to investigate the model's scalability, which yields ideal strong scaling for large problems. A two-dimensional Rushton impeller example uses a torque analysis to showcase the mesh convergence, and the corresponding velocity profile is in agreement with existing numerical results. In a three-dimensional pitched blade turbine case, the power number curve (Np vs Re) shows good agreement with experimental data for Reynolds number values from 1 to 2000. An energy balance analysis reports a numerical dissipation of 1% for Re=200 and of 10% for Re=2000. By exploiting modern hardware capabilities through matrix-free methods, the proposed model is a robust, accurate, and efficient framework suitable for simulating flows with rotating geometries.

1. Introduction

The paper motivates a mortar-based framework for rotating-device flows, combining weakly coupled rotor–stator domains with high-order, stabilized, matrix-free finite elements. This addresses interface complexity and computational cost while retaining accuracy and parallel efficiency.

Motivation: Rotating mixing devices are difficult to investigate experimentally because they combine complex geometries with turbulent flows.Large mesh deformation from rotation can also require remeshing.
Mortar Element Method: The Mortar Element Method subdivides the domain into non-overlapping subdomains and weakly enforces continuity across mortar elements.For rotating devices, the subdomains represent rotor and stator regions.
Computational challenge: MEM requires extra numerical effort to create mortar cells and compute additional boundary integrals.These costs make efficient parallelization important for large simulations.
Computational strategy: Matrix-free algorithms reduce memory-bandwidth limitations by exploiting floating-point operations, while high-order approximations improve accuracy on coarser meshes.The implementation uses on-the-fly vector–matrix multiplication in batches for modern HPC architectures.
Proposed framework: The proposed model combines MEM with a matrix-free finite-element solver using high-order continuous Galerkin Navier–Stokes discretization, ALE motion, and SUPG–PSPG stabilization.Interface continuity is imposed with DG-style boundary terms and SIPG treatment.

2. Governing equations

The governing formulation uses incompressible Navier–Stokes equations in an ALE frame, discretized by finite elements over rotor–stator subdomains. Mortar-interface continuity is imposed weakly with DG boundary terms, while SUPG and PSPG stabilize the formulation.

ALE Navier–Stokes equations: The ALE momentum equation includes fluid advection, mesh-velocity advection, pressure, viscous diffusion, and forcing.The incompressibility constraint is ∇·u = 0.
Finite-element discretization: Finite-element discretization applies test functions to the incompressibility and momentum equations and assembles their weighted-residual form element by element.The resulting nonlinear system is later linearized with Newton’s method.
Mortar decomposition: MEM partitions the domain into N non-overlapping subdomains; the rotor–stator configuration uses N = 2.Mortar cells are created from intersections between boundary-cell edges and the interface.
Interface coupling: Velocity and pressure continuity across the rotor–stator interface is enforced through DG-style boundary integrals using central pressure fluxes and SIPG viscous terms.Average and jump operators describe traces from the two interface sides.
Stabilization: SUPG and PSPG terms stabilize the weak form using the strong momentum residual and a stabilization parameter based on timestep and element size.The strong residual includes the ALE momentum equation terms.

3. Software implementation

The implementation integrates curved mortar interfaces into Lethe’s matrix-free finite-element infrastructure. It constructs interface operators, maps quadrature through rotating geometry, exchanges interface data, and supports geometric multigrid without assembling global matrices.

Software architecture: Lethe implements the mortar method using deal.II, p4est, Trilinos, and MPI, while supporting both matrix-based and matrix-free architectures.The reported simulations use the matrix-free implementation.
Mortar operators: The Mortar Manager stores mesh alignment, mortar-cell indices, quadrature weights, and normals for curved interface cells.The Coupling Operator handles interface-side data, penalties, and constraints; Coupling Evaluation computes mortar terms.
Rotor motion: Rotor-node positions are updated by applying a rotation matrix to the initial configuration.The three-dimensional rotation uses the axis vector, identity tensor, and skew-symmetric tensor.
Interface integration: Boundary integrals are approximated with Gauss quadrature after mapping a reference quadrature onto each mortar segment.The mapped Jacobian depends on interface radius, segment angle, and axial cell height.
Matrix-free evaluation: Matrix-free mortar integration evaluates interface variables, communicates data to ghost cells, and then integrates stored values.This organization supports repeated computations on cached data.
Geometric multigrid: Geometric multigrid coarsens mesh size and polynomial degree, while reconstructing mortar contributions without storing full matrices.The reconstruction is performed for diagonals and element contributions as needed by smoothers and coarse-grid solvers.

4. Numerical examples

The numerical studies verify convergence, examine scalability, and assess implementation assumptions for high-order rotor–stator simulations. The setup uses GMRES with geometric multigrid and inverse-diagonal smoothing, with BDF2 for transient cases.

Study design: The verification suite includes steady manufactured-solution and transient Taylor–Green vortex cases, followed by scalability and rotating-mixer studies.The mixer cases compare against existing numerical and experimental data.
Mortar reconstruction: The matrix-free implementation reconstructs mortar contributions for both the global matrix diagonal and the full matrix through local degree-of-freedom loops.The algorithms evaluate basis vectors at quadrature points before assembling the corresponding contributions.
Parallel coupling: The implementation exchanges evaluated interface values to ghost cells before integrating mortar boundary terms.This communication step is part of the matrix-free coupling evaluation.
Numerical setup: The linearized Navier–Stokes systems use GMRES, a GCMG preconditioner, and inverse-diagonal smoothing.Transient simulations use second-order implicit BDF2 time integration and equal-order QPQP elements with P ≥ 2.

4.1. Steady-state manufactured solution

The steady two-dimensional MMS verifies optimal high-order convergence and shows that the mortar interface matches the reference solution without adding error, even under rotor-mesh rotation.

Problem setup: The manufactured problem combines a plate with a hole and a circle, with a circular rotor-stator interface at radial distance 0.5.The domain is [−1, 1] × [−1, 1], and the rotor angular velocity is ω = 0.
Discretization: The study uses ν = 1, four refinement levels, and polynomial approximations from P = 2 to P = 5.The time-dependent term is neglected for this steady case.
Convergence: Velocity and pressure achieve the optimal convergence rates O(h^(P+1)) and O(h^P), respectively.Mesh size is characterized using average element area because the rotor-stator mesh is nonuniform.
Mortar accuracy: The subdivided rotor-stator domain and the no-mortar reference case produce identical L2 error norms.Both cases use the same parameters and mesh configuration apart from the mortar interface.
Rotating interface: Error norms remain constant across initial rotor rotations, confirming that mortar-cell misalignment does not cause loss of convergence properties.The tested rotation interval is [0, θmax], with θmax = π/8 on the coarsest mesh.

4.2. Taylor-Green vortex

The transient Taylor-Green vortex verifies the rotating-interface implementation against analytical behavior. Optimal convergence is retained across rotor angular velocities, with good numerical-analytical agreement.

Setup: The transient verification uses the two-dimensional Taylor-Green vortex with prescribed rotor angular velocity.The simulation uses ν = 1, ∆t = 10^-4 s, and tend = 0.1 s.
Analytical reference: The analytical velocity solution is u = exp(−2νt) cos(x) sin(y) and v = −exp(−2νt) sin(x) cos(y).These expressions provide the reference for the transient error analysis.
Viscid convergence: The Q3Q3 case achieves the expected velocity convergence rate O(h^(P+1)) for all tested rotor angular velocities.For P = 2, convergence deteriorates because of the odd-degree structure of the velocity Taylor expansion; Q5Q5 recovers the expected order.
Inviscid convergence: In the inviscid case, differences between rotating and stationary error norms remain on the order of 10^-4, with the same convergence trends for P = 2 and P = 3.The inviscid analysis assumes ν = 0 and Re →∞.
Solution comparison: Computed absolute errors range from 10^-4 to 10^-13 in the finest-mesh diagonal comparison with the analytical solution.Higher errors occur in the stator mesh, whose average elements are larger.

4.3. Performance and scalability tests

The three-dimensional tests assess solver performance, scalability, and mortar workload distribution. Matrix-free execution scales well for large problems, while mortar setup can become a substantial cost and source of imbalance.

Scalability setup: The three-dimensional MMS scalability study introduces a mortar interface by combining a cylinder with an extruded plate containing a hole.The study evaluates matrix-free MEM performance using multiple mesh configurations.
Solver setup: The steady-state tests use two coarse mesh configurations and polynomial degrees P = 2 and P = 3, with hp-multigrid used for Q3Q3.The coarsest hp-multigrid polynomial degree is reduced to P = 2.
Solver convergence: The number of linear solver iterations is more affected by cell aspect ratio than by the presence of the mortar interface.The l0 = 40 discretization has fewer iterations overall because it contains more uniformly shaped cells.
Preconditioning: ASM requires about one third fewer iterations than ID but can take up to 3.4 times more computational time.This reflects ASM’s more expensive setup and application despite its lower iteration count.
Scaling results: Weak-scaling time increases by approximately 41% from 163.8k to 1.3M cells and by approximately 105% from 1.3M to 10.5M cells.Each comparison increases both mesh size and node count by a factor of 8; communication costs affect weak scaling.
Mortar overhead: Mortar setup consumes around 20–50% of total simulation time and its workload imbalance grows with the number of processes.Mortar setup includes mortar-operator creation, boundary-integral computation, and rotor-mapping updates.
Load balancing: Using wm/wr = 20 reduces workload imbalance by up to five times with only small simulation-time variations.The load-balancing step itself takes slightly more time and is deemed unnecessary for the present mortar approach.

4.4. Two-dimensional Rushton mixer

The two-dimensional Rushton mixer demonstrates agreement with existing velocity results and mesh-convergent torque measurements. Refinement drives the boundary-torque sum toward zero in the laminar regime.

Numerical setup: The simulation runs for twenty revolutions to ensure the flow is fully established at tend = 20 s.Adaptive refinement is enabled while keeping the number of interface elements fixed.
Unbaffled case: The unbaffled MEM velocity profile shows good agreement with the published finite-volume solution.The comparison is made along the x-axis at t = 1 s.
Torque analysis: Torque measurements at the impeller and across all boundaries converge under four uniform refinement levels.The sum of boundary torques converges to zero, as expected for the laminar regime.
Mesh convergence: The outer-wall torque variation decreases from approximately 2.7% between l = 0 and l = 1 to 0.3% between l = 2 and l = 3.For biquadratic velocity and pressure fields, a second-order torque convergence rate is expected.

4.5. Three-dimensional pitched-blade turbine (PBT) mixer

The three-dimensional PBT mixer study evaluates mesh convergence, energy balance, and validation across Reynolds numbers using a matrix-free mortar-element solver. Results show converged torque, low numerical dissipation at Re = 200, and agreement with experimental power-number data, while higher-Re energy residuals remain more mesh sensitive.

Setup: The PBT validation varies fluid viscosity to impose Reynolds numbers over the experimental interval, using a fixed impeller speed of N = 16 Hz.The setup applies no-slip conditions to the outer cylinder, bottom wall, and baffles, with slip at the top boundary.
Mesh sensitivity: A mortar interface with diameter Dint = 1.2Di uses uniform cell sizes in the x-y plane and z-direction.The interface restriction limits further coarsening of the initial configuration.
Mesh sensitivity: 3.34% and 1.05% are the approximate torque variations from l = 0 to l = 1 for P = 2 and P = 3, respectively, at Re = 200.The l = 1 quadratic approximation yields converged torque values, motivating its use for the subsequent study.
Interpretation: Torque and energy-balance curves exhibit oscillatory pseudo-steady profiles caused by imposed rotation and blade-baffle interactions.The rotation-induced component asymptotically decreases with time and mesh refinement, while blade-baffle interaction produces an angle-dependent sinusoidal-like signal.
Mesh sensitivity: At Re = 2000, the coarse l = 0 Q2Q2 mesh diverges at t ≈ 0.05 s, while finer meshes produce closely agreeing torque curves.The higher-Re torque curves also exhibit more erratic oscillations due to increased turbulence.
Energy analysis: 1% numerical dissipation is reached on the finest Re = 200 meshes, compared with approximately 2% on the coarsest Q2Q2 mesh and ≈10% for Q1Q1 elements.At tend = 0.2 s, approximately 95% of the system’s kinetic energy is achieved.
Energy analysis: 10% and 3% numerical dissipation occur at Re = 2000 on l = 1 and l = 2 Q2Q2 meshes, while Q3Q3 at l = 1 yields around 5%.The coarsest mesh diverges, and the energy residual remains more mesh sensitive than the torque metric.
Validation: The power-number curve shows good agreement with experimental results, especially at higher Reynolds numbers, and is compared with prior Nitsche immersed-boundary results.The comparison includes setup and computational-time information for matrix-free MEM and matrix-based NIB configurations.

5. Conclusions and further research

The conclusions present a matrix-free high-order mortar-element framework for rotating mixers and summarize verification, scalability, validation, and energy-balance results. They identify equal-sized interface cells as a current limitation and propose broader geometry and multiphysics extensions.

Conclusions: The model combines a CG fluid discretization with a DG mortar interface, weak boundary-integral continuity, SIPG coupling, and ALE rotor-stator motion.The framework uses the Mortar Element Method within a matrix-free high-order formulation.
Verification and scalability: Two-dimensional steady and transient tests recover optimal velocity and pressure convergence rates despite rotation-induced interface misalignment.A three-dimensional steady test also reports good strong scalability for sufficiently large problems.
Verification and scalability: Improving mortar-cell load balancing did not significantly reduce computing time.This result is reported within the three-dimensional scalability study.
Validation: The converged meshes yield approximately 1% numerical dissipation at Re = 200 and 10% at Re = 2000.The PBT power-number curve agrees well with experimental results over the reported Reynolds-number range.
Further research: Future work targets devices with more complex geometries or flow patterns and coupling additional physics such as particle transport.These extensions are intended to broaden the model’s engineering applications.
Limitations: Equal-sized cells are currently required at the mortar interface, limiting mesh flexibility for more complex geometries.Future extensions may allow unstructured meshes and varying interface radii.