Potential in Open Rectangular Channel with Parallel Sides Grounded and Closed End at Potential (Subscript[V, o])

PROBLEM:

In two-dimensional Cartesian coordinates, Laplace’s Equation becomes

		(((∂^2)V)/(∂(x^2))) + (((∂^2)V)/(∂(y^2))) = 0

Consider an infinitely long, two-dimensional channel formed by two grounded (conducting) plates, parallel to the z-x plane and separated by a distance a, that extend to ± infinity in z-direction and from 0 to +∞ in the x-direction.  The channel between the two plates is terminated on the left by a non-conducting plate parallel to the y-z plane.  The potential on the end-plate is a function of y only and is given by 

		Vzero[y] =  Vo Sin[πy/a].  

Note that the Vzero[y] potentials on the left-hand plate are measured with respect to the grounded upper and lower plates.

(a) In the text cell below, explain why is this a two dimensional problem.  A symmetry argument is useful.

<Enter your explanation in this text cell>

(b) Use M to plot the function Vzero[y] from 0 to a; remember to plot a function you have to give any constants (here, a and Vo) reasonable values.  (M occasionally stumbles over variable names with subscripts.  As a precaution, I replace Vo with Vo in M code.)

(* Input code below in this input cell *)

(c) State the 2D boundary conditions that we will need to find V[x,y].

<List the boundary conditions in this text cell>

(d) No Brainer-to nail an important point:  Click inside the cell below (or select it by clicking on the bracket to the right) and execute it (Shift-Return);  Answer Boxes will appear; Click on the one you think is correct answer for this question:

	Question:  What is the boundary condition on the right hand boundary (as x → ∞)?

Button["1  V inside the channel goes to ∞ as x-> ∞.",{Print["  Wrong.  Sorry.   "]}]
Button["2 V inside the channel goes to 0 as x-> ∞.",{Print["  Correct!!    "]}]
Button["3  V inside the channel goes to Vo as x-> ∞.",{Print["  Wrong.  Sorry.   "]}]

(d-continued) Justify your answer.

<Enter your justification in this text cell>

A slight aside:  How would you “construct” the electrode structures to approximate the problem we are working on??

First, all the plate coordinates that are supposed to be infinite, for practical purposes, means "much larger than a".  So we cut two LARGE rectangles or squares of nice, flat sheet metal for these electrodes and mount them in parallel orientation.  They could be insulated but in final form grounded so they are both at V = 0.

To approximate the plate holding the y dependent Vzero[y] is trickier.  

Here is one scheme:  Consider a set of very long parallel wires insulated from each other and from ground; label the wires with some index, say j.  Let J (odd number) be the total number of wires.  The spacing of the wires would be a/J so that the gap between the parallel plates (at x = 0) would be filled in uniform increments of y.

To get the Vo Sin[((π y)/a)], we would apply the potential Vo on the center wire and connect resistors between the wires with magnitudes that would drop the potential on the wires in sync with the Sin[((π y)/a)] function as you increase and decrease j from J/2 (going from the center (y = a/2) to y = 0 and y = a).  

Finding the magnitudes of these R's is assigned to the "loudest mouth in the class".

Details on Separation of Variables:

Unlike the 1-D Laplace’s equation, the 2-D equation is not an ordinary differential equation.  That is, the derivatives are not total, but “partial”.  The bad news is that partial differential equations seldom have closed, analytic solutions.  The good news is that you can often construct an (infinite) series whose sum satisfies the boundary conditions.   

The most straightforward way to find solutions to a partial differential equation is to “separate the variables” in x and y.  To separate the variables, we assume that the solution of the differential equation can be written as the product of two parts, one that depends on x only and one that depends on y only.

V(x,y) = X(x) Y(y)

Solutions that cannot be written in this form can usually be written as the sum of solutions with this form.  When V(x,y) is substituted in Laplace’s equation,

(((∂^2)V)/(∂(x^2))) + (((∂^2)V)/(∂(y^2))) = (((∂^2)X)/(∂(x^2)))Y(y) + X(x)(((∂^2)Y)/(∂(y^2))) + (((∂^2)X)/(∂(x^2)))Y(y) + X(x)(((∂^2)Y)/(∂(y^2))) = (((∂^2)X)/(∂(x^2)))Y(y) + X(x)(((∂^2)Y)/(∂(y^2))) = 0

