where $u$ and $v$ are the components of the local fluid velocity in the $x$ and $y$ directions, $\tau_{xy}$ is the local "polymeric" shear stress,
and $\tau_{xx}$ and $\tau_{yy}$ are local, polymeric normal stresses.
The Reynolds number $Re$, Weissenberg number $Wi$, and viscosity ratio $\beta$ are defined according to
$$Re = \frac{\rho U L}{\eta},$$
$$Wi = \lambda \frac{U}{L},$$
and
$$ \beta = \frac{\eta_s}{\eta}. $$
Here, $\rho$, $\eta$, $\eta_s$, $U$, and $L$ are the over all fluid density, fluid viscosity, solvent viscosity, lid speed, and lid length, respectively. The time scale $\bar t$
for this problem is defined by the characteristic shear rate according to $\bar t = \frac{L}{U}$ and the inertial pressure scale was chosen,
equal to $\bar p = \rho U^2$. The boundary conditions are no-slip and no-flow at the cavity walls and the fluid is initially at rest.
Numerical Scheme:
Following Lee and Leveque (2003),1 a fractional step approach is used to solve the system of equations above, which decomposes the problem into the following steps:
Step 1. Solve the convection equation:
A Roe solver $\textemdash$ i.e., a locally linear, approximate Riemann solver based on the Godunov method $\textemdash$
is employed to evaulate $q = (u,v,\tau_{xx},\tau_{xy},\tau_{yy})$ satisfying:
$$ q_t + \hat A \cdot q_x + \hat B \cdot q_y = 0 . $$
where the Roe matrices $\hat A$ and $\hat B$ are approximate Jacobian matrices given by
Here, $\hat u$ and $\hat v$ are Roe averages (in this case, they are linear interpolations of cell-centered velocities) defined at the edge of each grid cell. For example, the Roe averages used to evaluate the matrix $A$ are given by
Dimensional splitting via the donor cell upwind method (DCU) is used to advanced the cell-centered velocities and stresses $Q=(U,V,T_{xx},T_{xy},T_{yy})$ forward in time via sweeps in the x-direction
where $F_{i-\frac{1}{2},j}$ is the numerical flux at the interface between cells $(i,j)$ and $(i-1,j)$ for the 1-dimensional problem in the x-direction and, similarly, $G_{i,j-\frac{1}{2}}$ is the flux at the interface between cells $(i,j)$ and $(i,j-1)$ for the 1D problem in the y-direction. In addition, monotenzied central flux limiters are used to achieve second order accuracy for this step where the solution is smooth.
Step 2. Solve the diffusion equation:
During this step, the following equation is solved
An alternating direction implicit (ADI) method is employed to update the velocities for diffusion.
Two difference equations are used to advance successive time steps of $\frac{\Delta t}{2}$. The first equation for the u-velocity is implict in the x-direction
and $\delta^2_x$ denotes the central difference of the 2nd partial derivative with respect to $x$. Similar equations are used to update the v-velocity.
Step 3. Update the edge velocities for diffusion:
The velocities at the edges of each grid cell are updated via linear interpolation according to
So far, the velocity field $(\widetilde u, \widetilde v)$ is not divergence free. In order to satisfy continuity, $(\widetilde u, \widetilde v)$ is projected into a divergence-free vector field by correcting the result for pressure-driven flow via
which is then discretized with Neumann boundary conditions to produce a system of linear equations, $Ax = b$.
However, the matrix $A$ is singular because the equation set has an inifinite number of solutions within an arbitrary reference pressure.
Hence, a ficticious source term $C_0 p^{n+1}$ has been added with proportionality constant $C_0$,
which is defined on the order of 1e-9 to render the influence of the source negligble, so that $A$ is non-singular.
Step 5. Update the edge velocities for pressure-driven flow:
The edge velocities are determined by using central differences for the pressure gradient via
Velocity and pressure distribution for $Re = 0.1$, $Wi = 0.1$, and $\beta = 0$ with $CFL \leq 0.5$ using the DCU method on a 51x51 cell grid over the time interval $[0, 1]$
Vector_field_DCU_Re_0_1_Wi_0_1_beta_0.mp4
Streamlines and normal stress distribution for $Re = 0.1$, $Wi = 0.1$, and $\beta = 0$ with $CFL \leq 0.5$ using the DCU method on a 51x51 cell grid over the time interval $[0, 1]$
Streamlines_tau_xx_DCU_Re_0_1_Wi_0_1_beta_0.mp4
References:
L.Lee and R.J.LeVeque, 2003. An immersed interface method for incompressible Navier
Stokes equations. SIAM J. Sci. Comput., 25, 832–856.
Clawpack Development Team (2023), Clawpack Version 5.9.2,
http://www.clawpack.org, doi: 10.5281/zenodo.10076317
R. J. LeVeque, 1997. Wave propagation algorithms for multi-dimensional
hyperbolic systems. J. Comput. Phys. 131, 327–353.
R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge
University Press, Cambridge, UK, 2002.