プレチルト角の決定

または，

max ^³ | ∆α

^(ν)₁

/α

^(ν)₁

| , | ∆α

^(ν)₂

/α

^(ν)₂

| , | ∆α

^(ν)₃

/α

₃^(ν)

| , | ∆γ

₃^(ν)

/γ

₃^(ν)

| ^´ < δ

_vis

(B.44)

ならば計算終了，そうでなければ反復継続とする．

− β

5

+ γ

2

√ β

1

√ β

₁

tan θ

o

H(θ

o

) cos

²

θ

o

#

+ 2α

1

β

6

sin 2θ

o

√ β

1

à 1

J(θ

o

) − 1 M (θ

o

)

!)

(B.49)

d

²

t(θ

o

)

dθ

²_o

= 1 2ǫ

o

∆ǫE

_b²

( β

2

β

3

"

− 4 cos 2θ

o

sin

²

2θ

_o

− − β

5

+ γ

2

√ β

1

√ β

1

G(θ

o

) − tan θ

o

( − G(θ

o

) sin 2θ

o

+ 2β

4

tan θ

o

) G

²

(θ

o

) cos

⁴

θ

o

− β

₅

+ γ

₂

√ β

₁

√ β

1

H(θ

_o

) − tan θ

_o

( − H(θ

_o

) sin 2θ

_o

+ 2β

₄

tan θ

_o

) H

²

(θ

o

) cos

⁴

θ

o

#

+ 4α

1

β

6

√ β

1

"

cos 2θ

_o

à 1

J(θ

o

) − 1 M (θ

o

)

!

+α

1

sin

²

2θ

o

à 1

J

²

(θ

_o

) − 1 M

²

(θ

_o

)

!#)

(B.50) β

1

< 0

のとき，

dt(θ

o

) dθ

o

= 1

ǫ

o

∆ǫE

_b²

( β

2

β

3

"

1

sin 2θ

o

− (γ

2

N (θ

o

) + β

5

) tan θ

o

(N

²

(θ

o

) − β

1

) cos

²

θ

o

#

− 2α

1

β

6

sin 2θ

o

P (θ

o

)

)

(B.51)

d

²

t(θ

o

)

dθ

²_o

= 1 ǫ

o

∆ǫE

_b²

( β

2

β

3

"

− 2 cos 2θ

o

sin

²

2θ

_o

− γ

2

(3β

4

tan

²

θ

o

+ β

7

)(N

²

(θ

o

) − β

1

) (N

²

(θ

o

) − β

1

)

²

cos

⁴

θ

o

+γ

2

N (θ

o

) tan θ

o

4β

4

N (θ

o

) tan θ

o

− (N

²

(θ

o

) − β

1

) sin 2θ

o

(N

²

(θ

o

) − β

1

)

²

cos

⁴

θ

o

+β

5

4β

4

N(θ

o

) tan

²

θ

o

− (N

²

(θ

o

) − β

1

)(2 sin

²

θ

o

+ 1) (N

²

(θ

o

) − β

1

)

²

cos

⁴

θ

o

#

− 4α

1

β

6

P (θ

o

) cos 2θ

o

+ 2α

1

sin

²

2θ

o

(α

1

cos 2θ

o

+ γ

2

) P

²

(θ

o

)

)

(B.52)

となる．G(θo

)，H(θ

o

)，J(θ

o

)，M (θ

o

)，N (θ

o

)，および P (θ

o

)

はそれぞれ

G(θ

o

) = β

4

tan

²

θ

o

+ β

7

+ ^q β

1

(B.53)

H(θ

o

) = β

4

tan

²

θ

o

+ β

7

− ^q β

1

(B.54)

J (θ

o

) = α

1

cos 2θ

o

+ γ

2

+ ^q β

1

(B.55) M(θ

o

) = α

1

cos 2θ

o

+ γ

2

− ^q β

1

(B.56) N (θ

o

) = β

4

tan

²

θ

o

+ β

7

(B.57) P (θ

o

) = (α

1

cos 2θ

o

+ γ

2

)

²

− β

1

(B.58)

である．

θ

o

= θ

^∗_oで

S

t

(θ

^∗_o

)

が最小になるとすれば，式

(B.46)

より，

S

t

(θ

_o^∗

) = minS

t

(θ

o

) = minS

t

(θ

^(ν)_o

+ ∆θ

_o^(ν)

)

= min

"

S

t

(θ

_o^(ν)

) + dS

t

(θ

o

) dθ

_o

¯

¯

¯

¯

¯

_θ

o=θ_o^(ν)

∆θ

^(ν)_o

+ 1

2

d

²

S

t

(θ

o

) dθ

_o²

¯

¯

¯

¯

¯

θo=θ⁽o^ν)

