Recently I bacame interested in Domain Decomposition algorithms for spectral methods, and in particular for Chebyshev collocation method.

There has been a lot of activity in this field since the second half of the 80’s. It will be hard to list all the papers, but maybe in some future post I will make a short list of some interesting ones.

Thanks to a funtastic scientific python blog I was able to get my hands on a masterfully written FD code implementing Multiplicative Schwarz method for Poisson equation on a square. With minor changes in the routine evaluating the 2D Laplacian, I was able to adapt the code for Chebyshev pseudospectral method. The results were wonderful!

I made two examples, and implemented them in another-chebpy. Both examples are in my google code git repository. One is solving Poisson eq. and other solves Biharmonic BVP (clamped BCs). I uploaded a small video on youtube that shows (rather slow) convergence for the case of biharmonic equation.

This can be accelerated by multigrid, which will be done in the future. Also parallelization is possible, but in this cae multicoloring (available here already) is required to decouple domains.

Following equations will help us in developing a collocation method using Bernstein polynomials for two-dimensional problems.

Let the following equation represent a surface in $\mathbb{R}^3$ defined by a tensor product of Bernstein basis functions:

$f(x,y)=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}B_{i,n}(x)B_{j,m}(y)$

We can define p-th order partial derivative in $x$ -axis direction in a following way

$\frac{\partial^pf(x,y)}{\partial x^p}=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}\frac{\partial^p B_{i,n}(x)}{\partial x^p}B_{j,m}(y)$

and an q-th order derivative in $y$ -axis direction

$\frac{\partial^qf(x,y)}{\partial y^q}=\sum_{i=0}^n \sum_{j=0}^m\beta_{i,j}B_{i,n}(x)\frac{\partial^qB_{j,m}(y)}{\partial y^q}$

Finally mixed derivative of order $p+q$ is defined as follows:

$\frac{\partial^{p+q}f(x,y)}{\partial x^p y^q}=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}\frac{\partial^pB_{i,n}(x)}{\partial x^p}\frac{\partial^qB_{j,m}(y)}{\partial y^q}$.

Laplacian operator can be written like this

$\Delta f(x,y)=\frac{\partial^2 f(x,y)}{\partial x^2 }+\frac{\partial^2 f(x,y)}{\partial y^2}$

$=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}B_{j,m}(y)\frac{\partial^2B_{i,n}(x)}{\partial x^2}+\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}B_{i,n}(x)\frac{\partial^2 B_{j,m}(y)}{\partial y^2}$,

which can further be simplified to

$\Delta f(x,y)=\sum_{i=0}^n \sum_{j=0}^m\beta_{i,j}\left[B_{j,m}(y)\frac{\partial^2 B_{i,n}(x)}{\partial x^2}+B_{i,n}(x)\frac{\partial^2 B_{j,m}(y)}{\partial y^2}\right]$.

“I want you, not to be able to make residual with the same basis functions, you used to create the solution.”

(P.S. This story is completely not true, and I urge my readers to find a better explanation of (Bubnov-) Galerkin method, if it exists.)

FLUENT users are aware of the “realizable” $k-\epsilon$ turbulence model. For more experienced modelers this means an implementation of a model presented in the paper by T.H. Shih, W.W. Liou,  A. Shabbir, Z. Yang and J. Zhu, “A New $k-\epsilon$ Eddy Viscosity Model for High Reynolds Number Turbulent Flows”, Computers and Fluids, Vol. 24 No. 3, pp. 227-238, 1995.

In this model, realizability is achieved by letting  $C_\mu$ coefficient appearing in the formulation of eddy viscosity, to be no longer a constant. Experiments also suggest it shouldn’t be a constant. While it has a value of 0.09 in inertial sublayer of BL flow, in homogeneous shear flow it is around 0,05. Therefore some way of changing the $C_\mu$ coefficient depending on the flow situation is desirable.

Let’s first say what realizability means mathematically, and then we’ll go into details about  $C_\mu$.

The realizability constraints any model should meet are following inequalities regarding Reynolds stress tensor components:

$\overline{u_\alpha^2}\textgreater0,\,(\alpha = 1,2,3)$

$\frac{\overline{u_\alpha u_\beta}^2}{\overline{u_\alpha^2} \, \overline{u_\beta ^2}} \leq 1, (\alpha = 1,2,3; \beta = 1,2,3)$

Eq. (2) is know as Schwarz’s inequality. So we need positivity of normal stresses (Eq.(1)), and we need Reynolds stress tensor components to obey Schwartz’s inequality.  Based on these realizability constraints W. Reynolds as well as Shih, Zhu and Lumley  proposed following form of an equation for $C_\mu$:

