2D LC Director Modeling Code
Vector Based

Version 1.0 July 09 1999

by
Salman Saeed

Whenever one wants to model the optics of a LC display one of the first things they have to do is obtain the configuration of the LC within the display. There are a number of programs available that will do this. Most are based on 1D code, i.e. the director is only allowed to fluctuate in one direction, typically around the z axis.

If however we allow the director to vary along, for example, the x axis also then we can model effects that fringing electric fields from pixeled electrodes have on the LC. The 2D vector based representation also does away with the stability problems of 1D theta-phi representations. In the theta-phi representation at high voltages the system becomes unstable as the director is upright in the cell and theta approaches 0. As many of the equations used have either the sine or cosine of theta as a divisor the system becomes unstable. Also at theta=0 the value of phi becomes ambiguous leading to more instability.

TECHNIQUE

As is the case in most director modelng programs we start off with the Frank-Oseen free energy equation. The difference here is that the director n is now represented in terms of a vector.

=Fs

The electric field however is where in 2 dimensions.

Of course n is also a vector and .
Now therefore we can write the electric energy part of the free energy as
Fe=.
is the dielectric tensor of the LC and is written as where j,k={x,z} in 2D.

The total free energy F is, and so the Euler-Lagrange equations become

where i={x,z}and Fg=Fs-Fe

Since the electric field depends on n, has to be solved for V after each time step. Values of V at discrete points on our grid are calculated using the method of finite differences.

Considering only the rotational viscosity and ignoring flow effects we can get the dynamics equation for the director.

The equation above can be solved by considering discrete time steps

Since the only purpose of the LaGrange multiplier is to maintain |n|=1 we can neglect l and calculate from in two steps.

and then normalize n

The algorithm to calculate the director configuration is shown in the flowchart below.

To summarize the method, we do the following

1. Set the initial director configuration based on the user input parameters.
2. Set the voltage through the system. Voltage is applied to one electrode and is linearly ramped down to zero at the second electrode.
3. Relax on the voltage using the initial director configuration until it has reached equilibrium. This leads to a faster relaxation process.
4. Calculate the director configuration given the new voltages at each grid point.
5. Check to see if the system has converged. If converged we end the program otherwise we repeat from step 4 and also calculate a new V given the new director configuration.

We shall use code generated by MAPLE as shown below to write our program.

Obtaining equations for use in the program via MAPLE

> # Maple Worksheet for 2D Vector representation of the Frank-Oseen free energy equation
> # Based on Work done by Hiroyuki Mori
> # Salman Saeed

> restart;

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> n:=vector([nx(x,z),ny(x,z),nz(x,z)]); Define n as a vector

> f1:=1/2*k11*(nxx+nzz)^2; First part of OF free energy eq after taking div of n.

> v:=[x,y,z];

> curl(n,v);

> f22:=1/2*k22*(dotprod(n,curl(n,v)))^2; Second part of OF free energy equation.

> f2:=1/2*k22*(-nx*nyz+ny*nxz-ny*nzx+nz*nyx)^2; Rewrite above for use later.

> a:=crossprod(n,curl(n,v));

> f31:=1/2*k33*dotprod(a, a); Third part of OF free energy equation

> f3:=1/2*k33*((ny*nyx-nz*nxz+nz*nzx)^2+(-nz*nyz-nx*nyx)^2+(nx*nxz-nx*nzx+ny*nyz)^2); Rewrite above

> q0*k22*(dotprod(n,curl(n,v))); Last part of OF free energy equation

> f4:=q0*k22*(-nx*nyz+ny*(nxz-nzx)+nz*nyx); Rewrite above

> Fs:=f1+f2+f3+f4; Finally the full OF free energy equation.

> E:=vector([Ex,Ey,Ez]); Define E as the electric field vector

> e:=linalg[matrix](3,3,[ep+de*nx*nx,de*nx*ny,de*nx*nz,de*ny*nx,ep+de*ny*ny,de*ny*nz,de*nz*nx,de*nz*ny,ep+de*nz*nz]); Define dielectric tensor

> Fe:=1/2*multiply(E,e0*multiply(e,E)); The electric field energy density

> Fg:=Fs-Fe; Gibbs free energy density

> diff(Fg,nx); For use in the Euler Lagrange equations we need Fg,nx Fg,ny, Fg,nz, Fg,nxx, Fg,nxz,Fg,nyx,Fg,nyz,Fg,nzx,Fg,nzz