(∆θ

^(ν)_o

)

²

#

(B.59)

となる．したがって，St

(θ

_o^∗

)

は

∆θ

^(ν)_o の

2

次式のうちの最小値となるから，θ^(ν+1)_o は

∆θ

_o^(ν)

= −

dS

_t

(θ

_o

) dθ

o

¯

¯

¯

¯

¯

θo=θ⁽o^ν)

d

²

S

t

(θ

o

) dθ

²_o

¯

¯

¯

¯

¯

θo=θo^(ν)

(B.60)

として，θ^(ν+1)_o

= θ

_o^(ν)

+ ∆θ

_o^(ν)で決まる．以上の手順を

| ∆θ

^(ν)_o

|

または

| ∆θ

^(ν)_o

/θ

^(ν)_o

|

が十 分に小さくなるまで繰り返す．具体的には，収束判定用の定数

δ

tiltを用意しておいて，

| ∆θ

^(ν)_o

| < δ

tilt または

| ∆θ

^(ν)_o

/θ

^(ν)_o

| < δ

tilt

(B.61)

ならば計算終了，そうでなければ反復継続とする．

参考文献

[1] M. Imai, H. Naito, M. Okuda, and A. Sugimura, Jpn. J. Appl. Phys. 33 , L 119 (1994).

[2] M. Imai, H. Naito, M. Okuda, and A. Sugimura, Jpn. J. Appl. Phys. 33 , 3482 (1994).

[3]

小国力, Fortran 95, C & Javaによる新数値計算法

–数値計算とデータ分析–,

サイ エンス社

(1997).

[4]

戸川隼人,新装版

UNIX

ワークステーションによる科学技術計算ハンドブック

[基

礎篇

C

言語版],サイエンス社

(1998).

®

­

©

付 録 C ª

プログラムのソースコード

 本論文で示した水平配向および垂直配向液晶セルに流れる過渡電流の数値計算，お よび垂直配向液晶セルに流れる過渡電流の実験結果と解析解との最小二乗フィッティン

グは，

Visual Basic

を用いて作成したプログラムにより行った．ここでは，プログラム

のソースコードを掲載する．ソースコード上の変数は，全て倍精度浮動小数点型で宣 言した．

C.1 液晶セルに流れる過渡電流の数値計算

式

(2.16)

，

(2.18)

，

(A.18)

，

(A.19)

，

(A.42)

，および

(A.45)

を数値的に解くことによ り，ダイレクタのツイスト変形も考慮して液晶セルに流れる過渡電流波形を計算できる．

no slip

境界条件において過渡電流を数値計算するプログラムのソースコードを以下に

示す．

no slip

境界条件では液晶セル基板上で流れは発生しない．この条件を満たすよう

に式

(A.18)

および

(A.19)

の

σ

x

(t)

，

σ

y

(t)

を計算しなければならないため，

free slip

境界 条件，流れを無視した場合と比較して計算過程が複雑になる．ここで，

σ

x

(t)

，

σ

y

(t)

は ともに時間

t

における

z

に対する積分定数である．

free slip

境界条件および流れを無視 した場合における過渡電流は，それぞれ

σ

x

(t) = σ

y

(t) = 0

および

v

x

(z, t) = v

y

(z, t) = 0

として計算すればよい．ここで，

v

x

(z, t)

および

v

y

(z, t)

は基板に対して

z

軸が垂直にな るように

3

次直交座標系を配置した際に，それぞれ

x

および

y

軸方向に発生する流速 である．また，ダイレクタの方位角

φ(z, t) = 0

とおいて計算することにより，ダイレ

クタのツイスト変形が発生しない場合の過渡電流波形を得ることができる．

