ME340A Elasticity Spring 2010 Wei Cai Stanford University
11
4/23/06 8:13 PM
C:\Documents and Settings\Wei Cai\My Documents\Courses...\S522a.m
%adpated from maple file S522 from J. R. Barber, Elasticity % http://www-personal.engin.umich.edu/~jbarber/elasticity/maple/S522 % % adapted by Wei Cai,
[email protected], for ME 340 Elasticity % Spring 2006, Stanford University % %This file gives the solution of the simply-supported beam problem, %following the strategy of Section 5.2.2. % clear all syms C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 C16 C17 C18 syms x y a b p %stress function phi = C1*x^2+C2*x*y+C3*y^2+C4*x^3+C5*x^2*y+C6*x*y^2+C7*y^3+C8*x^4+ ... C9*x^3*y+C10*x^2*y^2+C11*x*y^3+C12*y^4+C13*x^5+C14*x ^4*y+ ... C15*x^3*y^2+C16*x^2*y^3+C17*x*y^4+C18*y %stress field sxx = diff(diff(phi,y),y) syy = diff(diff(phi,x),x) sxy = -diff(diff(phi,x),y) %traction force on y=b %Ty on y=b t1 = subs(syy,y,b) t2 = subs(sxy,y,b) %Tx on y=b %traction force on y=-b t3 = subs(syy,y,-b) %Ty on y=-b t4 = subs(sxy,y,-b) %Tx on y=-b %traction force on x=a t5 = subs(sxx,x,a) %Tx on x=a t6 = subs(sxy,x,a) %Ty on x=a %traction force on x=-a t7 = subs(sxx,x,-a) %Tx on x=-a t8 = subs(sxy,x,-a) %Ty on x=-a %find coefficients of polynomials t1,t2,t3,t4 s1 = subs(diff(t1,x,3),x,0)/factorial(3) s2 = subs(diff(t1,x,2),x,0)/factorial(2) s3 = subs(diff(t1,x,1),x,0)/factorial(1) s4 = subs(t1,x,0) s5 s6 s7 s8
= = = =
subs(diff(t2,x,3),x,0)/factorial(3) subs(diff(t2,x,2),x,0)/factorial(2) subs(diff(t2,x,1),x,0)/factorial(1) subs(t2,x,0)
1 of 3
12
4/23 4/23/0 /06 6 8:13 8:13 PM PM
C:\D C:\Doc ocum umen ents ts and and Set Setti ting ngs\ s\We Wei i Cai\ Cai\My My Doc Docum umen ents ts\C \Cou ours rses es.. ...\ .\S5 S522 22a. a.m m
s9 s10 s11 s12
= = = =
subs(diff(t3,x,3) ,x,0)/factorial(3) subs(diff(t3,x,3),x,0)/factorial(3) subs(diff(t3,x,2),x,0)/factorial(2) subs(diff(t3,x,1),x,0)/factorial(1) subs(t3,x,0)
s13 s14 s15 s16
= = = =
subs(diff(t4,x,3),x,0)/factorial(3) subs(diff(t4,x,2),x,0)/factorial(2) subs(diff(t4,x,1),x,0)/factorial(1) subs(t4,x,0)
%The biharmonic equation is 4th order, so applying it to a 5th order polynomial %generates a first order polynomial. % biharm = diff(phi,x,4)+diff(phi,y,4)+2*diff(diff(phi,x,2),y,2) %coefficients of biharm b1 = subs(diff(biharm,x,1),{x,y},{0,0}) b2 = subs(diff(biharm,y,1),{x,y},{0,0}) b3 = subs(biharm,{x,y},{0,0}) %integrated force and torque on x=a Fx1 = simplify(int(t5, y, -b, b)) Fy1 = simplify(int(t6, y, -b, b)) M1 = simplify(int(t5*y simplify(int(t5*y, , y, -b, b)) % %integrated force and torque on x=-a Fx2 = simplify(int(t7, y, -b, b)) Fy2 = simplify(int(t8, y, -b, b)) M2 = simplify(int(t7*y simplify(int(t7*y, , y, -b, b)) %Solve all these equations together s = solve(s1, s2, s3, s4+p, s5, s6, s7, s8, ... s9, s10,s11,s12, s13,s14,s15,s16, ... b1, b2, b3, Fx1,M1, Fy1-p*a, ... 'C1','C2','C3','C4','C5','C6','C7','C8','C9',... 'C10','C11','C12','C13','C14','C15','C16','C17','C18') %construct cell arrays for future use coeffs = {C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,C18} solution = {s.C1,s.C2,s.C3,s.C4,s.C5,s.C6,s.C7,s.C8,s.C9, ... s.C10,s.C11,s.C12,s.C13,s.C14,s.C15,s.C16,s.C17,s.C18} %stress function and stress field solution phi2 = simplify(subs(phi, coeffs, solution)) sxx2 = syy2 = sxy2 =
simplify(subs(sxx , coeffs, solution)) simplify(subs(sxx, simplify(subs(syy, simplify(subs(syy , coeffs, solution)) simplify(subs(sxy, simplify(subs(sxy , coeffs, solution))
%print out results disp('ph disp('phi=') i=') pretty(p pretty(phi2) hi2) disp('sx disp('sxx=') x=') pretty(s pretty(sxx2) xx2)
2 of 3
13
4/23 4/23/0 /06 6 8:13 8:13 PM PM disp('sy disp('syy=') y=') disp('sx disp('sxy=') y=')
C:\D C:\Doc ocum umen ents ts and and Set Setti ting ngs\ s\We Wei i Cai\ Cai\My My Doc Docum umen ents ts\C \Cou ours rses es.. ...\ .\S5 S522 22a. a.m m
pretty(s pretty(syy2) yy2) pretty(s pretty(sxy2) xy2)
%plot stress fields bn = pn = an = xn = [-1:0.05:1]*an yn = [-1:0.1:1]*bn xx = ones(size(yn'))*xn yy = yn'*ones(size(xn)) phin = subs(phi2,{a,b subs(phi2,{a,b,p},{an, ,p},{an,bn,pn}) bn,pn}) sxxn = subs(sxx2,{a,b subs(sxx2,{a,b,p},{an, ,p},{an,bn,pn}) bn,pn}) syyn = subs(syy2,{ subs(syy2,{a,b,p},{ a,b,p},{an,bn,pn} an,bn,pn}) ) sxyn = subs(sxy2,{a,b subs(sxy2,{a,b,p},{an, ,p},{an,bn,pn}) bn,pn})
phin sxxn syyn sxyn
= = = =
subs(phin,{x,y subs(phin,{x,y},{xx,yy} },{xx,yy}) ) subs(sxxn,{x,y subs(sxxn,{x,y},{xx,yy} },{xx,yy}) ) subs(syyn,y subs(syyn,y,yy) ,yy) subs(sxyn,{x,y subs(sxyn,{x,y},{xx,yy} },{xx,yy}) )
figure(1) subplot( subplot(2,2, 2,2,1) 1) mesh(xn, mesh(xn,yn,p yn,phin) hin) title title('\ ('\ph phi') i') xlabe xlabel( l('x' 'x') ) ylabe ylabel(' l('y') y') subplot( subplot(2,2, 2,2,2) 2) mesh(xn, mesh(xn,yn,s yn,sxxn) xxn) title('\ title('\sigm sigma_{x a_{xx}') x}') xlabel(' xlabel('x') x') ylabel(' ylabel('y') y') subplot( subplot(2,2, 2,2,3) 3) mesh(xn, mesh(xn,yn,s yn,syyn) yyn) title('\ title('\sigm sigma_{y a_{yy}') y}') xlabel(' xlabel('x') x') ylabel(' ylabel('y') y') subplot( subplot(2,2, 2,2,4) 4) mesh(xn, mesh(xn,yn,s yn,sxyn) xyn) title('\ title('\sigm sigma_{x a_{xy}') y}') xlabel(' xlabel('x') x') ylabel(' ylabel('y') y') % results %phi= % % 2 3 2 2 3 2 3 2 2 3 5 % p (10 x b + 15 x y b - 2 y b + 5 y a - 5 x y + y ) % - 1/40 ---------------------- ------------------------- ------------% 3 % b %sxx= % % 2 2 2 2 % p y (-6 b + 15 a - 15 x + 10 y ) % - 1/20 ---------------------- ------------% 3 % b %syy= % % 3 2 3 % p (-2 b - 3 y b + y ) % 1/4 --------------------- -% 3 % b %sxy= % % 2 2 % p x (-b + y ) % - 3/4 -------------% 3 % b
3 of 3
σ
φ
xx
5
5
0
0
−5 1
−5 1
2 0 −1
−2
2 0
0 y
14
0 y
x
−1
−2
x
σ
σ
yy
xy
2 0 0
−0.5 −1 1
2 0
0 y
−1
−2
x
−2 1
2
0 y
0 −1
−2
x
Appendix
Potential Field of a Uniformly Charged Ellipsoid Wei Cai Department of Mechanical Engineering, Stanford University, CA 94305-4040 May 28, 2007
Contents 1 Problem Statement
1
2 Orthogonal Curvilinear Coordinates
2
3 Elliptic Coordinates
3
4 Alternative Definition of Elliptic Coordinates
4
5 Ellipsoidal Coordinates
5
6 Ellipsoidal Conductor
8
7 Ellipsoidal Shell
9
8 Uniformly Charged Ellipsoid
11
9 Special Cases
12
10 Application to Hertz Contact Problem
14
A Matlab Files for Analytic Derivation
16
1
Problem Statement
The purpose of this document is to discuss the derivation of a very useful result, which states that the potential field of a uniformly charged ellipsoid is a quadratic function inside the ellipsoid [1]. Specifically, consider an ellipsoid with uniform charge density ρ inside the following
region, x2 y2 z 2 + 2 + 2 =1 a2 b c
(1)
Let V 0 specify the volume of the ellipsoid. The potential field is defined as φ(x)
≡
V 0
ρ
|x − x | dV (x )
(2)
1
where x = (x,y,z) and x = (x , y , z ). If point x is inside the ellipsoid [2], ∞
φ(x,y,z) = π abc ρ
− 1
0
x2 a2 + s
−
y2 b2 + s
−
z2 c2 + s
where ϕ(s)
ds ϕ(s)
≡ (a2 + s)(b2 + s)(c2 + s)
(3)
(4)
If point x is outside the ellipsoid [1], ∞
φ(x,y,z) = π abc ρ
− 1
λ
x2 a2 + s
−
y2 b2 + s
−
z2 c2 + s
ds ϕ(s)
(5)
where λ is the greatest root of the equation f (s) = 0, where 2
f (s)
2
2
≡ a2x+ s + b2y+ s + c2z+ s − 1
(6)
The physical significance of function f (s) is that, for each s > 0, the equation f (s) = 0 defines an ellipsoid (larger than the original ellipsoid), which is an isosurface of the potential field φ(x) generated by the original ellipsoid (with uniform charge density). Therefore, all the value of s [0, ) corresponds to a family of ellipsoids, called the confocal family ; the potential field is a constant on each ellipsoid in the family.
∈
∞
Physical significance: As a consequence of the above result, the second spatial derivatives of the
potential field is a uniform second order tensor inside the ellipsoid. This is closely analogous to Eshelby’s Theorem , which states that the stress field inside an ellipsoidal inclusion is uniform [3]. Landau and Lifshitz also used the above result to show that normal stress distribution inside the contact area between two smooth elastic media (Hertzian contact problem) has the ellipsoidal shape [2].
2
Orthogonal Curvilinear Coordinates
The proof of this result is best discussed using ellipsoidal coordinates , which is an orthogonal curvilinear coordinate system. In this section, we first summarize the major results concerning orthogonal curvilinear coordinates. We will then introduce the elliptic coordinates (2D) and ellipsoidal coordinates (3D) in the following sections. Let the Cartesian coordinates be specified by (x1 , x2 , x3 ) = (x,y,z). An arbitrary differential length in space ds is specified by (ds)2 = (x1 )2 + (x2 )2 + (x3 )2 . Now consider a general orthogonal curvilinear coordinate system, (q 1 , q 2 , q 3 ), which are related to the Cartesian coordinates by q m = q m (x1 , x2 , x3 ),
xm = x m (q 1 , q 2 , q 3 )
(7)
An arbitrary differential length in space can be expressed by (ds)2 = (h1 dq 1 )2 + (h2 dq 2 )2 + (h3 dq 3 )2
(8)
2
x = const
y = c o ns t
! = c o ns t
" = c o ns t
x=
x=
-a
a
Figure 1: Contour lines in Cartesian and elliptic coordinates [5]. where hi are called scale factors and have the following expressions [4]. (h1 )2 = (h2 )2 = (h3 )2 =
∂x k ∂x k ∂q 1 ∂q 1 ∂x k ∂x k ∂q 2 ∂q 2 ∂x k ∂x k ∂q 3 ∂q 3
(9)
The index notation is used here and k is a dummy index that is summed from 1 to 3. Let ek be ˆk be the basis vectors of the the basis vectors of the Cartesian coordinates (x1 , x2 , x3 ) and let e curvilinear coordinates (q 1 , q 2 , q 3 ). The gradient of a scalar field φ(q 1 , q 2 , q 3 ) is, ∂φ 1 ∂φ 1 ∂φ ∇φ = ˆe1 h11 ∂q + ˆ + ˆe3 e2 h2 ∂q 2 h3 ∂q 3 1
(10)
The Laplacian of a scalar field φ(q 1 , q 2 , q 3 ) is,
∇
1 ∂ 2 φ = h1 h2 h3 ∂q 1
h2 h3 ∂φ h1 ∂q 1
∂ + ∂q 2
h3 h1 ∂φ h2 ∂q 2
∂ + ∂q 3
h1 h2 ∂φ h3 ∂q 3
(11)
In 2-dimension, the Laplacian of a scalar field φ(q 1 , q 2 ) reduces to the following.
∇ 3
2
1 ∂ φ = h1 h2 ∂q 1
h2 ∂φ h1 ∂q 1
∂ + ∂q 2
h1 ∂φ h2 ∂q 2
Elliptic Coordinates
The most common definition of elliptic coordinate (µ, ν ) is [5], x = a cosh µ cos ν y = a sinh µ sin ν 3
(12)
With this definition, we can show that x2 x2 + = cos2 ν + sin2 ν = 1 2 2 2 2 a cosh µ a sinh µ x2 x2 = cosh2 µ sinh2 µ = 1 2 2 2 2 a cos ν a sin ν
−
(13)
−
Therefore, the contour lines of µ = const form a set of ellipses, and the contour lines of ν = const form a set of hyperbolas, as shown in Fig. 1. This figure also shows that in the limit of a 0, or when the distance from the origin is much greater than a, the elliptic coordinates becomes very close to polar coordinates.
→
The Jacobian matrix between Cartesian coordinates (x, y) and elliptic coordinates (µ, ν ) is J
≡
∂x ∂µ ∂x ∂µ
∂x ∂ν ∂x ∂ν
=
a sinh µ cos ν a cosh µ sin ν
−a cosh µ sin ν a sinh µ cos ν
(14)
The elliptic coordinates is an orthogonal coordinate system because the two columns of matrix J are orthogonal to each other. The scale factors are hµ = hν =
∂x ∂µ
2
∂x ∂ν
2
2
+
∂y ∂µ
2
+
∂y ∂ν
(15)
It is easy to show that hµ = hν = a
sinh2 µ + sin2 ν =
det(J )
(16)
Therefore, the Laplacian of a scalar field φ is 2
∇φ
= =
4
∂ 2 ∂ 2 + φ(x, y) ∂x 2 ∂y 2 1 ∂ 2 ∂ 2 + φ(µ, ν ) a2 (sinh2 µ + sin2 ν ) ∂µ 2 ∂ν 2
(17)
Alternative Definition of Elliptic Coordinates
An alternative definition of elliptic coordinates makes it more natural to generalize the concept to ellipsoidal coordinates in 3D. Consider an ellipse x2 y2 + 2 =1 a2 b
(18)
We will assume a > b without loss of generality. Now consider a family of curves defined by f (s) When s >
≡
x2 y2 + a2 + s b2 + s
− 1 = 0
(19)
−b2, it defines an ellipse. When −b2 > s > −a2, it defines a hyperbola. 4
For a given (x, y), let (µ, ν ) be the two largest roots of the equation f (s) = 0. There is a one-to-one correspondence between (x, y) and (µ, ν ), if we assume x > 0, y > 0. Specifically, (a2 + µ)(a2 + ν ) a2 b2 (b2 + µ)(b2 + ν ) b2 a2
x2 =
−
y2 =
−
(20)
This relationship can be verified by plugging it into the definition of f (s), x2 y2 a2 + ν b2 + ν f (µ) = + 1= 2 + 1=0 a2 + µ b2 + µ a b2 b2 a2 x2 y2 a2 + µ b2 + µ f (ν ) = + 1= 2 + 1=0 a2 + ν b2 + ν a b2 b2 a2 The four components of the Jacobian matrix can be obtained.
−
−
−
−
∂x ∂µ
1 a2 + ν = 2 (a2 + µ)(a2 b2 )
∂y ∂µ
1 b2 + ν = 2 (b2 + µ)(b2 a2 )
−
−
1/2
1/2
− − − −
∂x 1 a2 + µ = ∂ν 2 (a2 + ν )(a2
(21)
1/2
− b2)
∂y 1 b2 + µ = ∂ν 2 (b2 + ν )(b2 a2 )
−
The (µ, ν ) coordinate system is orthogonal because
1/2
(22)
∂x ∂x ∂y ∂y + =0 ∂µ ∂ν ∂µ ∂ν
(23)
≡ (a2 + s)(b2 + s). The scale factors can be expressed as µ − ν
Define function ϕ(s) hµ = hν =
1 2 1 2
−
ϕ(µ) ν µ ϕ(ν )
(24)
Therefore, the Laplacian of a scalar field φ(µ, ν ) is 2
∇φ
−
=
1 ∂ hµ hν ∂µ
=
4 ϕ(µ)ϕ(ν ) ∂ µ ν ∂µ
=
4
µ
− ν
hν ∂φ hµ ∂µ
ϕ(µ)
∂ ∂µ
+
∂ ∂ν
hµ ∂φ hν ∂ν
ϕ(µ) ∂φ ϕ(ν ) ∂µ
+
∂φ ∂µ
+
ϕ(µ)
∂ ∂ν
ϕ(ν )
ϕ(µ) ∂φ ϕ(ν ) ∂ν
∂ ∂ν
ϕ(ν )
∂φ ∂ν
(25)
For derivation details see elliptic coord.m.
5
Ellipsoidal Coordinates
Generalizing the elliptic coordinates defined above, we obtain the ellipsoidal coordinates [6]. Consider an ellipsoid, Consider an ellipse x2 y2 z 2 + 2 + 2 =1 a2 b c
(26) 5
! =
c o n st
" =
# =
c o n st
c o n st
Figure 2: Isosurfaces of ellipsoidal coordinates [7]. We will assume a > b > c without loss of generality. Now consider a family of curves defined by 2
f (s)
2
2
≡ a2x+ s + b2y+ s + c2z+ s − 1 = 0
(27)
For λ > c2 , f (λ) = 0 defines an ellipsoid. When c2 > µ > b2 , f (µ) = 0 defines a one-sheet hyperbola. When b2 > ν > a2 , f (ν ) = 0 defines a two-sheet hyperbola, as shown in Fig. 2. For a given (x,y,z), let (λ,µ,ν ) be the three largest roots of equation f (s) = 0. There is a one-to-one correspondence between (x,y,z) and (λ,µ,ν ), if we assume x > 0, y > 0, z > 0.
−
2
x
=
y2 = z2 =
−
−
−
−
(a2 + λ)(a2 + µ)(a2 + ν ) (a2 b2 )(a2 c2 ) (b2 + λ)(b2 + µ)(b2 + ν ) (b2 a2 )(b2 c2 ) (c2 + λ)(c2 + µ)(c2 + ν ) (c2 a2 )(c2 b2 )
−
−
−
−
−
−
(28)
(29)
The following limit applies, λ>
−c2 > µ > −b2 > ν > −a2
(30)
The nine components of the Jacobian matrix can be obtained.
6
∂x ∂λ
1 (a2 + µ)(a2 + ν ) = 2 (a2 + λ)(a2 b2 )(a2 c2 )
1/2
∂x ∂µ
1 (a2 + λ)(a2 + ν ) = 2 (a2 + µ)(a2 b2 )(a2 c2 )
1/2
∂x 1 (a2 + λ)(a2 + µ) = ∂ν 2 (a2 + ν )(a2 b2 )(a2 c2 )
1/2
∂y ∂λ
1 (b2 + µ)(b2 + ν ) = 2 (b2 + λ)(b2 a2 )(b2 c2 )
1/2
∂y ∂µ
1 (b2 + λ)(b2 + ν ) = 2 (b2 + µ)(b2 a2 )(b2 c2 )
−
−
−
−
−
−
−
−
−
−
−
−
∂z ∂λ
1 (c2 + µ)(c2 + ν ) = 2 (c2 + λ)(c2 a2 )(c2 b2 )
∂z ∂µ
1 (c2 + λ)(c2 + ν ) = 2 (c2 + µ)(c2 a2 )(c2 b2 )
−
−
−
−
−
(31)
1/2
1/2
∂y 1 (b2 + λ)(b2 + µ) = ∂ν 2 (b2 + ν )(b2 a2 )(b2 c2 )
−
(32)
1/2
1/2
1/2
∂z 1 (c2 + λ)(c2 + µ) = ∂ν 2 (c2 + ν )(c2 a2 )(c2 b2 ) The (λ,µ,ν ) coordinate system is orthogonal because ∂x ∂x ∂y ∂y ∂z ∂z + + = 0 ∂λ ∂µ ∂λ ∂µ ∂λ ∂µ ∂x ∂x ∂y ∂y ∂z ∂z + + = 0 ∂µ ∂ν ∂µ ∂ν ∂µ ∂ν ∂x ∂x ∂y ∂y ∂z ∂z + + = 0 ∂λ ∂ν ∂λ ∂ν ∂λ ∂ν Define function ϕ(s) (a2 + s)(b2 + s)(c2 + s). The scale factors can be expressed as hλ =
1 2
hµ =
1 2
hν =
1 2
≡− − − (λ
µ)(λ ϕ(λ)
− ν )
(µ
λ)(µ ϕ(µ)
− ν )
(ν
λ)(ν ϕ(ν )
− µ)
The Laplacian of a scalar field φ(λ,µ,ν ) is, 1 ∂ hµ hν ∂φ ∂ 2 φ = + hλ hµ hν ∂λ hλ ∂λ ∂µ
∇
=
4 ϕ(λ) ∂ (λ µ)(λ ν ) ∂λ +
4 ϕ(ν ) ∂ (ν λ)(ν µ) ∂ν
−
−
ϕ(λ)
∂φ ∂λ
ϕ(ν )
+
∂φ ∂ν
7
+
(34)
(35)
− − − − hν hλ ∂φ hµ ∂µ
(33)
∂ ∂ν
hλ hµ ∂φ hν ∂ν
4 ϕ(µ) ∂ (µ λ)(µ ν ) ∂µ
ϕ(µ)
∂φ ∂µ
(36)
For derivation details see ellipsoidal coord.m.
6
Ellipsoidal Conductor
The discussion from here on follows Kellogg’s book closely [1], with some variables renamed to follow the notation here. Suppose we would like to solve for the potential function in space produced by an ellipsoidal conductor that contains surface charges [1]. For a perfect conductor, the potential on its surface (as well as the interior) is a constant. Therefore, we are trying to solve the Poisson’s equation,
∇2φ(x) = 0
(37)
subject to the boundary condition that φ(x) = φ 0 when point x is on the ellipsoidal surface, x2 y2 z 2 + 2 + 2 =1 a2 b c
(38)
and that φ(x) = 0 as x . Introducing the ellipsoidal coordinates (λ,µ,ν ) as defined in the previous section. The surface of the (original) ellipsoid is simply the isosurface of λ = 0. The limit of x corresponds to the limit of λ . Therefore, when φ is expressed in term of the ellipsoidal coordinates, i.e. φ(λ,µ,ν ), the boundary condition is very simple,
| | → ∞
| | → ∞
→ ∞
φ(λ = 0, µ , ν ) = φ0 φ(λ
→ ∞, µ , ν)
= 0
(39)
Notice that the Laplacian in the Possion’s equation ( 2 φ = 0) in the elliptical coordinates is defined in Eq. (36). A natural trial solution is a function φ(λ) that only depends on λ, but not on µ or ν . In this case,
∇
∇
4 ϕ(λ) ∂ 2 φ = (λ µ)(λ ν ) ∂λ
−
ϕ(λ)
∂φ ∂λ ∂φ ∂λ
= =
−
ϕ(λ)
∂φ ∂λ
=0
(40)
− E 2 − 2 E ϕ(λ)
∞
φ(λ) =
λ
E ds 2 ϕ(s)
(41)
where E is a constant. The potential field inside the conductor is a constant and equals to the potential on the surface (λ = 0), which is, ∞
φ0 =
E ds 2 ϕ(s)
0
(42)
The surface charge of the conductor σ(x) can be obtained from the following relationship. ∂φ(x) = ∂n +
−4πσ(x)
(43) 8
where
∂ ∂n +
is the gradient along the surface normal. x) − 4π1 ∂φ( ∂n +
σ = =
−
1 4π
Notice that at λ = 0, h λ = 1 4π
σ =
1 ∂φ(λ) hλ ∂λ
(44)
λ=0
µν/ϕ(λ)/2. Therefore,
2 ϕ(λ) E E = µν 2ϕ(λ) 4π µν
√
√
(45)
On the surface of the ellipsoid, λ = 0, µν = a 2 b2 c2
x2 y 2 z 2 + 4 + 4 a4 b c
(46)
The equation of the plane tangent to ellipsoid at point (x,y,z) is (X
− x) xa2 + (Y − y) yb2 + (Z − z) zc2 = 0
(47)
The surface normal of the tangent plane is
x2 y 2 z 2 + 4 + 4 a4 b c
x y z , , / n = a2 b2 c2
(48)
The shortest distance from the origin to this plane is p =
x2 a2
+
x2 a4
y2 b2
+
+
y2 b4
z2 c2
+
z2 c4
=
1
x2 a4
+
y2 b4
+
z2 c4
=
abc √ µν
(49)
Therefore, the surface charge density is σ =
E E p = 4π abc 4π abc
1
x2 a4
+
y2 b4
+
(50)
z2 c4
This result is related to the problem of a uniformly charged ellipsoid. As will be shown in the following section, the above expression of σ is exactly the amount of charge contained in a thin shell between two similar ellipsoids, in the limit of shell thickness going to zero.
7
Ellipsoidal Shell
Consider a set of similar ellipsoids, x2 y2 z2 x2 y 2 z 2 + + = 1 , equivalently 2 + 2 + 2 = u 2 2 2 2 (au) (bu) (cu) a b c
(51)
whose semi-axes are au, bu and cu. They are simply the original ellipsoid scaled by a factor u in all three directions. Notice that this family of ellipsoids are different from the family of ellipsoids defined by f (λ) = 0 (whose shapes are not similar to each other). For 0 < u < 1, these ellipsoids 9
are smaller than the original ellipsoid (while for 0 < λ < the ellipsoids defined by f (λ) = 0 are all larger than the original ellipsoid). Consider an ellipsoidal shell contained between two ellipsoids defined by u1 and u 2 = u 1 + ∆u. Let the volume density of the charge distribution inside the shell to be a constant ρ. In the limit of ∆u 0, the shell reduces to a surface with a surface charge density σ. Obviously, the surface charge density is proportional to the local thickness of the shell, ∆h, i.e.,
∞
→
σ = ρ ∆h
(52)
Let (ux,uy,uz) be a point on the surface of the ellipsoid defined by u. Let p(x,y,z,u) be the shortest distance from the origin to the plane tangent to the ellipsoid at point (ux,uy,uz). Then, u p(x,y,z,u) = 2 y2 x2 + + zc4 4 a b4 ∆h(x,y,z) = σ =
x2 a4
x2 a4
∆u +
y2 b4
+
z2 c4
+
z2 c4
ρ ∆u +
y2 b4
(53)
Compare this with Eq. (50), we can conclude that the surface charge of the shell is the same as the equilibrium surface charge of a ellipsoidal conductor. The correspondence is made complete if we set u1 = 1, u 2 = 1 + ∆u, and E = ρ∆u 4πabc E = 4πabcρ∆u (54) This means that the potential field produced by this ellipsoidal shell (u1 = 1, u 2 = 1 + ∆u) is ∞
φ(x) = φ(λ) = 2πabcρ ∆u
ds ϕ(s)
λ
(55)
The potential field inside the ellipsoidal shell is a constant and equals to the potential on the surface (λ = 0), which is, ∞
φ0 = 2πabcρ ∆u
ds ϕ(s)
0
(56)
In summary, Eq. (55) describes the potential generated by an ellipsoidal shell with uniform density ρ, whose boundary is the original ellipsoid and a similar ellipsoid scaled by a factor (1 + ∆u). This result can be generalized to an ellipsoidal shell between u1 = u and u2 = u + ∆u for arbitrary u. The potential field at a point x = (x,y,z) outside this shell is, 2
φ(x) = 2πabcρu ∆u
∞
ds ϕ(u, s)
λ(u)
(57)
where λ(u) is the greatest root of equation f (u, s) = 0 for given (x,y,z) and u. f (u, s) and ϕ(u, s) are generalization of the previously defined functions f (s) and ϕ(s). f (u, s) ϕ(u, s)
≡ ≡
x2 y2 z2 + + a2 u2 + λ b2 u2 + λ c2 u2 + λ (a2 u2 + s)(b2 u2 + s)(c2 u2 + s)
− 1
(58) (59)
The factor u2 in Eq. (57) accounts for the fact that surface area of the scaled shell and hence its total charge content is u 2 times those of the original shell (u = 1). 10
8
Unif Unifor ormly mly Char Charge ged d Ell Ellip ipso soid id
A uniformly charged ellipsoid can be considered as a collection of many layers of ellipsoid shells considered above. For a point x outside the ellipsoid, its potential value should be an integral of u u from 0 to 1, 1
∞
ds du ϕ(u, s)
≡ ≡ U e (x) = 2πabcρ
u2
0
λ(u)
λ(u)/u2 and t and t
Define new variables v
1
U e (x) = 2πabcρ
∞
0
v
Perform integration by parts on 1
∞
u
0
v
(60)
s/u2 . Then φ Then φ((u, s) = u 2 φ(t).
dt du ϕ(t)
u
(61)
du, du,
dt u2 du = du = 2 ϕ(t)
∞
dt ϕ(t)
v
1
+
0
1 2
1
1 dv du ϕ(v ) du
u2
0
(62)
Notice that v is the greatest root of equation, x2 y2 z2 + + = u 2 a2 + v b2 + v c2 + v For u or u = = 1, v = λ = λ,, while for u 1
(63)
→ 0, v 0, v → ∞. Hence
− u2 2
1 2
∞
v
1
0
u2
dt ϕ(t)
=
0
1 2
∞
dt ϕ(t)
λ
1 dv du = ϕ(v ) du =
λ
1 2
u2
∞
1 2
∞
λ
(64)
1 dv ϕ(v ) x2 y2 z2 + + a2 + v b2 + v c2 + v
Therefore, the potential outside the ellipsoid is ∞
U e (x) = π abc ρ
− 1
λ
x2 y2 z2 + + a2 + v b2 + v c2 + v
where λ is the greatest root of equation
1 dv ϕ(v )
x2 y2 z2 + + =1 a2 + λ b2 + λ c2 + λ
1 dv ϕ(v )
(65)
(66)
(67)
as previously defined. Let us now consider the potential field at a point x = (x,y,z) x,y,z) inside the ellipsoid. Here we have to cut the ellipsoid into two parts. Point x is on the outside of part 1, but is on the inside of part 2. Let u Let u 0 correspond to the ellipsoid that pass through point x , i.e. x2 y2 z 2 + 2 + 2 = u 20 2 a b c
11
(68)
Therefore, part 1 contains the ellipsoidal shells with 0 < 0 < u < u0 and part 2 contains the ellipsoidal shells with u0 < u < 1. The potent potential ial at point x from an ellipsoidal shell in part 1 has the same expres expressio sion n as the above above.. But potential potential at point point x from an ellipsoidal shell in part 2 equals to the potential value at the shell surface (because point x is inside inside the shell). Therefore, Therefore, the total potential at an interior point x is,
u0
U i (x) = 2πabcρ 2 πabcρ
∞
u
0
1
dt du + ϕ(t)
v
∞
dt du ϕ(t)
u
u0
0
(69)
Perform integration by parts on the first integral, u0
∞
u
0
v
Notice that when u = 0, v 0, v = u0
∞
u
0
v
u2
dt du = du = ϕ(t)
2
u0
− − ∞
v
dt ϕ(t)
0
u0
1 + 2
1 dv du ϕ(v ) du
u2
0
(70)
∞, and when u = u = u 0 , v = 0. Hence
dt du = ϕ(t)
u20 2
=
u20 2
∞
0
∞
0
dt ϕ(t) dt ϕ(t)
1 2 1 2
∞
0
∞
0
1 dv ϕ(v ) x2 y2 z2 + + a2 + v b2 + v c2 + v
u2
1 dv ϕ(v )
The second integral in Eq. (69) can be carried out because the inner integral is a constant. 1
dt 1 u20 du = du = 2 ϕ(t)
∞
−
u
u0
0
∞
0
dt ϕ(t)
(71)
Therefore, the total potential at an interior point x is, i s, ∞
U i (x) = π abc ρ
− 1
0
x2 y2 z2 + + a2 + v b2 + v c2 + v
1 dv ϕ(v )
(72)
Notice that the only difference between Eq. (66) and Eq. (72) is in the lower limit of integration (λ versus versus 0). Because Because the range of integration integration is constant, constant, the potential field inside the ellipsoid, U i (x), is simply a quadratic function of space.
9
Spec Specia iall Case Casess
9.1
Unifor Uniformly mly Charged Charged Sphere Sphere
The situation of a = a = b b = = c c describes describes a sphere with radius a radius a.. When the sphere has a uniform charge density ρ density ρ,, the potential distribution inside the sphere is, φ(x,y,z) x,y,z) = π a3 ρ (A
− B x2 − B y2 − B z2) = π a3 ρ (A − B r2)
(73)
where ∞
A =
0
∞
B =
0
ds 2 = a (a2 + s)3/2 ds 2 = 3a3 (a2 + s)5/2
(74)
12
(75)
Therefore φ(r ) =
2πρ (3a (3a2 3
− r2 )
(76)
The potential value on the sphere surface is φ(r = a = a)) =
2πρ (3a (3a2 3
− a2) = 4πρ a2 3
(77)
Notice that Q = 4π a3 ρ/3 ρ/3 is the total charge contained in the sphere. Hence φ(r = a = a)) =
Q a
(78)
This is equivalen equivalentt to the potential potential produced by a point charge charge at origin. origin. This confirms confirms Newton’s Theorem, which states that the potential field outside a uniformly charged sphere is equivalent to that produced by a p oint charge charge located at the center center of the sphere. sphere. The potential potential field outside outside the sphere (r (r > a) is φ(r ) = 9.2 9.2
Q r
(79)
Charg Charged ed Meta Metall Disc Disc
Consider the case of a = b and c 0. In this limit, limit, the region region inside inside the ellips ellipsoid oid reduces reduces to a circle of radius a in the x the x--y plane (z (z = 0). The potential distribution inside this area is,
→
φ(x, y) =
(x
−
ρ dx dy dz x )2 + (y (y y )2 + z2
−
Because the integration limit for z for z is φ(x, y) = 2ρ c
± c 1 − (x /a) /a)2 − (y /a) /a)2 ,
x 2 +y 2 ≤a2
− (x
dx dy x )2 + (y (y
− y )2
− − x a
1
2
y a
2
(80)
Using the results obtained above, the potential must be a quadratic function inside the area of radius a, φ(x, y) = π a2 c ρ (A B x2 B y2 ) ∞ ds π A = = (a2 + s) s a 0 ∞ ds π B = = 2a3 (a2 + s)2 s 0
−
−
√
√
(81)
(82)
In other words, we have obtained the following relationship,
x 2 +y 2 ≤a2
where r =
− (x
dx dy x )2 + (y (y
2
−y )
− − 1
x a
2
y a
2
=
π 2 a (A 2
− B x2 − B y2)
=
π2 (2a (2a2 4a
− r2 )
(83)
x2 + y2 . The potential value is maximum at the circumference of the circle, r = a = a,,
φ(r = a = a)) = 2ρ c
π2 2 π Q a = 4a 3a
13
(84)
where Q = 34 π a2 c ρ is the total charge in the ellipsoid. The results obtained here can be used to model a very thin circular metal (conductor) plate with radius a. The total charge in the metal plate is Q, whereas the equilibrium charge distribution on the surface is 3Q σ(x, y) = 2π a2
− − x a
1
2
y a
2
(85)
The potential on the metal surface is a constant φ0 =
10
πQ 3a
(86)
Application to Hertz Contact Problem
For simplicity, let us consider a rigid sphere of radius R indenting an elastic half space. The discussion here follows closely that in Landau and Lifshitz [2]. Choose the coordinate system such that the elastic half space occupies the z 0 domain and the z = 0 plane is the surface of the half space. The Boussinesq solution tells us about the surface displacement of the elastic half space in response to a point force of magnitude F acting at the origin in the z direction.
≤
−
uz =
−
F (1 ν 2 ) 1 π E r
−
(87)
where E = 2µ(1 + ν ) is the Young’s modulus of the elastic half space. The shape of the indentor can be described by function u0 (x, y) =
x2 y2 + 2R 2R
(88)
Let d be the indentation depth. Hence inside the region of contact (S ), the surface displacement of the half space is uz (x, y) =
−d + u0(x, y) = −
x2 y2 d+ + 2R 2R
(89)
Let pz (x, y) be the normal stress on the surface of the half space. Inside the region of contact, pz (x, y) < 0; outside the region of contact, p z (x, y) = 0. The total indenting force F is, F =
−
pz (x, y) dxdy
(90)
S
Using the Boussinesq solution, the surface stress p z (x, y) and the surface displacement u z (x, y) are related by, uz (x, y) =
1
− ν 2
π E
pz (x , y ) (x x )2 + (y
S
−
− y )2
dx dy
(91)
Therefore, our task is to find a function p z (x, y) that satisfies the following condition, pz (x , y ) (x x )2 + (y
S
−
π E dx dy = 1 ν 2 y )2
−
−
−
14
x2 y2 d+ + 2R 2R
(92)
By symmetry, we expect the contact area to be a circle, and let a be its radius. Motivated by Eq. (83), we postulate the following form for pz (x, y), pz (x, y) = p0
−
− − 2
x a
1
y a
2
(93)
The constant p0 is related with F by F = p0
− − x a
1
S
p0
2
y a
2
dxdy
2 π a2 = p0 3 3 F = 2 π a2
(94)
Plug the expression of p z (x, y) into Eq. (92), and using Eq. (83), we have π2 p0 (2a2 r 2 ) = 4a 3 F π 2 (2a2 r 2 ) = 2 2π a 4a 3 (1 ν 2 ) F ( 2a2 + r 2 ) = 3 8Ea
−
−
− −
−
−
π E 1 ν 2 π E 1 ν 2
− −
r2 d+ 2R r2 d+ 2R
− − 2
r −d + 2R
(95)
Therefore, 1 R
=
a =
3 (1 ν 2 ) F 4 E a3 3 (1 ν 2 ) F R 4 E
−
−
1/3
(96)
The indentation depth is d =
3 (1 ν 2 ) F 2 3 (1 ν 2 ) F 3 (1 ν 2 ) 2a = = 8 E a3 4Ea 4 E
−
−
−
2/3
F 2/3 R1/3
(97)
Acknowledgement
I want to thank Prof. David Barnett for useful discussions and lending me Kellogg’s book to read.
15
A
Matlab Files for Analytic Derivation
% File: elliptic_coord.m % Purpose: analytic derivation of the properties of elliptic coordinates syms x y a b mu nu x = sqrt( (a^2+mu)*(a^2+nu)/(a^2-b^2) ); y = sqrt( (b^2+mu)*(b^2+nu)/(b^2-a^2) ); dxdmu dxdnu dydmu dydnu
= = = =
simplify(diff(x,mu)); simplify(diff(x,nu)); simplify(diff(y,mu)); simplify(diff(y,nu));
disp(’dxdmu=’); disp(’dxdnu=’); disp(’dydmu=’); disp(’dydnu=’);
pretty(dxdmu); pretty(dxdnu); pretty(dydmu); pretty(dydnu);
%check orthogonality simplify(dxdmu*dxdnu + dydmu*dydnu) h_mu = simplify( sqrt( (dxdmu)^2 + (dydmu)^2 ) ); h_nu = simplify( sqrt( (dxdnu)^2 + (dydnu)^2 ) ); disp(’h_mu=’); pretty(h_mu); disp(’h_nu=’); pretty(h_nu);
16
% File: ellipsoidal_coord.m % Purpose: analytic derivation of the properties of ellipsoidal coordinates syms x y z a b c lm mu nu x = sqrt( (a^2+lm)*(a^2+mu)*(a^2+nu)/(a^2-b^2)/(a^2-c^2) ); y = sqrt( (b^2+lm)*(b^2+mu)*(b^2+nu)/(b^2-a^2)/(b^2-c^2) ); z = sqrt( (c^2+lm)*(c^2+mu)*(c^2+nu)/(c^2-a^2)/(c^2-b^2) ); dxdlm dxdmu dxdnu dydlm dydmu dydnu dzdlm dzdmu dzdnu
= = = = = = = = =
simplify(diff(x,lm)); simplify(diff(x,mu)); simplify(diff(x,nu)); simplify(diff(y,lm)); simplify(diff(y,mu)); simplify(diff(y,nu)); simplify(diff(z,lm)); simplify(diff(z,mu)); simplify(diff(z,nu));
disp(’dxdlm=’); disp(’dxdmu=’); disp(’dxdnu=’); disp(’dydlm=’); disp(’dydmu=’); disp(’dydnu=’); disp(’dzdlm=’); disp(’dzdmu=’); disp(’dzdnu=’);
pretty(dxdlm); pretty(dxdmu); pretty(dxdnu); pretty(dydlm); pretty(dydmu); pretty(dydnu); pretty(dzdlm); pretty(dzdmu); pretty(dzdnu);
%check orthogonality [ simplify(dxdlm*dxdmu + dydlm*dydmu + dzdlm*dzdmu) simplify(dxdmu*dxdnu + dydmu*dydnu + dzdmu*dzdnu) simplify(dxdlm*dxdnu + dydlm*dydnu + dzdlm*dzdnu) ] h_lm = simplify( sqrt( (dxdlm)^2 + (dydlm)^2 + (dzdlm)^2 ) ); h_mu = simplify( sqrt( (dxdmu)^2 + (dydmu)^2 + (dzdmu)^2 ) ); h_nu = simplify( sqrt( (dxdnu)^2 + (dydnu)^2 + (dzdnu)^2 ) ); disp(’h_lm=’); pretty(h_lm); disp(’h_mu=’); pretty(h_mu); disp(’h_nu=’); pretty(h_nu);
17