> c1:=-k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*diff(ny(x,z),z)+1/2*k33*(-2*(-nz(x,z)*diff(ny(x,z),z)-nx(x,z)*diff(ny(x,z),x))*diff(ny(x,z),x)+2*(nx(x,z)*diff(nx(x,z),z)-nx(x,z)*diff(nz(x,z),x)+ny(x,z)*diff(ny(x,z),z))*(diff(nx(x,z),z)-diff(nz(x,z),x)))-q0*k22*diff(ny(x,z),z)-1/2*Ex*e0*(2*de*nx(x,z)*Ex+de*ny(x,z)*Ey+de*nz(x,z)*Ez)-1/2*Ey*e0*de*ny(x,z)*Ex-1/2*Ez*e0*de*nz(x,z)*Ex;

> diff(Fg,nxx);

> a1:=k11*(diff(nx(x,z),x)+diff(nz(x,z),z));

> diff(a1,x);

> diff(Fg,nxz);

> b1:=k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*ny(x,z)+1/2*k33*(-2*(ny(x,z)*diff(ny(x,z),x)-nz(x,z)*diff(nx(x,z),z)+nz(x,z)*diff(nz(x,z),x))*nz(x,z)+2*(nx(x,z)*diff(nx(x,z),z)-nx(x,z)*diff(nz(x,z),x)+ny(x,z)*diff(ny(x,z),z))*nx(x,z))+q0*k22*ny(x,z);

> EL_nx:=c1-diff(a1,x)-diff(b1,z); Euler lagrange Equation ELnx

> EL_nx_C:=simplify(%);

>

>

> readlib(C);

> C(EL_nx_C,ansi); Convert to C to be used in main program

C/expression:   Unknown function:   ny   will be left as is.

C/expression:   Unknown function:   diff   will be left as is.

C/expression:   Unknown function:   nx   will be left as is.

      s1 = -k22*pow(ny(x,z),2.0)*diff(diff(nx(x,z),z),z)+2.0*k22*pow(diff(ny(x,

z),z),2.0)*nx(x,z)+k33*nx(x,z)*pow(diff(ny(x,z),x),2.0)-k33*nx(x,z)*pow(diff(nx

C/expression:   Unknown function:   nz   will be left as is.

(x,z),z),2.0)+k22*pow(ny(x,z),2.0)*diff(diff(nz(x,z),x),z)+k33*nx(x,z)*pow(diff

(nz(x,z),x),2.0)-k33*pow(nz(x,z),2.0)*diff(diff(nx(x,z),z),z)+k33*pow(nz(x,z),

2.0)*diff(diff(nz(x,z),x),z)-k33*pow(nx(x,z),2.0)*diff(diff(nx(x,z),z),z)+k33*

pow(nx(x,z),2.0)*diff(diff(nz(x,z),x),z)-k33*nx(x,z)*pow(diff(ny(x,z),z),2.0)+

k22*ny(x,z)*nx(x,z)*diff(diff(ny(x,z),z),z)-k33*nx(x,z)*ny(x,z)*diff(diff(ny(x,

z),z),z)+k33*diff(nz(x,z),z)*ny(x,z)*diff(ny(x,z),x)+2.0*k33*nz(x,z)*diff(nz(x,

z),z)*diff(nz(x,z),x);

      t0 = s1-2.0*k33*nz(x,z)*diff(nz(x,z),z)*diff(nx(x,z),z)+k33*nz(x,z)*ny(x,

z)*diff(diff(ny(x,z),x),z)-k22*ny(x,z)*nz(x,z)*diff(diff(ny(x,z),x),z)-k22*ny(x

,z)*diff(nz(x,z),z)*diff(ny(x,z),x)-2.0*q0*k22*diff(ny(x,z),z)-Ey*e0*de*ny(x,z)

*Ex-Ez*e0*de*nz(x,z)*Ex-2.0*k22*diff(ny(x,z),z)*ny(x,z)*diff(nx(x,z),z)+3.0*k22

*diff(ny(x,z),z)*ny(x,z)*diff(nz(x,z),x)-2.0*k22*diff(ny(x,z),z)*nz(x,z)*diff(

ny(x,z),x)-k11*diff(diff(nx(x,z),x),x)-k11*diff(diff(nz(x,z),x),z)+2.0*k33*diff

(ny(x,z),x)*nz(x,z)*diff(ny(x,z),z)-k33*ny(x,z)*diff(ny(x,z),z)*diff(nz(x,z),x)