(By construction dX/dy = dY/dx = 0.)  If we divide through by the product, X(x)Y(y) we get two terms—one depending only on x and the other depending only on y.  

(1/(X(x)))(((∂^2)X)/(∂(x^2))) + (1/(Y(x)))(((∂^2)Y)/(∂(y^2))) = 0

[Problems may occur at points where X(x)=0 and/or Y(y)=0.  We can often work around such points if they are isolated.]  Since x and y are independent variables, it is unreasonable to expect that a function that depends only on x to always be equal (with opposite sign) to a function that depends on y—unless the term in x is constant and the term in y equals (–1) times that constant.  One term will be positive and the other negative.  The sign seriously impacts the solutions X(x) and Y(y), and the wrong choice makes it impossible to satisfy the boundary conditions.  In practice, you quickly run through the solutions for both sets of signs, then choose the solutions that work for your problem.  

Since the choice of signs is important, we will set one constant to be something definitely positive, say, the square of a real number, (λ^2).  The other constant will be –(λ^2).  For instance, 

(1/(X(x)))(((∂^2)X)/(∂(x^2))) = (λ^2)       (1/(Y(x)))(((∂^2)Y)/(∂(y^2)))= -(λ^2)

This procedure has reduced the partial differential equation above to two, independent, ordinary differential equations.  We can solve ordinary differential equations.  Any pair of solutions X(x) and Y(y) that satisfy the above ordinary differential equations will satisfy the partial differential equations above.  The easiest way to solve the above equations is “by inspection.”  Since each equation is second order, each equation has two independent solutions.

(e) In the space below, solve the two ordinary differential equations, either by inspection, or by using M.

<If you solve the equations by inspection, enter your solutions in this text cell>

Using M.  You will need to assume that λ is real and is greater than or equal to zero.  (You can omit the possibility that λ = 0 for now, because the corresponding solution would be trivial and uninteresting.)

(*If you use M to find the solutions, enter your code below in this cell, below $Assumptions  *)

$Assumptions → λ ∈ Reals && λ > 0&&y ∈ Reals &&  a ∈ Reals && a > 0  && 0 ⩽ y ≲ a   && (Subscript[V, o]) ∈ Reals && (Subscript[V, o]) > 0; 
    (* Leave assumptions *)

To this point, we have four constants of integration and the unknown λ that are constrained by the boundary conditions.  

We require V[x →∞,y] = 0.  If C[1]≠0 in X[x], the Exp[λx] term will go to infinity and violate this constraint.  Therefore C[1]=0 in X[x] no matter what λ is. (Note that the Exp[+λx] terms are needed to satisfy boundary conditions in other problems.  Some problems require both the Exp[-λx] and Exp[+λx] terms.)   

We require V[x,y=0] = 0.  If If C[1]≠0 in Y[y], the Cos[λy] term will violate this constraint, so C[1] in this equation must also be zero. Our solution must take the form:

X(x) = X2 Exp[–λx], and 
Y(y) = Y2 Sin[λy]

But we also expect V[x,y=a] = 0.  The Sin[λx] term will allow this only if λ = nπ/a, where n is a positive integer (n≠0). This constraint poses a problem in that we have an infinite number of λ’s.  It also poses an opportunity, in that we may need an infinite number of them to satisfy our remaining boundary condition:  Vzero[y] = (Subscript[V, o]) Sin[πy/a].  

We proceed by assuming that Vzero[y] can be expressed as the sum of functions that satisfy the other boundary conditions; that is, for some set of coefficients {(Subscript[C, n]) , n = 1,2,3, ....}

Vzero[y] = (∑_(n=1)^∞)(Subscript[C, n]) Sin[((n π y)/a)]

This particular series expansion is known as a Fourier sine series.  The existence and convergence of this kind of series is discussed in detail in most Applied or Engineering Mathematics texts.  The constants {(Subscript[C, n])} can be determined using the “Fourier trick”.  The Fourier trick takes advantage of the fact that the set of functions {sin[nπy/a]} is orthogonal over the interval {0,a}.  That is, the integral of the product of two members of this set over the interval {0,a} is zero unless they have the same n-value.  This easily shown using the trig identity sinθ sinϕ = (1/2)[cos(θ-ϕ) - sin(θ+ϕ)].

