Numerical method for the kinetic model governing the probability distribution function of liquid crystalline polymers

Abstract

We first present the kinetic model that governs the orientational probability distribution function for liquid crystalline polymers. Then we summarize the numerical method used to solve the model equations.

Keywords:

• Liquid crystalline polymers
• kinetic model
• numerical method

1. Kinetic model of liquid crystalline polymers

We consider viscous fluid suspensions of liquid crystalline polymers which are either macromolecular rods or platelets. Individual particles are advected by a superposition of a macroscopic velocity $\mathbf{v}$ of the complex fluid mixture and the rotational velocity given by the Jeffery orbit equation: (1) $\dot{\mathbf{m}}=\mathrm{\Omega }\cdot \mathbf{m}+a[\mathbf{D}\cdot \mathbf{m}-\mathbf{D}:\mathbf{m}\mathbf{m}\mathbf{m}],$(1) where $\mathbf{m}$ is the orientation of the individual molecular, a is the particle aspect ratio parameter, $\mathbf{D}$ and Ω are rate-of-strain and vorticity tensors, respectively. The standard Maier–Saupe or Onsager excluded volume potential [1, 2] (2) $U=-\frac{3N}{2}⟨\mathbf{mm}⟩:\mathbf{mm},$(2) gives rise to nematic equilibria in the semi-dilute limit, where N is the strength of the nematic order. The average $⟨\cdot ⟩$ is over the unit sphere in orientational space, (3) $⟨(\cdot )⟩={\int }_{||\mathbf{m}||=1}(\cdot )f(\mathbf{m},\mathbf{x},t)\mathrm{d}\mathbf{m},$(3) from which various moments of the probability distribution function (PDF) $f(\mathbf{x},\mathbf{m},t)$ can be defined.

The Smoluchowski or Fokker-Planck equation governing f is given by [3–8] (4) $\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }f=D_s\mathrm{\nabla }\cdot (\mathrm{\nabla f}+f\mathrm{\nabla U})+D_r\mathcal{R}\cdot (\mathcal{R}f+f\mathcal{R}U)-\mathcal{R}\cdot (\mathbf{m}\times \dot{\mathbf{m}}f),$(4) where $\mathcal{R}$ is the rotational gradient operator, $D_s$ and $D_r$ are translational and rotational diffusivity coefficients, respectively. Note that the Maier–Saupe potential introduces nonlinear integro-differential terms into the kinetic Equation (4(4) $\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }f=D_s\mathrm{\nabla }\cdot (\mathrm{\nabla f}+f\mathrm{\nabla U})+D_r\mathcal{R}\cdot (\mathcal{R}f+f\mathcal{R}U)-\mathcal{R}\cdot (\mathbf{m}\times \dot{\mathbf{m}}f),$(4) ) above. For the macroscopic velocity $\mathbf{v}$ and the density of the material ρ, we have the mass and momentum conservation, respectively, (5) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \mathbf{v})=0,\\ & \rho [\frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }\mathbf{v}]=\mathrm{\nabla }\cdot (-p\mathbf{I}+\tau ),\end{array}$(5) where p is the hydrostatic pressure, $\begin{array}{rl}\tau & =\frac{2}{\mathrm{Re}}\mathbf{D}+G[\mathbf{M}-\frac{\mathbf{I}}{3}-N{\mathbf{M}}^2+N\mathbf{M}:⟨\mathbf{mmmm}⟩]+\frac{1}{R{e}_2}[\mathbf{D}\cdot \mathbf{M}+\mathbf{M}\cdot \mathbf{D}]\\ & +\frac{1}{R{e}_3}⟨\mathbf{mmmm}⟩:\mathbf{D}\end{array}$ is the extra stress tensor including contributions from the viscous solvent and the elasticity of the molecular rod ensemble.

2. Numerical methods