-Ex*Ex*e0*de*nx(x,z);

> diff(Fg,ny);

> c2:=k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*(diff(nx(x,z),z)-diff(nz(x,z),x))+1/2*k33*(2*(ny(x,z)*diff(ny(x,z),x)-nz(x,z)*diff(nx(x,z),z)+nz(x,z)*diff(nz(x,z),x))*diff(ny(x,z),x)+2*(nx(x,z)*diff(nx(x,z),z)-nx(x,z)*diff(nz(x,z),x)+ny(x,z)*diff(ny(x,z),z))*diff(ny(x,z),z))+q0*k22*(diff(nx(x,z),z)-diff(nz(x,z),x))-1/2*Ex*e0*de*nx(x,z)*Ey-1/2*Ey*e0*(de*nx(x,z)*Ex+2*de*ny(x,z)*Ey+de*nz(x,z)*Ez)-1/2*Ez*e0*de*nz(x,z)*Ey;

> diff(Fg,nyx);

> a2:=k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*nz(x,z)+1/2*k33*(2*(ny(x,z)*diff(ny(x,z),x)-nz(x,z)*diff(nx(x,z),z)+nz(x,z)*diff(nz(x,z),x))*ny(x,z)-2*(-nz(x,z)*diff(ny(x,z),z)-nx(x,z)*diff(ny(x,z),x))*nx(x,z))+q0*k22*nz(x,z);

> diff(Fg,nyz);

> b2:=-k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*nx(x,z)+1/2*k33*(-2*(-nz(x,z)*diff(ny(x,z),z)-nx(x,z)*diff(ny(x,z),x))*nz(x,z)+2*(nx(x,z)*diff(nx(x,z),z)-nx(x,z)*diff(nz(x,z),x)+ny(x,z)*diff(ny(x,z),z))*ny(x,z))-q0*k22*nx(x,z);

> EL_ny:=c2-diff(a2,x)-diff(b2,z);

> EL_ny_C:=simplify(%);

> C(EL_ny_C,ansi);

      s2 = 2.0*q0*k22*diff(nx(x,z),z)-2.0*q0*k22*diff(nz(x,z),x)-k33*ny(x,z)*nz

(x,z)*diff(diff(nz(x,z),x),x)-2.0*k33*nx(x,z)*nz(x,z)*diff(diff(ny(x,z),x),z)

-2.0*k33*nx(x,z)*diff(nx(x,z),x)*diff(ny(x,z),x)-k33*diff(nx(x,z),x)*nz(x,z)*

diff(ny(x,z),z)+k22*nx(x,z)*ny(x,z)*diff(diff(nx(x,z),z),z)-2.0*k22*nx(x,z)*

diff(ny(x,z),z)*diff(nx(x,z),z)+k22*nx(x,z)*diff(ny(x,z),z)*diff(nz(x,z),x)-4.0

*k22*ny(x,z)*diff(nx(x,z),z)*diff(nz(x,z),x);

      s1 = s2+k22*nz(x,z)*diff(ny(x,z),x)*diff(nx(x,z),z)-2.0*k22*nz(x,z)*diff(

ny(x,z),x)*diff(nz(x,z),x)-k33*nz(x,z)*diff(ny(x,z),x)*diff(nx(x,z),z)-k33*nx(x

,z)*diff(ny(x,z),z)*diff(nz(x,z),x)-Ey*Ey*e0*de*ny(x,z)+k22*nz(x,z)*diff(nx(x,z

),x)*diff(ny(x,z),z)+2.0*k22*nz(x,z)*nx(x,z)*diff(diff(ny(x,z),x),z)-k22*nz(x,z

)*ny(x,z)*diff(diff(nx(x,z),x),z)+k22*nz(x,z)*ny(x,z)*diff(diff(nz(x,z),x),x)-

k22*nx(x,z)*ny(x,z)*diff(diff(nz(x,z),x),z);

      s2 = s1+k22*nx(x,z)*diff(nz(x,z),z)*diff(ny(x,z),x)-2.0*k33*nz(x,z)*diff(

nz(x,z),z)*diff(ny(x,z),z)+k33*ny(x,z)*nz(x,z)*diff(diff(nx(x,z),x),z)-k33*diff

(nz(x,z),z)*nx(x,z)*diff(ny(x,z),x)-k33*ny(x,z)*nx(x,z)*diff(diff(nx(x,z),z),z)

