ELECTROSTATICS

Remove["Global`*"]

Highlights

Abstract

Introduction

Needs["Graphics`PlotField`"];

Off[General::"spell", General::"spell1"]; Needs["Graphics`Arrow`"]; tera = (10^12); micro = (1/(10^6)); giga = (10^9); nano = (1/(10^9)); mega = (10^6); pico = (1/(10^12)); kilo = (10^3); femto = (1/(10^15)); milli = (1/(10.^3)); atto = (1/(10^18)); NSub
       = {c -> 2.99792458 (10^8), e -> (1.6021773/(10^19)), 
        h -> (6.6260755/(10.^34)), hbar -> (1.05457/(10.^34)), 
        K -> (1.380658/(10^23)), me -> (9.109389/(10^31)), 
        mu0 -> ((4. π)/(10^7)), eps0 -> (8.85419/(10^12))}; ESi = 0.179 tera; rhoSi = 2330; nuSi = 0.28; ESiO2 = 0.06 tera; rhoSiO2 = 2200; nuSiO2 = 0.3; ESi3N4 = 0.15 tera; rhoSi3N4 = 3100; nuSi3N4 = 0.24; EAu = 0.08 tera; rhoAu = 19300; nuAu = 0.42; ENi = 0.207 tera; rhoNi = 8900; nuNi = 0.41; EW = 0.345 tera; rhoW = 19300; nuW = 0.28; EPt = 0.17 tera; rhoPt = 21400; nuPt = 0.39; EIr = 0.52 tera; rhoIr = 22500; nuIr = 0.3; EHOPG = 0.2 tera; rhoHOPG = 2250; nuHOPG = 0.3; EPol = 0.001 tera; rhoPol = 1000; nuPol = 0.3; xx = 0.8; EPtIr = xx 0.8 EPt + (1 - xx) EIr; rhoPtIr = xx rhoPt + (1 - xx) rhoIr; nuPtIr = 0.3; \(opt = {PlotRange -> All, AxesOrigin -> {0, 0}, 
        Frame -> True, GridLines -> Automatic, 
        DefaultFont -> {"Times-Roman", 8}, 
        PlotStyle 
          -> {Dashing[{.0001}], Dashing[{.0500}], Dashing[{.0200}], 
            Dashing[{.0001}]}}; 

    \(opt1 = {DefaultFont -> {"Times-Roman", 8}, 
        FrameLabel -> {"x", "y"}, Frame -> True}; \)\)

Isolated Point Charge

Clear[vpoint]; vpoint[x_, y_] := (1/(Sqrt[(x^2) + (y^2)]));

discGraphics = 
      Graphics[Disk[{0, 0}, .1], AxesOrigin -> {0, 0}, 
        AspectRatio -> Automatic, PlotLabel -> "a", opt1];