To solve the above kinetic model for liquid crystalline polymers numerically, the probability distribution function $f(\mathbf{x},\mathbf{m},t)$ is first expanded using spherical harmonic functions (6) $f(\mathbf{x},\mathbf{m},t)\approx \sum _{l=0}^L\sum _{m=-l}^l{a}_l^m(\mathbf{x},t)Y_l^m(\mathbf{m}),$(6) where, (7) $Y_l^m(\mathbf{m})=\sqrt{\frac{(2l-1)(l-m)!}{4\pi (l+m)!}}{P}_l^m(\cos \theta )e^{i\mathrm{m\varphi }},\mathbf{m}=\left(\begin{array}{c}\sin \mathrm{\theta cos}\varphi \\ \sin \mathrm{\theta sin}\varphi \\ \cos \theta \end{array}\right)$(7) and $P_l^m$'s are associated with Legendre polynomials. Based on the spherical harmonic expansion of f, the potential (2(2) $U=-\frac{3N}{2}⟨\mathbf{mm}⟩:\mathbf{mm},$(2) ) has simple expansion (8) $U=-2\mathrm{\pi N}{a}_0^0{Y}_0^0-\frac{4\pi }{5}N\sum _{p=-2}^2{a}_2^p{Y}_2^p,$(8) involving only the zeroth and second moments of f. After complicated derivations, the Smoluchowski Equation (4(4) $\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }f=D_s\mathrm{\nabla }\cdot (\mathrm{\nabla f}+f\mathrm{\nabla U})+D_r\mathcal{R}\cdot (\mathcal{R}f+f\mathcal{R}U)-\mathcal{R}\cdot (\mathbf{m}\times \dot{\mathbf{m}}f),$(4) ) can be transformed to a large system of PDEs for the amplitudes $a_l^m(\mathbf{x},t)$ for all spherical harmonic modes. If we take the truncation in (6(6) $f(\mathbf{x},\mathbf{m},t)\approx \sum _{l=0}^L\sum _{m=-l}^l{a}_l^m(\mathbf{x},t)Y_l^m(\mathbf{m}),$(6) ) at L = 10, then 121 PDEs will be produced in the system. Coupled with Navier–Stokes Equation (5(5) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \mathbf{v})=0,\\ & \rho [\frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }\mathbf{v}]=\mathrm{\nabla }\cdot (-p\mathbf{I}+\tau ),\end{array}$(5) ), these PDEs can be discretized in space and integrated in time as shown below.