+k33*ny(x,z)*nx(x,z)*diff(diff(nz(x,z),x),z)-k33*pow(nx(x,z),2.0)*diff(diff(ny(

x,z),x),x)-k33*pow(nz(x,z),2.0)*diff(diff(ny(x,z),z),z)-k33*pow(ny(x,z),2.0)*

diff(diff(ny(x,z),x),x)-k33*ny(x,z)*pow(diff(nz(x,z),x),2.0);

      t0 = s2-k22*pow(nx(x,z),2.0)*diff(diff(ny(x,z),z),z)-k33*ny(x,z)*pow(diff

(nx(x,z),z),2.0)-k33*pow(ny(x,z),2.0)*diff(diff(ny(x,z),z),z)+2.0*k22*ny(x,z)*

pow(diff(nz(x,z),x),2.0)+2.0*k22*ny(x,z)*pow(diff(nx(x,z),z),2.0)-k33*ny(x,z)*

pow(diff(ny(x,z),z),2.0)-k33*ny(x,z)*pow(diff(ny(x,z),x),2.0)-k22*pow(nz(x,z),

2.0)*diff(diff(ny(x,z),x),x)+2.0*k33*ny(x,z)*diff(nx(x,z),z)*diff(nz(x,z),x)-Ex

*e0*de*nx(x,z)*Ey-Ez*e0*de*nz(x,z)*Ey;

> diff(Fg,nz);

> c3:=k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*diff(ny(x,z),x)+1/2*k33*(2*(ny(x,z)*diff(ny(x,z),x)-nz(x,z)*diff(nx(x,z),z)+nz(x,z)*diff(nz(x,z),x))*(-diff(nx(x,z),z)+diff(nz(x,z),x))-2*(-nz(x,z)*diff(ny(x,z),z)-nx(x,z)*diff(ny(x,z),x))*diff(ny(x,z),z))+q0*k22*diff(ny(x,z),x)-1/2*Ex*e0*de*nx(x,z)*Ez-1/2*Ey*e0*de*ny(x,z)*Ez-1/2*Ez*e0*(de*nx(x,z)*Ex+de*ny(x,z)*Ey+2*de*nz(x,z)*Ez);

>

> diff(Fg,nzx);

> a3:=-k22*(-nx(x,z)*diff(ny(x,z),z)+ny(x,z)*diff(nx(x,z),z)-ny(x,z)*diff(nz(x,z),x)+nz(x,z)*diff(ny(x,z),x))*ny(x,z)+1/2*k33*(2*(ny(x,z)*diff(ny(x,z),x)-nz(x,z)*diff(nx(x,z),z)+nz(x,z)*diff(nz(x,z),x))*nz(x,z)-2*(nx(x,z)*diff(nx(x,z),z)-nx(x,z)*diff(nz(x,z),x)+ny(x,z)*diff(ny(x,z),z))*nx(x,z))-q0*k22*ny(x,z);

> diff(Fg,nzz);

> b3:=k11*(diff(nx(x,z),x)+diff(nz(x,z),z));

> EL_nz:=c3-diff(a3,x)-diff(b3,z);

> EL_nz_C:=simplify(%);

> C(EL_nz_C,ansi);

      s1 = -Ex*e0*de*nx(x,z)*Ez-Ey*e0*de*ny(x,z)*Ez+2.0*q0*k22*diff(ny(x,z),x)-

k11*diff(diff(nx(x,z),x),z)-k11*diff(diff(nz(x,z),z),z)+k33*pow(nz(x,z),2.0)*

diff(diff(nx(x,z),x),z)+2.0*k22*pow(diff(ny(x,z),x),2.0)*nz(x,z)+k33*nz(x,z)*

pow(diff(nx(x,z),z),2.0)-k33*nz(x,z)*pow(diff(nz(x,z),x),2.0)+k33*nz(x,z)*pow(

diff(ny(x,z),z),2.0)+k22*pow(ny(x,z),2.0)*diff(diff(nx(x,z),x),z)-k22*pow(ny(x,

z),2.0)*diff(diff(nz(x,z),x),x)-k33*nz(x,z)*pow(diff(ny(x,z),x),2.0)-k33*pow(nz

(x,z),2.0)*diff(diff(nz(x,z),x),x)+k33*pow(nx(x,z),2.0)*diff(diff(nx(x,z),x),z)

;

      s2 = s1-k33*pow(nx(x,z),2.0)*diff(diff(nz(x,z),x),x)-2.0*k22*diff(ny(x,z)

,x)*nx(x,z)*diff(ny(x,z),z)+3.0*k22*diff(ny(x,z),x)*ny(x,z)*diff(nx(x,z),z)-2.0