fieldGraphics = 
      PlotGradientField[
        -vpoint[x, y], {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
        ScaleFunction -> (2 + Log[1 + #1]&), 
        ScaleFactor -> 0.5, PlotPoints -> 8, 
        DisplayFunction -> Identity];

discfieldGraphics = 
      Show[{discGraphics, fieldGraphics}, DisplayFunction -> Identity];

potentialGraphics = 
      ContourPlot[vpoint[x, y], {x, -2.5, 2.5}, {y, -2.5, 2.5}, 
        PlotPoints -> 20, Contours -> 7, PlotLabel -> "b", 
        ContourShading -> False, DisplayFunction -> Identity, 
        Evaluate[opt1]];

Fig1 = Show[GraphicsArray[{discfieldGraphics, potentialGraphics}]];

Point Charge and Plane

case1 = {d -> 1}; Qi = -1; si = -2 d /. case1; Clear[vreal, vimage]; vreal[x_, y_] := (1/(Sqrt[(x^2) + (y^2)])); vimage[x_, y_] := -(1/(Sqrt[(x^2) + (((y - si))^2)]));

discGraphics = 
      Graphics[{Disk[{0, 0}, .1], 
          Rectangle[{-2, -1.1}, {2, -1}]}, 
        AspectRatio -> Automatic, PlotLabel -> "a", opt1];

fieldGraphics = 
      PlotGradientField[
        -(vreal[x, y] + vimage[x, y]), {x, -1.5, 1.5}, {y, 
          -1, 1}, ScaleFunction -> (2 + Log[1 + #1]&), 
        ScaleFactor -> 0.5, PlotPoints -> 8, 
        DisplayFunction -> Identity];

discfieldGraphics = 
      Show[{discGraphics, fieldGraphics}, DisplayFunction -> Identity];

potentialGraphics = 
      ContourPlot[
        vreal[x, y] + vimage[x, y], {x, -2.5, 2.5}, {y, -3.5, 1.5}, 
        PlotPoints -> 20, Contours -> 10, 
        PlotLabel -> "b", ContourShading -> False, 
        DisplayFunction -> Identity, Evaluate[opt1]];

Fig2 = Show[GraphicsArray[{discfieldGraphics, potentialGraphics}]];

Point Charge and Sphere

case2 = {R -> 1., d -> 2.}; Qi = -(R/d) /. case2; si = d - ((R^2)/d) /. case2; Clear[vreal, vimage]; vreal[x_, y_] := (1/(Sqrt[(x^2) + (y^2)])); vimage[x_, y_] := (Qi/(Sqrt[(x^2) + (((y + si))^2)]));

discGraphics = 
      Graphics[{Disk[{0, 0}, .1], Disk[{0, -d /. case2}, 1]}, 
        AspectRatio -> Automatic, PlotLabel -> "a", opt1];

fieldGraphics = 
      PlotGradientField[
        -(vreal[x, y] + vimage[x, y]), {x, -2, 2}, {y, -3.5, 
          1.5}, ScaleFunction -> (2 + Log[1 + #1]&), 
        ScaleFactor -> 0.5, PlotPoints -> 8, 
        DisplayFunction -> Identity];

discfieldGraphics = 
      Show[{discGraphics, fieldGraphics}, DisplayFunction -> Identity];

potentialGraphics = 
      ContourPlot[
        vreal[x, y] + vimage[x, y], {x, -2.5, 2.5}, {y, -3.5, 1.5}, 
        PlotPoints -> 30, Contours -> 10, 
        PlotLabel -> "b", ContourShading -> False, 
        DisplayFunction -> Identity, Evaluate[opt1]];

Fig3 = Show[GraphicsArray[{discfieldGraphics, potentialGraphics}]];

Isolated Sphere

Clear[Csp]; Csp[R_] := 4 π eps0 R;

case3 = {R -> 10 nano}; Print["Csp = ", N[Csp[R]/atto] /. case3 /. NSub, " [aF]"]

Sphere and Plane

case4 = {R -> 1., d -> 2.}; nn = 5;

Position of Charges Inside the Sphere

h = Array[q, nn]; h[[1]] = 0; h[[2]] = ((R^2)/(2 d)); tbl1a = 
      Simplify[Table[
          h[[n]] = 
            ((R^2)/(2 d - 
                h[[n - 1]])), {n, 3, 
            nn}]]; p11 = Table[h[[n]], {n, 1, nn}]; 
    p1 = Take[p11, 3]

Magnitude of Charges Inside the Sphere

g = Array[q, nn]; g[[1]] = 1; g[[2]] = r; tbl2a = 
      Simplify[Table[
          g[[n]] = 
            ((r g[[n - 1]])/(1 - 
                ((r g[[n - 1]])/g)
                    [[n - 2]])), {n, 3, 
            nn}]]; c11 = 
      Table[g[[n]], {n, 1, nn}] /. 
        r -> (R/(2 d)); c1 = Take[c11, 3]

Position of Charges Outside the Sphere

p22 = 
      Table[2 d - h[[n]], {n, 1, nn}]; 
    p2 = Take[p22, 3]

Magnitude of Charges Outside the Sphere

c22 = 
      Table[-g[[n]], {n, 1, nn}] /. 
        r -> (R/(2 d)); c2 = Take[c22, 3]

Table of Coefficients

Clear[tbl1]; scx = N[c11 /. case4, 3]; sc4 = Take[scx, 4]; spx = N[p11 /. case4, 3]; sp4 = Take[spx, 4]; pcx = N[c22 /. case4, 3]; pc4 = Take[pcx, 4]; ppx = N[p22 /. case4, 3]; pp4 = Take[ppx, 4];

tbl1 = 
      TableForm[{sc4, sp4, pc4, pp4}, 
        TableHeadings 
          -> {{"c1", "p1", "c2", "p2"}, {"1", 
              "2", "3", "4"}}]

Potential and Field

Clear[vsp]; vsp[x_, y_] := 
      (∑_(i = 1)^nn)(
          scx[[i
                ]]/(Sqrt[(x^2) + 
                  (((y + spx[[i]]))^2)
                    ]) + pcx[[i
                ]]/(Sqrt[(x^2) + 
                  (((y + ppx[[i]]))^2)
                    ]));

discGraphics = 
      Graphics[{Disk[{0, 0}, 1], Rectangle[{-2, -2.1}, {2, -2}]}, 
        AspectRatio -> Automatic, PlotLabel -> "a", opt1];

fieldGraphics = 
      PlotGradientField[-vsp[x, y], {x, -1.5, 1.5}, {y, -2, 1.5}, 
        ScaleFunction -> (2 + Log[1 + #1]&), 
        ScaleFactor -> 0.5, PlotPoints -> 8, 
        DisplayFunction -> Identity];

discfieldGraphics = 
      Show[{discGraphics, fieldGraphics}, DisplayFunction -> Identity];

potentialGraphics = 
      ContourPlot[vsp[x, y], {x, -2.5, 2.5}, {y, -3.5, 1.5}, 
        PlotPoints -> 30, Contours -> 10, 
        PlotLabel -> "b", ContourShading -> False, 
        DisplayFunction -> Identity, Evaluate[opt1]];

Fig4 = Show[GraphicsArray[{discfieldGraphics, potentialGraphics}]];

Potential Along the Axis of Symmetry

RR = R /. case4; dd = d /. case4; Fig5 = 
      Plot[vsp[0, y], {y, -dd, -RR}, 
        FrameLabel -> {"y", "V"}, Evaluate[opt]];

Capacitance

Clear[Csp]; Csp := 
      Function[{R, 
          s}, {alpha = Log[1 + (s/R) + (Sqrt[((s^2)/(R^2)) + ((2 s)/R)])]; 
            cap2 = 1 + (∑_(n = 2)^100) Sinh[alpha]/Sinh[n alpha]; 
            N[4. π eps0 R cap2] //. NSub}[[1
          ]]];

Clear[Cspa]; Cspa[R_, s_] := N[2 π eps0 R (2 + (R/s))] /. NSub;

Example

case5 = {R -> 10 nano, s -> 50 nano}; RR = R /. case5; ss = s /. case5; Timing[Print["Csp  = ", Csp[RR, ss]/atto, " [aF]"]]Timing[Print["Cspa = ", Cspa[RR, ss]/atto, " [aF]"]]

Fig6 = 
      Plot[Cspa[RR, s nano]/Csp[RR, s nano], {s, (RR/nano), (ss/nano)}, 
        FrameLabel -> {"s in nm", "Cspa over Csp"}, 
        Evaluate[opt]];

Two Spheres

Capacitance: Exact Solution

Clear[R1, R2, s]; Css := 
      Function[{R1, R2, 
          s}, {temp1 = 4. π eps0 R1 //. NSub; 
            temp2 = 4. π eps0 R2 //. NSub; R1sq = R1 R1; 
            R2sq = R2 R2; q = 1.; da = s + R1 + R2; d1 = da; qt = 0.; 
            qqt = 0.; 
            Do[{qq = -((R2 q)/da); dd = d1 - (R2sq/da); qt = qt + q; 
                qqt = qqt + qq; q = -((R1 qq)/dd); 
                da = d1 - (R1sq/dd); }, {i, 1, 100}]; c11 = qt temp1; 
            c12 = qqt temp1; q = 1.; da = s + R1 + R2; qt = 0.; qqt = 0.; 
            Do[{qq = -((R1 q)/da); dd = d1 - (R1sq/da); qt = qt + q; 
                q = -((R2 qq)/dd); da = d1 - (R2sq/dd); }, {i, 1, 
                100}]; c22 = qt temp2; c21 = c12; 
            N[((c11 c22 - c12 c21)/(c11 + c12 + c21 + c22))]}
            [[1]]]

Example

case6
       = {R1 -> 100 nano, R2 -> 10 nano, s -> 1 nano}; RR1 = R1 /. case6; RR2 = R2 /. case6; ss = s /. case6; Print["Css = ", Css[RR1, RR2, ss]/atto, " [aF]"]

Matrix Elements

tbl2 = 
      TableForm[{{N[(c11/atto)], N[(c12/atto)]}, {N[(c21/atto)], N[(c22/atto)]}}, 
        TableHeadings 
          -> {{"j=1", "j=2"}, {"i=1", "i=2"}}]

Capacitance: Approximate Solution

Clear[Cssa]; Cssa = 
      Function[{R1, R2, 
          s}, {R11 = Max[{R1, R2}]; R22 = Min[{R1, R2}]; 
            rb = (1./(1. + (R22/R11))); N[rb Csp[R22, rb s]]}
            [[1]]];

Example

Timing[Print["Css  = ", Css[RR1, RR2, ss]/atto, " [aF]"]]Timing[Print["Cssa = ", Cssa[RR1, RR2, ss]/atto, " [aF]"]]

Fig7 = 
      Plot[Cssa[RR1, RR2, s nano]/Css[RR1, RR2, s nano], {s, .1, 100.}, 
        FrameLabel -> {"s in nm", "Cssa over Css"}, 
        PlotRange -> {.9, 1.}, Evaluate[opt]];

Electrostatic Force

Clear[F]; F[R_, s_, v_] = -(1/2) (∂_s) Cspa[R, s] (v^2);

Fig8 = 
      Plot[F[100 nano, s nano, 1]/nano, {s, 10, 100}, 
        FrameLabel -> {"s in nm", "F in nN"}, Evaluate[opt]];