２１章：軸対称円筒座標の偏微分方程式（ラプラスの方程式）

作成2013.01.18
    　偏微分方程式（ラプラスの方程式）は温度分布や電圧分布で使用されます。
    　１２章で直交座標2次元のラプラスの方程式の計算を行いましたが、3次元での計算は面倒なものがあります。
    　ここでは、軸対称円筒座標の偏微分方程式（ラプラスの方程式）の計算を行ってみたいと思います。

  1. 計算方法（復習です）
    　ラプラスの方程式は

    で与えられる偏微分方程式です。今任意の関数φ(x,y)を考え、微小量hだけ位置を変化させた場合

    が成立します。（12.2）式から（12.5）式を加算すると

    となります。ここで、（12.1）式と（12.6）式から

    の関係式が得られます。しかし、（12.7）式は近似式であり正確ではありません。このため新しいφをφ'とし

    （12.8）式の計算を繰り返すことにより、誤差を低減する必要があります。また（12.8）式が成立するには、xの変化量hとｙの変化量hが等しい必要があります。


  2. 境界条件
    　（12.8）式の計算を実行するには、計算領域より小さい方向と大きい方向にひとつ余分な値を設定する必要があります。
    　この設定が境界条件です。境界条件の設定無しには（12.8）式の実行はできません。（12.1）式の解は無限に存在し、境界条件を設定しないと値が決定できません。
    　偏微分方程式はひとつの条件を規定しますが、変数が多数あるため厳密に境界条件を規定する必要があります。偏微分方程式が難解なのは、厳密に境界条件を規定する必要があるためです。


  3. 軸対称円筒座標の偏微分方程式（ラプラスの方程式）
    　軸対称円筒座標の偏微分方程式（ラプラスの方程式）は

    で与えられる偏微分方程式です。今任意の関数φ(x,y)を考え、微小量hだけ位置を変化させた場合

    が成立します。（21.2）式から（21.5）式を加算すると

    また、（21.2）式と（21.3）式を減算すると

    となります。（21.1）式に（21.6）式と（21.7）式を代入して整理すると

    の関係式が得られます。しかし、（21.8）式は近似式であり正確ではありません。このため新しいφをφ'とし

    （21.9）式の計算を繰り返すことにより、誤差を低減する必要があります。また（21.9）式が成立するには、ｒの変化量hとｚ の変化量hが等しい必要があります。。


  4. 境界条件
    　（21.9）式の計算を実行するには、計算領域より小さい方向と大きい方向にひとつ余分な値を設定する必要があります。
    　この設定が境界条件です。境界条件の設定無しには（21.9）式の実行はできません。（21.1）式の解は無限に存在し、境界条件を設定しないと値が決定できません。
    　偏微分方程式はひとつの条件を規定しますが、変数が多数あるため厳密に境界条件を規定する必要があります。偏微分方程式が難解なのは、厳密に境界条件を規定する必要があるためです。
    　なお、（21.9）式はｒ＝0が特異点となり、ｒ＝0をさけて計算領域を設定する必要があります。
  5. 軸対称偏微分方程式.xls(フリーソフト)のダウンロード
    1. ダウンロード
      　下記の軸対称偏微分方程式.xls(フリーソフト)]をダウンロードしてください。

      　ダウンロード後はダブルクリックで解凍してから使用してください。
      　       [軸対称偏微分方程式.xls(フリーソフト)]をダウンロードする。


    2. 軸対称偏微分方程式.xls(フリーソフト)
    3. ファイル構成
      （１） 軸対称偏微分方程式.xls　：フリーソフトです。
      （２）シート「Sheet1」：計算条件の設定を行います。
      （２）シート「Sheet2」：計算結果の出力されます。

    4. 注意事項
      （１）ファイルの保存場所の制限はありません。

    5. 標準的な実行方法
      （１）「軸対称偏微分方程式.xls」をダブルクリックで起動します。
      　　　（マクロを有効にして開いてください！！）
      （２）行数（ｊ）、列数(j)、変化幅、収束判定値、最大演算回数、演算領域指定を設定します。
      （３）境界条件表を設定します。（最外周は全て 境界条件とする必要があります。）
      （４）iがｚ方向、jが半径ｒ方向です。
      （５）境界条件表の最大は100×100です。
      （６）演算領域の初期値は境界条件(i=0,j=0)と同じ値に自動設定されます。
      （７）「計算実行」ボタンを押します。
      （８）計算結果がシート「Sheet2」に表示されます。


  6. プログラムリスト
    　VBAのプログラムリストは以下のようになります。
    Dim T(101, 101), DT(101, 101), C(101, 101)
    Sub Main()
    'Sheet2クリア
    Sheets("Sheet2").Select
    Cells.Select
    Selection.Clear
    '初期入力
    Ni = Sheets("Sheet1").Cells(4, 3).Value
    Nj = Sheets("Sheet1").Cells(5, 3).Value
    h = Sheets("Sheet1").Cells(6, 3).Value
    e = Sheets("Sheet1").Cells(7, 3).Value
    Nc = Sheets("Sheet1").Cells(8, 3).Value
    Ce = Sheets("Sheet1").Cells(9, 3).Value
    '計算部分
    For i = 0 To Ni - 1 Step 1
    For j = 0 To Nj - 1 Step 1
    F = Sheets("Sheet1").Cells(15 + i, 2 + j).Value
    If F = Ce Then
    C(i, j) = 1
    T(i, j) = T(0, 0)
    Else
    C(i, j) = 0
    T(i, j) = F
    End If
    Next j
    Next i

    For k = 0 To Nc Step 1
    Nk = 0
    Se = 0
    For i = 0 To Ni - 1 Step 1
    For j = 0 To Nj - 1 Step 1
    If C(i, j) = 1 Then
    DT(i, j) = (T(i, j + 1) - T(i, j - 1)) / (8 * j) + (T(i - 1, j) + T(i + 1, j) + T(i, j - 1) + T(i, j + 1)) / 4
    Nk = Nk + 1
    Se = Se + Abs(DT(i, j) - T(i, j))
    Else
    End If
    Next j
    Next i

    For i = 0 To Ni - 1 Step 1
    For j = 0 To Nj - 1 Step 1
    If C(i, j) = 1 Then T(i, j) = DT(i, j)
    Next j
    Next i
    '収束状況出力
    Sheets("Sheet2").Cells(5 + k, 1).Value = k
    Sheets("Sheet2").Cells(5 + k, 2).Value = Se / Nk
    If Se / Nk < e Then Exit For
    Next k
    '結果出力
    Sheets("Sheet2").Cells(1, 1).Formula = "計算結果"
    Sheets("Sheet2").Cells(4, 1).Formula = "計算回数"
    Sheets("Sheet2").Cells(4, 2).Formula = "変化量"
    For j = 0 To Nj - 1 Step 1
    Sheets("Sheet2").Cells(4, 4 + j).Value = j
    Next j

    For i = 0 To Ni - 1 Step 1
    Sheets("Sheet2").Cells(5 + i, 3).Value = i
    For j = 0 To Nj - 1 Step 1
    Sheets("Sheet2").Cells(5 + i, 4 + j).Value = T(i, j)
    Next j
    Next i
    End Sub


  7. サンプル条件の実行結果
    偏微分方程式1（ラプラスの方程式）（温度分布の計算）
    計算条件
    項目記号値
    行数（ｊ）Ni21
    列数(j)Nj21
    変化幅h1
    収束判定値e0.01
    最大演算回数Nc400
    演算領域指定Ce-99

    *境界条件表(最大100×100まで）
    -ｊ＝0ｊ＝1ｊ＝2ｊ＝3ｊ＝4ｊ＝5ｊ＝6ｊ＝7ｊ＝20ｊ＝8ｊ＝9ｊ＝10ｊ＝11ｊ＝12ｊ＝13ｊ＝14ｊ＝15ｊ＝16ｊ＝17ｊ＝18ｊ＝20
    i＝0000000000000000000000
    i＝10-99-99-99-99-99-99-990-99-10-10-10-10-99-99-99-99-99-990
    i＝20-99-99-99-99-99-99-990-99-10-10-10-10-99-99-99-99-99-990
    i＝30-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝40-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝50-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝60-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝70-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝80-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝90-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝100-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝110-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝120-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝130-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝140-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝150-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝160-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝170-99-99-99-99-99-99-990-99-99-99-99-99-99-99-99-99-99-990
    i＝180-99-99-99-99-99-99-990-9910101010-99-99-99-99-99-990
    i＝190-99-99-99-99-99-99-990-9910101010-99-99-99-99-99-990
    i＝20000000000000000000000

    計算結果
    *計算回数毎の変化量が表示されます
    計算回数変化量
    00.116639
    10.080584
    20.068019
    30.059061
    40.053723
    50.049315
    60.046253
    70.043537
    80.041493

    計算結果表
    -01234567208910111213141516171820
    0000000000000000000000
    10-0.22-0.33-0.46-0.66-0.98-1.53-2.580-4.77-10-10-10-10-4.29-2.06-1.07-0.59-0.33-0.180
    20-0.4-0.59-0.82-1.15-1.65-2.46-3.790-6.06-10-10-10-10-5.4-3-1.71-0.99-0.57-0.310
    30-0.51-0.75-1.02-1.4-1.94-2.72-3.830-5.33-6.98-7.42-7.31-6.62-4.61-2.97-1.87-1.15-0.69-0.380
    40-0.55-0.8-1.08-1.43-1.92-2.55-3.360-4.28-5.1-5.41-5.28-4.69-3.61-2.56-1.72-1.12-0.7-0.40
    50-0.52-0.76-1-1.31-1.69-2.17-2.720-3.28-3.73-3.9-3.78-3.36-2.72-2.04-1.45-0.98-0.63-0.370
    60-0.46-0.65-0.85-1.09-1.38-1.71-2.080-2.41-2.66-2.75-2.65-2.38-1.98-1.54-1.13-0.79-0.52-0.310
    70-0.35-0.51-0.65-0.83-1.02-1.25-1.470-1.68-1.82-1.86-1.79-1.61-1.36-1.08-0.81-0.58-0.39-0.230
    80-0.24-0.34-0.44-0.55-0.67-0.8-0.940-1.05-1.13-1.15-1.1-1-0.85-0.68-0.52-0.38-0.25-0.160
    90-0.12-0.17-0.22-0.27-0.33-0.39-0.450-0.51-0.54-0.54-0.52-0.47-0.41-0.33-0.25-0.18-0.13-0.080
    10000000000000000000000
    1100.120.170.220.270.330.390.4500.510.540.540.520.470.410.330.250.180.130.080
    1200.240.340.440.550.670.80.9401.051.131.151.110.850.680.520.380.250.160
    1300.350.510.650.831.021.251.4701.681.821.861.791.611.361.080.810.580.390.230
    1400.460.650.851.091.381.712.0802.412.662.752.652.381.981.541.130.790.520.310
    1500.520.7611.311.692.172.7203.283.733.93.783.362.722.041.450.980.630.370
    1600.550.81.081.431.922.553.3604.285.15.415.284.693.612.561.721.120.70.40
    1700.510.751.021.41.942.723.8305.336.987.427.316.624.612.971.871.150.690.380
    1800.40.590.821.151.652.463.7906.06101010105.431.710.990.570.310
    1900.220.330.460.660.981.532.5804.77101010104.292.061.070.590.330.180
    20000000000000000000000

    計算結果グラフ

    　図21-1において、横軸が半径ｒ方向、縦軸が高さｚ方向となります。左側が 円筒の内側となります。


  8. まとめ
    　軸対称円筒電極において、電圧分布は円筒の内側に向かうことがわかります。 これから、円筒の内側の方が電界が強くなることがわかります。

２２章：フロ−ティング電極の電位分布に行く。



トップページに戻る。