*k22*diff(ny(x,z),x)*ny(x,z)*diff(nz(x,z),x)-k33*ny(x,z)*diff(ny(x,z),x)*diff(

nx(x,z),z)+2.0*k33*diff(ny(x,z),z)*nx(x,z)*diff(ny(x,z),x)-Ez*Ez*e0*de*nz(x,z);

      t0 = s2-k22*ny(x,z)*diff(nx(x,z),x)*diff(ny(x,z),z)-k22*ny(x,z)*nx(x,z)*

diff(diff(ny(x,z),x),z)+k22*ny(x,z)*nz(x,z)*diff(diff(ny(x,z),x),x)+2.0*k33*nx(

x,z)*diff(nx(x,z),x)*diff(nx(x,z),z)-2.0*k33*nx(x,z)*diff(nx(x,z),x)*diff(nz(x,

z),x)+k33*nx(x,z)*ny(x,z)*diff(diff(ny(x,z),x),z)+k33*diff(nx(x,z),x)*ny(x,z)*

diff(ny(x,z),z)-k33*nz(x,z)*ny(x,z)*diff(diff(ny(x,z),x),x);

> gradV:=grad(V(x,z),[x,y,z]);

> e:=linalg[matrix](3,3,[e11(x,z),e12(x,z),e13(x,z),e21(x,z),e22(x,z),e23(x,z),e31(x,z),e32(x,z),e33(x,z)]);

> a4:=multiply(e,gradV);

> eq0:=diverge(a4,[x,y,z]);

> V2x2:=(vip1j+vim1j-2*vij)/dx^2;

> V2xz:=(vip1jp1+vim1jm1-vim1jp1-vip1jm1)/(4*dx*dz);

> V2z2:=(vijp1+vijm1-2*vij)/dz^2;

> Vx:=(vip1j-vim1j)/(2*dx);

> Vz:=(vijp1-vijm1)/(2*dz);

> de11dx:=(e11ip1j-e11im1j)/(2*dx);

> de13dx:=(e13ip1j-e13im1j)/(2*dx);

> de31dz:=(e31ijp1-e31ijm1)/(2*dz);

> de33dz:=(e33ijp1-e33ijm1)/(2*dz);

> EL_V:=de11dx*Vx+e11[i,j]*V2x2+de13dx*Vz+e13[i,j]*V2xz+de31dz*Vx+e31[i,j]*V2xz+de33dz*Vz+e33[i,j]*V2z2=0;

> EL_V2:=solve(EL_V,vij);

> C(EL_V2,ansi);

      s1 = -1.0/8.0;

      s4 = -dx*dz*e13ip1j*vijp1+dx*dz*e13ip1j*vijm1+dx*dz*e13im1j*vijp1-dx*dz*

e13im1j*vijm1-dx*dz*e31ijp1*vip1j+dx*dz*e31ijp1*vim1j+dx*dz*e31ijm1*vip1j-dx*dz

*e31ijm1*vim1j-e31[i-1][j-1]*dx*dz*vip1jp1-e31[i-1][j-1]*dx*dz*vim1jm1-dx*dx*

e33ijm1*vijm1+dx*dx*e33ijm1*vijp1+dx*dx*e33ijp1*vijm1-dx*dx*e33ijp1*vijp1;

      s3 = s4-dz*dz*e11im1j*vim1j+dz*dz*e11im1j*vip1j+dz*dz*e11ip1j*vim1j-dz*dz

*e11ip1j*vip1j-4.0*e11[i-1][j-1]*dz*dz*vip1j-4.0*e33[i-1][j-1]*dx*dx*vijm1-4.0*

e33[i-1][j-1]*dx*dx*vijp1+e31[i-1][j-1]*dx*dz*vim1jp1+e31[i-1][j-1]*dx*dz*

vip1jm1-e13[i-1][j-1]*dx*dz*vip1jp1-e13[i-1][j-1]*dx*dz*vim1jm1+e13[i-1][j-1]*

dx*dz*vim1jp1+e13[i-1][j-1]*dx*dz*vip1jm1-4.0*e11[i-1][j-1]*dz*dz*vim1j;

      s4 = 1/(e11[i-1][j-1]*dz*dz+e33[i-1][j-1]*dx*dx);

      s2 = s3*s4;

      t0 = s1*s2;

>