We choose two-dimensional physical x-y space in our simulations. The above system of PDEs is discretized in space using either the second- or fourth-order finite difference methods (for convenience, the PDF f, instead of the amplitude $a_l^m$, is used in the following). (9) $\begin{array}{rl}& \frac{\mathrm{d}{f}_{i,j}}{\mathrm{d}t}+\mathbf{v}\cdot {\mathrm{\nabla }}_{\mathbf{h}}{f}_{i,j}\\ & =D_s({\mathrm{\nabla }}_{\mathbf{h}}^2{f}_{i,j}+{\mathrm{\nabla }}_h\cdot (f_{i,j}{\mathrm{\nabla }}_{\mathbf{h}}{U}_{i,j}))+D_r\mathcal{R}\cdot (\mathcal{R}{f}_{i,j}+f_{i,j}\mathcal{R}{U}_{i,j})-\mathcal{R}\cdot (\mathbf{m}\times \dot{\mathbf{m}}{f}_{i,j}),\end{array}$(9) where $f_{i,j}$ is the numerical value at the grid point $(x_i,y_j)$ and ${\mathrm{\nabla }}_{\mathbf{h}}$ is the standard finite difference operator. In time, the 1st order semi-implicit method can be used to integrate the system, and then the spectral deferred correction method [9, 10] can be used to increase the order in time. For some detail, we first split the right side (9(9) $\begin{array}{rl}& \frac{\mathrm{d}{f}_{i,j}}{\mathrm{d}t}+\mathbf{v}\cdot {\mathrm{\nabla }}_{\mathbf{h}}{f}_{i,j}\\ & =D_s({\mathrm{\nabla }}_{\mathbf{h}}^2{f}_{i,j}+{\mathrm{\nabla }}_h\cdot (f_{i,j}{\mathrm{\nabla }}_{\mathbf{h}}{U}_{i,j}))+D_r\mathcal{R}\cdot (\mathcal{R}{f}_{i,j}+f_{i,j}\mathcal{R}{U}_{i,j})-\mathcal{R}\cdot (\mathbf{m}\times \dot{\mathbf{m}}{f}_{i,j}),\end{array}$(9) ) as a sum of an explicit and an implicit term: (10) $\frac{\mathrm{d}{f}_{i,j}}{\mathrm{d}t}=F_I(f_{i,j})+F_E(f_{i,j}),$(10) where $F_I$ and $F_E$ are two operators corresponding to the parts in the equation treated implicitly and explicitly, respectively (11) $F_I(f_{i,j})=D_s\left({\mathrm{\nabla }}_{\mathbf{h}}^2{f}_{i,j}\right)+D_r\mathcal{R}\cdot (\mathcal{R}{f}_{i,j}),F_E(f_{i,j})=L_E(f_{i,j})+N_E(f_{i,j}),$(11) where $L_E$ and $N_E$ are two operators corresponding to linear and nonlinear terms, respectively (12) $\begin{array}{rl}{L}_E(f_{i,j})& =-(\mathbf{v}\cdot {\mathrm{\nabla }}_{\mathbf{h}}{f}_{i,j}-\mathcal{R}\cdot (\mathbf{m}\times \dot{\mathbf{m}}{f}_{i,j}),\\ N_E(f_{i,j})& =D_s{\mathrm{\nabla }}_{\mathbf{h}}\cdot (f_{i,j}{\mathrm{\nabla }}_{\mathbf{h}}{U}_{i,j})+D_r\mathcal{R}\cdot (f_{i,j}\mathcal{R}{U}_{i,j}).\end{array}$(12) In the spectral deferred correction method, we need a start method that produces a lower order numerical solutions, and a correction method that is used to increase the time integration order. For the start method, the semi-implicit method (13) $\frac{f_{i,j}^{n+1}-f_{i,j}^n}{\mathrm{\Delta }t}=F_I(f_{i,j}^{n+1})+F_E(f_{i,j}^n)$(13) can be used. When transformed to spherical harmonic modes, it can be written as a large system of discretized nonlinear Helmholtz equations (14) $\begin{array}{rl}& [(1+\mathrm{\Delta }t{D}_{r}l(l+1))\mathbf{I}-\mathrm{\Delta t}{D}_s{\mathrm{\nabla }}_{\mathbf{h}}^2](a_l^m{)}_{i,j}^{n+1}=(a_l^m{)}_{i,j}^n+\mathrm{\Delta t}{F}_E(f_{i,j}^n),\\ & i=1,2,\dots N_x,j=1,2,\dots N_y,l=0,1,\dots L,m=-l,-l+1,\dots ,l-1,l.\end{array}$(14) We remark that, in the approach above (14(14) $\begin{array}{rl}& [(1+\mathrm{\Delta }t{D}_{r}l(l+1))\mathbf{I}-\mathrm{\Delta t}{D}_s{\mathrm{\nabla }}_{\mathbf{h}}^2](a_l^m{)}_{i,j}^{n+1}=(a_l^m{)}_{i,j}^n+\mathrm{\Delta t}{F}_E(f_{i,j}^n),\\ & i=1,2,\dots N_x,j=1,2,\dots N_y,l=0,1,\dots L,m=-l,-l+1,\dots ,l-1,l.\end{array}$(14) ), a total $(L+1{)}^2{N}_x{N}_y$ of Helmholtz equations need to be solved to get the 1st order numerical solutions. To overcome the computational cost, we can use another splitting of the equation for the correction method, namely, (15) $\frac{\mathrm{d}{f}_{i,j}}{\mathrm{d}t}={\tilde{F}}_I(f_{i,j})+{\tilde{F}}_E(f_{i,j}),\text{where},{\tilde{F}}_I(f_{i,j})=D_r\mathcal{R}\cdot (\mathcal{R}{f}_{i,j}),$(15) and all other terms are included in the explicit operator ${\tilde{F}}_E(f_{i,j})$. Then the correction method (16) $\frac{f_{i,j}^{n+1}-f_{i,j}^n}{\mathrm{\Delta }t}={\tilde{F}}_I(f_{i,j}^{n+1})+{\tilde{F}}_E(f_{i,j}^n),$(16) when transformed to spherical harmonic modes, becomes (17) $[1+\mathrm{\Delta }t{D}_{r}l(l+1)](a_l^m{)}_{i,j}^{n+1}=(a_l^m{)}_{i,j}^n+\mathrm{\Delta t}{\tilde{F}}_E(f_{i,j}^n)$(17) which is much easier to solve.

To calculate the explicit terms $F_E(f_{i,j})$ or ${\tilde{F}}_E(f_{i,j})$, the following formula involving multiplications of two spherical harmonic functions and the dot product of their rotational gradients are derived as follows: (18) $\begin{array}{rl}{\text{Y}}_j^p{\text{Y}}_l^m& =(-1{)}^{m+p}\sum _{k=|l-j|}^{l+j}{\alpha }_{\mathrm{lmjpk}}{\text{Y}}_k^{m+p}\\ \mathcal{R}{\text{Y}}_j^p\cdot \mathcal{R}{\text{Y}}_l^m& =(-1{)}^{m+p}\sum _{k=|l-j|}^{l+j}{\beta }_{\mathrm{lmjpk}}{\text{Y}}_k^{m+p}\\ {\alpha }_{\mathrm{lmjpk}}& ={\int }_{||\mathbf{m}||=1}{\text{Y}}_j^p{\text{Y}}_l^m{\text{Y}}_k^{-m-p}\mathrm{d}\mathbf{m}\\ {\beta }_{\mathrm{lmjpk}}& =\frac{1}{2}(j^2+j+l^2+l-k^2-k){\alpha }_{\mathrm{lmjpk}}\end{array}\}$(18) For the incompressible Navier-Stokes equation, we utilize the velocity-pressure formulation [11]. (19) $\begin{array}{rl}\frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }t}+(\mathbf{v}\cdot \mathrm{\nabla })\mathbf{v}+\mathrm{\nabla }p& =\eta \mathrm{\Delta }\mathbf{v}+\mathrm{\nabla }\cdot \tau ,\end{array}$(19) (20) $\begin{array}{rl}\mathrm{\Delta }p+\mathrm{\nabla u}\cdot \frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }x}+\mathrm{\nabla v}\cdot \frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }y}& =\alpha (\mathbf{x})\mathrm{\nabla }\cdot \mathbf{v}+\mathrm{\nabla }\cdot (\mathrm{\nabla }\cdot \tau ).\end{array}$(20) We remark that the artificial divergence damping term $\alpha (x)\mathrm{\nabla }\cdot \mathbf{v}$ (with $\alpha (\mathbf{x})\ge 0$) is added in the above equation to stabilize the algorithm.