”Form1”

Option Explicit ’変数の宣言を強制 Private Sub Form_Load()

Text1.Text = ""

Text2.Text = ""

Text3.Text = ""

Text4.Text = ""

Text5.Text = ""

Text6.Text = ""

Text7.Text = ""

Text8.Text = ""

Text9.Text = ""

Text10.Text = ""

Text11.Text = ""

Text12.Text = ""

Text13.Text = ""

Text14.Text = ""

Text15.Text = ""

Text16.Text = ""

Text17.Text = ""

Text18.Text = ""

Text19.Text = ""

Text20.Text = ""

Label1.Caption = "α1 (Pa s)" ’Leslie粘性係数 Label2.Caption = "α2 (Pa s)"

Label3.Caption = "α3 (Pa s)"

Label4.Caption = "α4 (Pa s)"

Label5.Caption = "α5 (Pa s)"

Label6.Caption = "θ0_0 (deg)" ’z=0でのプレチルト角 Label7.Caption = "θ0_L (deg)" ’z=Lでのプレチルト角

Label8.Caption = "ε⊥" ’ダイレクタに垂直方向の誘電定数 Label9.Caption = "Δε" ’誘電率異方性

Label10.Caption = "S (cm^2)" ’電極面積 Label11.Caption = "L (μm)" ’セル厚 Label12.Caption = "V (V)" ’印加電圧 Label13.Caption = "k11 (pN)" ’スプレイ弾性定数 Label14.Caption = "k22 (pN)" ’ツイスト弾性定数 Label15.Caption = "k33 (pN)" ’ベンド弾性定数

Label16.Caption = "n_o" ’ダイレクタに垂直方向の屈折率 Label17.Caption = "n_e" ’ダイレクタに水平方向の屈折率 Label18.Caption = "λ(nm)" ’入射光波長

Label19.Caption = "Δt (s)" ’微小時間幅 Label20.Caption = "セル分割数"

Command1.Caption = "START"

π = 3.14159265358979 ’円周率