$C_\mu=\frac{1}{A_0+A_sU^{(*)}\frac{k}{\epsilon}}$

Shih, Zhu, Lumley give this formulation for the constituents of the previous equation

$U^{(*)}=\sqrt{S_{ij}S_{ij}+\tilde{\Omega}_{ij}\tilde{\Omega}_{ij}}$,

Where

$\tilde{\Omega}_{ij}=\Omega_{ij}-2\epsilon_{ijk}\omega_k$

$\Omega_{ij}=\overline{\Omega_{ij}}-\epsilon_{ijk}\omega_k$

where $\overline{\Omega_{ij}}$ is the mean rate-of-rotation tensor viewed in a rotating reference frame with the angular velocity $\omega_k$. So this model can be applied in cases like rotating pipe etc.

$A_s=\sqrt{6}cos\phi$, $\phi=\frac{1}{3} arccos(\sqrt{6}W)$, $W=\frac{S_{ij}S_{jk}S_{ki}}{\tilde{S}^3}$, $\tilde{S}=\sqrt{S_{ij}S_{ij}}$.

Only undetermined is coefficient $A_0$. As far as I’m aware all authors treat it as a constant, and calibrate it in a way that it returns desired values in some cases (like 0.09 in inertial sublayer in BL flow). In original reference authors came to the value of $A_0=4.0$, while in FLUENT this model is implemented with value $A_0=4.04$.  I look forward to some study where $A_0$ might not be a constant.

Now going to transport equations for $k$ and $\epsilon$ . The first one ($k$-equation) remains the same as in standard $k-\epsilon$ model, while there are some changes with equation for dissipation. Here I’ll show them both for the sake of completeness

$\frac{\partial}{\partial t} (\rho k) + \frac{\partial}{\partial x_j} (\rho k u_j) =$

$\frac{\partial}{\partial x_j} \left [ \left(\mu + \frac{\mu_t}{\sigma_k}\right) \frac{\partial k} {\partial x_j} \right ] + P_k + P_b - \rho \epsilon - Y_M + S_k$

$\frac{\partial}{\partial t} (\rho \epsilon) + \frac{\partial}{\partial x_j} (\rho \epsilon u_j) =$

$\frac{\partial}{\partial x_j} \left[ \left(\mu + \frac{\mu_t}{\sigma_{\epsilon}}\right) \frac{\partial \epsilon}{\partial x_j} \right ] + \rho \, C_1 S \epsilon - \rho \, C_2 \frac{{\epsilon}^2} {k + \sqrt{\nu \epsilon}} + C_{1 \epsilon}\frac{\epsilon}{k} C_{3 \epsilon} P_b + S_{\epsilon}$

Where

$C_1 = \max\left[0.43, \frac{\eta}{\eta + 5}\right] , \;\;\;\;\; \eta = S \frac{k}{\epsilon}, \;\;\;\;\; S =\sqrt{2 S_{ij} S_{ij}}$

Production of turbulence kinetic energy due to velocity gradients $P_k$, as well as production of TKE due to buoyancy $P_b$, are treated in the same way as in standard $k-\epsilon$ model.

When we find both $k$ and $\epsilon$ fields we move to computation of eddy viscosity in the standard form for $k-\epsilon$ model

$\mu_t=\rho C_\mu \frac{k^2}{\epsilon}$

The purpose of this note is to point out a few things I found useful for stable implementation of the model.

Argument of the arccos function should be in interval [-1,1], therefore I suggest implementation of $\phi$ in a following way:

$\phi=\frac{1}{3}, \;\;arccos(max(-1, min(\sqrt{6}W, 1)))$

Otherwise you might get  a few NaN’s in the $C_\mu$ field which usually ruin the solution.

Experiments suggest there is no reason to let $C_\mu$ be larger than 0.09, so we will limit it from above

$C_\mu =min(C_\mu, 0.09)$

It would be useful to see what are the lowest values in the $C\mu$ field and check with some experimental data, just to make sure we are on the right track.

Finally constants of the model are:

$C_{1 \epsilon} = 1.44, \;\; C_2 = 1.9, \;\; \sigma_k = 1.0, \;\; \sigma_{\epsilon} = 1.2$

I wish you happy computing!

I was interested to see how mlab handles plots from my Chebyshev pseudospectral code. The result was quite surprising. The light effects give nice impression of this object’s position in space. The colors you might chose from are the same as in Matplotlib, but here they change smoothly – the one thing I was missing there.
I also posted this figure to figshare here.