To integrate the velocity equation, we again use semi-implicit method (21) $\frac{{\mathbf{v}}^{n+1}-{\mathbf{v}}^n}{\mathrm{\Delta }t}+({\mathbf{v}}^n\cdot {\mathrm{\nabla }}_{\mathbf{h}}){\mathbf{v}}^n+{\mathrm{\nabla }}_{\mathbf{h}}{p}^n=\eta {\mathrm{\nabla }}_{\mathbf{h}}^2{\mathbf{v}}^{n+1}+{\mathrm{\nabla }}_{\mathbf{h}}\cdot {\tau }^n.$(21) Then the spectral deferred correction method can be used to increase the integration order.

As a summary, in each time step of the whole algorithm, we first solve the elliptic Equation (20(20) $\begin{array}{rl}\mathrm{\Delta }p+\mathrm{\nabla u}\cdot \frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }x}+\mathrm{\nabla v}\cdot \frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }y}& =\alpha (\mathbf{x})\mathrm{\nabla }\cdot \mathbf{v}+\mathrm{\nabla }\cdot (\mathrm{\nabla }\cdot \tau ).\end{array}$(20) ) for the pressure p. Then spectral deferred correction method is used to solve the velocity field $\mathbf{v}$ through (21(21) $\frac{{\mathbf{v}}^{n+1}-{\mathbf{v}}^n}{\mathrm{\Delta }t}+({\mathbf{v}}^n\cdot {\mathrm{\nabla }}_{\mathbf{h}}){\mathbf{v}}^n+{\mathrm{\nabla }}_{\mathbf{h}}{p}^n=\eta {\mathrm{\nabla }}_{\mathbf{h}}^2{\mathbf{v}}^{n+1}+{\mathrm{\nabla }}_{\mathbf{h}}\cdot {\tau }^n.$(21) ). Finally, the Smoluchowski equation for the probability distribution function f is solved as described above.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Notes on contributors

Ruhai Zhou

Dr. Ruhai Zhou obtained his Ph.D. degree from Department of Mathematics, University of New Mexico in 2001. He did the postdoc work in Department of Mathematics, University of North Carolina at Chapel Hill from 2001 to 2004. Since 2004, he has been an Assistant Professor, Associate Professor, and Professor at Old Dominion University.