ε0 = 1 / (4 * π* (2.99792458) ^ 2 * 1000000000#) ’真空中の誘電率(F/m) End Sub

Private Sub Command1_Click()

α1 = Text1.Text ’Leslie粘性係数(Pa s)

α2 = Text2.Text α3 = Text3.Text α4 = Text4.Text α5 = Text5.Text

θ00 = Text6.Text * π / 180 ’z=0でのプレチルト角(rad) θ0L = Text7.Text * π / 180 ’z=Lでのプレチルト角(rad)

ε⊥ = Text8.Text ’ダイレクタに垂直方向の誘電定数

Δε = Text9.Text ’誘電率異方性

S = Text10.Text * 0.0001 ’電極面積(m^2) L = Text11.Text * 0.000001 ’セル厚(m)

V = Text12.Text ’印加電圧(V)

k11 = Text13.Text * 0.000000000001 ’スプレイ弾性定数(N) k22 = Text14.Text * 0.000000000001 ’ツイスト弾性定数(N) k33 = Text15.Text * 0.000000000001 ’ベンド弾性定数(N)

n_o = Text16.Text ’ダイレクタに垂直方向の屈折率

n_e = Text17.Text ’ダイレクタに水平方向の屈折率

λ = Text18.Text * 0.000000001 ’入射光波長(m)

Δt = Text19.Text ’微小時間幅(s)

z = Text20.Text ’セル分割数

ε0Δε = ε0 * Δε

γ1 = α3 - α2 ’回転粘性率(Pa s) γ2 = α2 + α3

Δz = L / z ’セル厚分割幅(m)

rate = 2 * π * Δz / λ

o_r = Cos(n_o * rate) ’常光exp(i*n_o*rate)の実数成分 o_i = Sin(n_o * rate) ’常光exp(i*n_o*rate)の虚数成分 α = 45 ’偏光子とx軸との角度(deg)

’極角分布の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "z_θ.xy " For Output As #1 Print #1, "z", "θ"

Print #1, ""

Print #1, ""

’方位角分布の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "z_φ.xy " For Output As #2 Print #2, "z", "φ"

Print #2, ""

Print #2, ""

’電場分布の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "z_E.xy " For Output As #3 Print #3, "z", "E"

Print #3, ""

Print #3, ""

’x方向の流速分布の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "z_vx.xy " For Output As #4 Print #4, "z", "vx"

Print #4, ""

Print #4, ""

’y方向の流速分布の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "z_vy.xy " For Output As #5 Print #5, "z", "vy"

Print #5, ""

Print #5, ""

’過渡電流の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "t_I.xy " For Output As #6 Print #6, "t", "I"

Print #6, ""

Print #6, ""

’透過光強度の計算結果書き込み用ファイルのオープン

Open "D:\vb_NLC\V" + Text12.Text + "t_intensity.xy " For Output As #7 Print #7, "t", "intensity"

Print #7, ""

Print #7, ""

’t=0sについての計算//////////////////////////////////////////////////////////////////////

’初期配向の計算//////////////////////////////////////////////////////////////////////////

For iz = 0 To z

θ(iz) = (θ0L - θ00) * iz / z + θ00 ’ダイレクタの極角

φ(iz) = 0 ’ダイレクタの方位角

Next iz

For orientation = 0 To 1000000 Caption = orientation

Command1.Caption = orientation For iz = 1 To z - 1

dθdt2 = (k11 - k33) * (θ(iz + 1) ^ 2 - 2 * θ(iz + 1) * θ(iz - 1) _ + θ(iz - 1) ^ 2) / (2 * Δz) ^ 2

dθdt3 = -(2 * k22 * Sin(θ(iz)) ^ 2 + k33 * (2 * Cos(θ(iz)) ^ 2 - 1)) _

* (φ(iz + 1) ^ 2 - 2 * φ(iz + 1) * φ(iz - 1) + φ(iz - 1) ^ 2) _ / (2 * Δz) ^ 2

θ(iz) = (θ(iz + 1) + θ(iz - 1) + (dθdt2 + dθdt3) * Sin(θ(iz)) _

* Cos(θ(iz)) * Δz ^ 2 / (k11 * Sin(θ(iz)) ^ 2 _ + k33 * Cos(θ(iz)) ^ 2)) / 2

dφdt2 = 2 * (2 * k22 * Sin(θ(iz)) ^ 2 + k33 * (2 * Cos(θ(iz)) ^ 2 - 1)) _

* Cos(θ(iz)) / Sin(θ(iz))

dφdt3 = (θ(iz + 1) - θ(iz - 1)) * (φ(iz + 1) - φ(iz - 1)) _ / (2 * Δz) ^ 2

φ(iz) = (φ(iz + 1) + φ(iz - 1) + dφdt2 * dφdt3 * Δz ^ 2 _ / (k22 * Sin(θ(iz)) ^ 2 + k33 * Cos(θ(iz)) ^ 2)) / 2 Next iz

Next orientation

’////////////////////////////////////////////////////////////////////////////////////////

’光強度を計算する行列(z=0)

p1ra = 1 ’1_1成分の実数成分

p1rb = 0 ’1_2成分の実数成分

p1rc = 0 ’2_1成分の実数成分

p1rd = 1 ’2_2成分の実数成分

p1ia = 0 ’1_1成分の虚数成分

p1ib = 0 ’1_2成分の虚数成分

p1ic = 0 ’2_1成分の虚数成分

p1id = 0 ’2_2成分の虚数成分

vx = 0 ’x方向の流速

vy = 0 ’y方向の流速

εrev_t_1 = 0 ’誘電率の逆数の和 For iz = 0 To z

ε(iz) = ε0 * (ε⊥ + Δε * Cos(θ(iz)) ^ 2) ’ネマティック液晶の誘電率 εrev_t_1 = εrev_t_1 + (1 / ε(iz))

If iz = 0 Then

’光を常光成分と異常光成分に分割する行列R_0の成分 rs = -Cos(φ(0))

rc = Sin(φ(0)) Else

’光を常光成分と異常光成分に分割する行列R_iの成分 rs = Sin(φ(iz) - φ(iz - 1))

rc = Cos(φ(iz) - φ(iz - 1)) End If

’異常光成分の全光線が感じる屈折率

n_e_gam = n_o * n_e / Sqr((n_o * Sin(θ(iz))) ^ 2 + (n_e * Cos(θ(iz))) ^ 2) e_r = Cos(n_e_gam * rate) ’異常光exp(i*n_e_gam*rate)の実数成分 e_i = Sin(n_e_gam * rate) ’異常光exp(i*n_e_gam*rate)の虚数成分

’第i番目の分割点の行列R_iと第i-1番目の分割点に至るまでの光強度を表す行列との積 p2ra = rc * p1ra + rs * p1rc ’1_1成分の実数成分

p2rb = rc * p1rb + rs * p1rd ’1_2成分の実数成分 p2rc = -rs * p1ra + rc * p1rc ’2_1成分の実数成分 p2rd = -rs * p1rb + rc * p1rd ’2_2成分の実数成分 p2ia = rc * p1ia + rs * p1ic ’1_1成分の虚数成分 p2ib = rc * p1ib + rs * p1id ’1_2成分の虚数成分 p2ic = -rs * p1ia + rc * p1ic ’2_1成分の虚数成分 p2id = -rs * p1ib + rc * p1id ’2_2成分の虚数成分 If iz = z Then

’第z番目の分割点に至るまでの光強度を表す行列 p1ra = p2ra ’1_1成分の実数成分

p1rb = p2rb ’1_2成分の実数成分 p1rc = p2rc ’2_1成分の実数成分 p1rd = p2rd ’2_2成分の実数成分 p1ia = p2ia ’1_1成分の虚数成分 p1ib = p2ib ’1_2成分の虚数成分 p1ic = p2ic ’2_1成分の虚数成分 p1id = p2id ’2_2成分の虚数成分 Else

’第i番目の分割点に至るまでの光強度を表す行列

p1ra = o_r * p2ra - o_i * p2ia ’1_1成分の実数成分 p1rb = o_r * p2rb - o_i * p2ib ’1_2成分の実数成分 p1rc = e_r * p2rc - e_i * p2ic ’2_1成分の実数成分 p1rd = e_r * p2rd - e_i * p2id ’2_2成分の実数成分 p1ia = o_r * p2ia + o_i * p2ra ’1_1成分の虚数成分 p1ib = o_r * p2ib + o_i * p2rb ’1_2成分の虚数成分 p1ic = e_r * p2ic + e_i * p2rc ’2_1成分の虚数成分 p1id = e_r * p2id + e_i * p2rd ’2_2成分の虚数成分 End If

Print #4, iz * Δz, vx ’x方向の流速分布の計算結果書き込み Print #5, iz * Δz, vy ’y方向の流速分布の計算結果書き込み Next iz

Q = V / (εrev_t_1 * Δz) ’セルに蓄えられいる単位面積あたりの電荷 For iz = 0 To z

E(iz) = Q / ε(iz) ’ポアソン方程式 If iz = 0 Or iz = z Then

dθdt(iz) = 0 dφdt(iz) = 0 Else

dθdt1 = (k11 * Sin(θ(iz)) ^ 2 + k33 * Cos(θ(iz)) ^ 2) * (θ(iz + 1) _

- 2 * θ(iz) + θ(iz - 1)) / Δz ^ 2

dθdt2 = (k11 - k33) * (θ(iz + 1) ^ 2 - 2 * θ(iz + 1) * θ(iz - 1) _ + θ(iz - 1) ^ 2) / (2 * Δz) ^ 2

dθdt3 = -(2 * k22 * Sin(θ(iz)) ^ 2 + k33 * (2 * Cos(θ(iz)) ^ 2 - 1)) _

* (φ(iz + 1) ^ 2 - 2 * φ(iz + 1) * φ(iz - 1) + φ(iz - 1) ^ 2) _ / (2 * Δz) ^ 2

dθdt(iz) = (dθdt1 + (dθdt2 + dθdt3 - ε0Δε * E(iz) ^ 2) _

* Sin(θ(iz)) * Cos(θ(iz))) / γ1

dφdt1 = (k22 * Sin(θ(iz)) ^ 2 + k33 * Cos(θ(iz)) ^ 2) _

* (φ(iz + 1) - 2 * φ(iz) + φ(iz - 1)) / Δz ^ 2

dφdt2 = 2 * (2 * k22 * Sin(θ(iz)) ^ 2 + k33 * (2 * Cos(θ(iz)) ^ 2 - 1)) _

* Cos(θ(iz)) / Sin(θ(iz))

dφdt3 = (θ(iz + 1) - θ(iz - 1)) * (φ(iz + 1) - φ(iz - 1)) _ / (2 * Δz) ^ 2

dφdt4 = dφdt2 * dφdt3

dφdt(iz) = (dφdt1 + dφdt4) / γ1 End If

Print #1, iz * Δz, θ(iz) * 180 / π ’極角分布の計算結果書き込み Print #2, iz * Δz, φ(iz) * 180 / π ’方位角分布の計算結果書き込み

Print #3, iz * Δz, E(iz) ’電場分布の計算結果書き込み

Next iz

rs = Sin(φ(z)) rc = Cos(φ(z))

pra = rs * p1ra + rc * p1rc prb = rs * p1rb + rc * p1rd prc = -rc * p1ra + rs * p1rc prd = -rc * p1rb + rs * p1rd pia = rs * p1ia + rc * p1ic pib = rs * p1ib + rc * p1id pic = -rc * p1ia + rs * p1ic pid = -rc * p1ib + rs * p1id αs = Sin(α * π / 180) αc = Cos(α * π / 180)

pr = αc * (-αs * pra + αc * prc) + αs * (-αs * prb + αc * prd) pi = αc * (-αs * pia + αc * pic) + αs * (-αs * pib + αc * pid) intensity = pr ^ 2 + pi ^ 2 ’透過光強度

Print #7, 0, intensity

’////////////////////////////////////////////////////////////////////////////////////////

’t=0s以降の計算//////////////////////////////////////////////////////////////////////////

For it = 1 To 100000

If Cos(it * 2 * π / 10) = 1 Then ’時間ステップの経過を10ステップごとに表示 Caption = it

Command1.Caption = it Debug.Print it End If

t = Δt * it p1ra = 1 p1rb = 0 p1rc = 0 p1rd = 1 p1ia = 0 p1ib = 0 p1ic = 0 p1id = 0

εrev = 0 ’誘電率の逆数の和

For iz = 0 To z

θ(iz) = θ(iz) + Δt * dθdt(iz) φ(iz) = φ(iz) + Δt * dφdt(iz) sinθ(iz) = Sin(θ(iz))

cosθ(iz) = Cos(θ(iz)) sinφ(iz) = Sin(φ(iz)) cosφ(iz) = Cos(φ(iz))

ε(iz) = ε0 * (ε⊥ + Δε * cosθ(iz) ^ 2)

a(iz) = α1 * (sinθ(iz) * cosθ(iz)) ^ 2 + (-γ2 * (2 * cosθ(iz) ^ 2 - 1) _ + α3 + α4 + α5) / 2

b(iz) = γ2 * cosθ(iz) ^ 2 - α3

c(iz) = ((α5 - α2) * cosθ(iz) ^ 2 + α4) / 2 εrev = εrev + (1 / ε(iz))

If iz = 0 Then rs = -Cos(φ(0)) rc = Sin(φ(0))

Elsers = Sin(φ(iz) - φ(iz - 1)) rc = Cos(φ(iz) - φ(iz - 1))

End If

n_e_gam = n_o * n_e / Sqr((n_o * sinθ(iz)) ^ 2 + (n_e * cosθ(iz)) ^ 2) e_r = Cos(n_e_gam * rate)

e_i = Sin(n_e_gam * rate) p2ra = rc * p1ra + rs * p1rc p2rb = rc * p1rb + rs * p1rd p2rc = -rs * p1ra + rc * p1rc p2rd = -rs * p1rb + rc * p1rd p2ia = rc * p1ia + rs * p1ic p2ib = rc * p1ib + rs * p1id p2ic = -rs * p1ia + rc * p1ic p2id = -rs * p1ib + rc * p1id If iz = z Then

p1ra = p2ra p1rb = p2rb p1rc = p2rc p1rd = p2rd p1ia = p2ia p1ib = p2ib p1ic = p2ic p1id = p2id Else

p1ra = o_r * p2ra - o_i * p2ia p1rb = o_r * p2rb - o_i * p2ib p1rc = e_r * p2rc - e_i * p2ic p1rd = e_r * p2rd - e_i * p2id p1ia = o_r * p2ia + o_i * p2ra p1ib = o_r * p2ib + o_i * p2rb p1ic = e_r * p2ic + e_i * p2rc p1id = e_r * p2id + e_i * p2rd End If

Next iz rs = sinφ(z) rc = cosφ(z)

pra = rs * p1ra + rc * p1rc prb = rs * p1rb + rc * p1rd prc = -rc * p1ra + rs * p1rc prd = -rc * p1rb + rs * p1rd pia = rs * p1ia + rc * p1ic pib = rs * p1ib + rc * p1id pic = -rc * p1ia + rs * p1ic pid = -rc * p1ib + rs * p1id αs = Sin(α * π / 180) αc = Cos(α * π / 180)

pr = αc * (-αs * pra + αc * prc) + αs * (-αs * prb + αc * prd) pi = αc * (-αs * pia + αc * pic) + αs * (-αs * pib + αc * pid) intensity = pr ^ 2 + pi ^ 2

Q = V / (εrev * Δz)

I = S * V * ((1 / εrev) - (1 / εrev_t_1)) / (Δz * Δt) ’過渡電流 εrev_t_1 = εrev

If Cos(it * 2 * π / 10) = 1 Then

Print #6, t - Δt, I ’過渡電流の計算結果書き込み

Print #7, t, intensity ’透過光強度の計算結果書き込み End If

For iz = 0 To z E(iz) = Q / ε(iz)

If Cos(it * 2 * π / 10000) = 1 Then Print #1, iz * Δz, θ(iz) * 180 / π Print #2, iz * Δz, φ(iz) * 180 / π Print #3, iz * Δz, E(iz)

End If Next iz

’no-slip境界条件を満たす積分定数σx，σyを計算///////////////////////////////////////

σx_1 = 0 σy_1 = 0 Δσ = 100

Do Until Δσ < 0.00001 σ1 = 0

σ2 = 0 σ3 = 0 σ4 = 0 σ5 = 0 σ6 = 0 σ7 = 0 For iz = 0 To z

σ1 = σ1 + b(iz) * sinφ(iz) * dθdt(iz) / a(iz)

σ2 = σ2 + b(iz) * cosφ(iz) * dθdt(iz) / a(iz)

σ3 = σ3 + sinθ(iz) * cosθ(iz) * sinφ(iz) * dφdt(iz) / c(iz) σ4 = σ4 + sinθ(iz) * cosθ(iz) * cosφ(iz) * dφdt(iz) / c(iz) σ5 = σ5 + (sinφ(iz) ^ 2 / a(iz) + cosφ(iz) ^ 2 / c(iz)) σ6 = σ6 + (cosφ(iz) ^ 2 / a(iz) + sinφ(iz) ^ 2 / c(iz)) σ7 = σ7 + (1 / a(iz) - 1 / c(iz)) * sinφ(iz) * cosφ(iz) Next iz

σx = (σ2 * σ5 - σ1 * σ7 - α2 * (σ3 * σ5 + σ4 * σ7)) _ / (σ5 * σ6 - σ7 ^ 2)

σy = (σ1 * σ6 - σ2 * σ7 + α2 * (σ4 * σ6 + σ3 * σ7)) _ / (σ5 * σ6 - σ7 ^ 2)

For iz = 0 To z

If iz = 0 Or iz = z Then dθdt(iz) = 0 dφdt(iz) = 0

Elseθvxvy = (σx * cosφ(iz) + σy * sinφ(iz) - b(iz) * dθdt(iz)) / a(iz) dθdt1 = (k11 * sinθ(iz) ^ 2 + k33 * cosθ(iz) ^ 2) _

* (θ(iz + 1) - 2 * θ(iz) + θ(iz - 1)) / Δz ^ 2

dθdt2 = (k11 - k33) * (θ(iz + 1) ^ 2 - 2 * θ(iz + 1) * θ(iz - 1) _ + θ(iz - 1) ^ 2) / (2 * Δz) ^ 2

dθdt3 = -(2 * k22 * sinθ(iz) ^ 2 + k33 * (2 * cosθ(iz) ^ 2 - 1)) _

* (φ(iz + 1) ^ 2 - 2 * φ(iz + 1) * φ(iz - 1) _ + φ(iz - 1) ^ 2) / (2 * Δz) ^ 2

dθdt4 = (dθdt2 + dθdt3 - ε0Δε * E(iz) ^ 2) * sinθ(iz) * cosθ(iz) dθdt5 = (γ1 - γ2 * (2 * cosθ(iz) ^ 2 - 1)) *θvxvy / 2

dθdt(iz) = (dθdt1 + dθdt4 + dθdt5) / γ1

φvxvy = (-σx * sinφ(iz) + σy * cosφ(iz) - α2 * sinθ(iz) _

* cosθ(iz) * dφdt(iz)) / c(iz)

dφdt1 = (k22 * sinθ(iz) ^ 2 + k33 * cosθ(iz) ^ 2) _

* (φ(iz + 1) - 2 * φ(iz) + φ(iz - 1)) / Δz ^ 2

dφdt2 = 2 * (2 * k22 * sinθ(iz) ^ 2 + k33 * (2 * cosθ(iz) ^ 2 - 1)) _

* cosθ(iz) / sinθ(iz)

dφdt3 = (θ(iz + 1) - θ(iz - 1)) * (φ(iz + 1) - φ(iz - 1)) _ / (2 * Δz) ^ 2

dφdt4 = dφdt2 * dφdt3

dφdt5 = (γ1 - γ2) * cosθ(iz) * φvxvy / (2 * sinθ(iz)) dφdt(iz) = (dφdt1 + dφdt4 + dφdt5) / γ1

End If Next iz

Δσ = Abs(σx - σx_1)

If Abs(σx - σx_1) < Abs(σy - σy_1) Then Δσ = Abs(σy - σy_1)

End If σx_1 = σx σy_1 = σy Loop

’////////////////////////////////////////////////////////////////////////////////////

’流速の計算//////////////////////////////////////////////////////////////////////////

If Cos(it * 2 * π / 10000) = 1 Then For iz = 0 To z

If iz = 0 Or iz = z Then vx = 0

vy = 0 Else

vx1 = σx * (cosφ(iz - 1) ^ 2 / a(iz - 1) + sinφ(iz - 1) ^ 2 _ / c(iz - 1))

vx2 = σy * (1 / a(iz - 1) - 1 / c(iz - 1)) _

* sinφ(iz - 1) * cosφ(iz - 1)

vx3 = -b(iz - 1) * cosφ(iz - 1) * dθdt(iz - 1) / a(iz - 1) vx4 = α2 * sinθ(iz - 1) * cosθ(iz - 1) * sinφ(iz - 1) _

* dφdt(iz - 1) / c(iz - 1) vx = vx + Δz * (vx1 + vx2 + vx3 + vx4) vy1 = σx * (1 / a(iz - 1) - 1 / c(iz - 1)) _

* sinφ(iz - 1) * cosφ(iz - 1)

vy2 = σy * (sinφ(iz - 1) ^ 2 / a(iz - 1) + cosφ(iz - 1) ^ 2 _ / c(iz - 1))

vy3 = -b(iz - 1) * sinφ(iz - 1) * dθdt(iz - 1) / a(iz - 1) vy4 = -α2 * sinθ(iz - 1) * cosθ(iz - 1) * cosφ(iz - 1) _

* dφdt(iz - 1) / c(iz - 1) vy = vy + Δz * (vy1 + vy2 + vy3 + vy4) End If

Print #4, iz * Δz, vx ’x方向の流速分布の計算結果書き込み

Print #5, iz * Δz, vy ’y方向の流速分布の計算結果書き込み Next iz

End If

’////////////////////////////////////////////////////////////////////////////////////

Next it

’////////////////////////////////////////////////////////////////////////////////////////

Close #1 Close #2 Close #3 Close #4 Close #5 Close #6 Close #7

Command1.Caption = "END"

End Sub

C.2 垂直配向液晶セルに流れる過渡電流の実験結果と解析

ドキュメント内 負の誘電率異方性を有するネマティック液晶の過渡応答に関する研究 (ページ 134-145)