(∫_0^a)Sin[((n π y)/a)]Sin[((n'π y)/a)] dy  = (1/2)(∫_0^a)Cos[(((n-n')π y)/a)] dy - (1/2)(∫_0^a)Sin[(((n+n')π y)/a)] dy

Now the integral of any sine or cosine function over any number of complete cycles is equal to zero.  Therefore both integrals equal zero for every choice of {n,n'} except for n = n', when Cos[0] = 1.  Then the integral equals a/2.  In summary:

(∫_0^a)Sin[((n π y)/a)]Sin[((n'π y)/a)] dy  =({{{, ({{0, n≠n' }, {a/2, n=n'}})}})

This identity allows us to compute the coefficients {(Subscript[C, n])}.  If we multiply each side of the series expansion for V[x] by sin[n'π y/a] and Integrate from 0 to a, we have:

(∫_0^a)V[y] Sin[((n'πy)/a)]dy = (∑_(n=1)^∞)(Subscript[C, n]) (∫_0^a)Sin[((nπy)/a)]Sin[((n'πy)/a)]dy=((a (Subscript[C, n']))/2)
(Subscript[C, n'])= (2/a)(∫_0^a)V[y]Sin[((n'πy)/a)]dy

Determining Cn for n = 1,2,3...  produces an infinite series that converges to the solution of the 2-D Laplace’s equation that meets all the boundary conditions. (In general, the series may fail to converge at a finite number of points.  We ignore these points.) The solution will include both X[x] and Y[y] factors, so that:

V[x,y] = X[x] Y[y] = (∑_(n=1)^∞)(Subscript[C, n]) Exp[-((n π x)/a)] Sin[((n π y)/a)]

We need to determine (Subscript[C, n]) for for the remaining boundary condition Vzero[y] = (Subscript[V, o]) Sin[πy/a].  Since Vzero[y] includes only one term, and that term corresponds to the n=1 term in our infinite series.  Direct integration shows that (Subscript[C, 1]) = (Subscript[V, o]).  More complicated cases benefit from M.

(f) Use M to confirm that (Subscript[C, 1]) = Vo.

(* Input code below in this input cell *)
ClearAll["`*"] (* Leave the ClearAll["`*"] statement and assumptions *)
$Assumptions n ∈ Integers && n > 0 && x ∈ Reals && x > 0 && y ∈ Reals && y > 0 && a ∈ Reals && a > 0;

More complicated solutions require a number of (Subscript[C, n]) (= CC[n] below) values.  In general, all four solutions to Laplace’s equation [two for X(x) and two for Y(y)] will be needed.  While this can get complicated, it can be shown that virtually any function f(x) can be written as the sum of an infinite sine series and an infinite cosine series (with possible exceptions at isolated points).  That is, the set of functions {CC[n] Sin[πx/a] + DD[n] Cos[πx/a]} is complete.  

The same is true of the set of functions {A[n] Exp[πy/a] + B[n] Exp[-πy/a]} is complete:  any function g(y) can be written in terms of an infinite series of terms with this form.  

When a large number of coefficients must be calculated, the Table[] command is useful.  For instance, the CC[n] values for n = 1...10 are tabulated below.  I use CC[n] because the variable C is protected in M.

(* Execute code below *)
ClearAll["`*"]
$Assumptions -> n ∈ Integers && y ∈ Reals && a ∈ Reals;
CC[n_] =  ((2 Vo)/a)(∫_0^a)Sin[((π y)/a)]Sin[((n π y)/a)]\[DifferentialD]y
Table[CC[n],{n, 1, 10}]

The values for n > 1 are all zero, as expected, but the general function given by M for aa[n] fails for n = 1, where the denominator of  ((2 Vo Sin[n π])/(π-(n^2) π)) is zero.  One can show by L’Hospital’s rule that it yields the correct value, (Subscript[V, o]), in the limit as n → 1.   M can also evaluate the limit for you:

(* Execute code below *)
CC[1] = Limit[CC[n], n-> 1]

(g)  Write down (as input) your final solution for the given Vzero[y] = Vo Sin[((π y)/a)] channel problem.  Then test that this solution indeed satisfies Laplace’s Equation and the Boundary Conditions you determined above in (c).

(* Input code below in this input cell *)
ClearAll["`*"]  (* Leave ClearAll["`*"] statement *)

So the solution given by the Fourier trick has the required Laplacian (zero, the first zero above) and meets all three boundary conditions for the grounded surfaces without condition.  The condition on the surface with potential distribution Vo Sin[πy/a] is also met, with a condition, namely  a > 0.  Since the problem requires a to be positive and nonzero, we’re cool.

(h) Plot V[x,y] using the Plot3D command and the contour plot.  For plotting, you will need to set Vo and a to some reasonable values.  I chose Vo = a = 1.   Comment on the results.  For instance, do the boundary conditions look like they are being satisfied?  (Remember: you can use your mouse to rotate the Plot3D graph; you can click on the contours to see their magnitude.)

(* Input code below in this input cell *)
Clear[Vo,a]; (* Leave Clear[] statement *)

<Enter comments in this text cell>

(i) Calculate the x- and y- components of the electric field corresponding to this geometry.

(* Input code below in this input cell below the Clear statement*)
Clear[Vo,a]; (* Leave Clear[] statement *)

(j)  Create a VectorPlot of the resulting electric field.  You will have to give values to  Vo and a.  While the VectorPlot is useful, it is limited by the space required for the longer vector arrows.  It is often difficult to visualize the magnitude and direction of the field in regions of strong and weak fields on the same plot.  

In the past, we have superimposed the ContourPlot over a StreamPlot of the field get a more complete picture of the field.  The background color and spacing of the potential contours reflects the magnitude of the field, while the StreamPlot indicates its direction.  In addition to the vector plot, superimpose the corresponding ContourPlot and StreamPlot.  Comment.

(* Input code below in this input cell *)

<Enter your comments in this text cell>

Both the vector plot and the contour plot suggest that the magnitude of the electric field just next to the left hand plate (where x ≈ 0) is constant:  the length of the vector arrows and the spacing of the potential contours are uniform to the eye—as near as I can see.  The constant magnitude can be confirmed by calculation, using the field components for x = 0.   

For a quick look at this, we can use a DensityPlot of the Magnitude of EE (this is given by Norm[EE[x,y])  I superimpose the VectorPlot and the DensityPlot below.  You can see the bands of color are surprisingly uniform for a given value of y.  Note that the vertical rows of EE arrows, although differing in direction, have similar lengths.

(* Execute code below *)
density = DensityPlot[Norm[EE[x,y]],{x,0,a},{y,0,a},AspectRatio->1,ColorFunction->"SunsetColors"];
Show[density,vec]

(k) Charge Density.  For –∞ < z < ∞, the charge density on the left hand surface of the channel (at x = 0) is a function of y only (σleft[y]).  Find σleft[y] in units of C/(m^2).  (You will need to introduce (Subscript[ϵ, 0]) = 8.85×(10^(-12)) F/m).  Note the magnitude of the values of σ displayed.  As before, I assumed Vo = a = 0.  

Plot your resulting σleft[y].

(* Input code below in this input cell *)

(Subscript[ϵ, 0]) = 8.85×(10^(-12)); a=1; (Subscript[V, o]) = 1;   (* (Subscript[ϵ, 0]) in F/m; a in m; (Subscript[V, o]) in V *)
EE[x_,y_] = {(E^(-π x)) π Sin[π y],-(E^(-π x)) π Cos[π y]}
σleft[y_] = (Subscript[ϵ, 0]) EE[0,y][[1]]
Plot[σleft[y],{y,0, a},AxesLabel->{"y", "σleft (C-(m^(-2)))"}, PlotRange-> {0, 3 (10^(-11))}]

(l) Now find the charge density on the bottom surface of the channel (at y = 0) [For 0 < x < ∞, this will be a function of x only (σbottom[x]). Again, find σbottom[x] in units of C/(m^2) (so you will need to introduce (Subscript[ϵ, 0]) = 8.85×(10^(-12)) F/m again).  Note the magnitude of the values of σ displayed.  Plot your result for σbottom[x].  Comment on your results.

(* Input code below in this input cell *)

<Enter comments in this text cell>

(m) I claim that the net charge of the system as a whole is zero.  

We have a system where –∞ > z > ∞ so the total charge on the left hand surface would be infinite; likewise for the charges on the horizontal plates.  Nevertheless, it is fair to compare the charges per unit length, λ, in the z-direction.  

So first, calculate the charges/unit length (λ) in the z direction by integrating σleft over y and σbottom over x.   Compare them.  Are you happy with the result?

(* Input code below in this input cell *)

<Enter your comparison in this text cell>