Projekt 6: modelowanie pola elektrycznego i magnetycznego w otoczeniu obracającej się naładowanej sfery

(1)

Projekt 6: modelowanie pola elektrycznego i magnetycznego w otoczeniu obracającej się naładowanej sfery

10 czerwca 2021

1 Wstęp

Na powierzchni sfery o promieniu R znajduje się ładunek o gęstości powierzchniowej σ. Sfera obraca się wokół osi leżącej w płaszczyźnie xy (rys.1). Na zajęciach wyznaczymy rozkład potencjału i pola elektrycznego oraz pola magnetycznego w jej otoczeniu.

Rysunek 1: Środek naładowanej sfery znajduje się w początku układu współrzędnych. Sfera obraca się wokół osi leżącej w płaszczyźnie xy.

1.1 Pole elektryczne

Skorzystamy z zasady superpozycji

V (⃗ r) = 1 4πε ₀

Z ρ(⃗ r ^′ )

|⃗r − ⃗r ^′ | d ³ r ^′ (1)

Gęstość ładunku deﬁnujemy w postaci ρ(⃗ r ^′ ) = σ · δ(r ^′ − R) (σ = const) i wprowadzamy współrzędne sferyczne

V (⃗ r) = 1 4πε 0

Z π 0

dθ ^′ Z 2π 0

dϕ ^′ Z ∞ 0

dr ^′ sinθ ^′ r ^′2 σδ(r ^′ − R)

|⃗r − ⃗r ^′ | (2)

= 1

4πε ₀ Z π 0

dθ ^′ Z 2π 0

dϕ ^′ sinθ ^′ R ² σ

|⃗r − ⃗R ^′ | (3)

(2)

gdzie:

⃗ r →

 



 

x = r sin θ cos φ y = r sin θ sin φ z = r cos θ

R ⃗ ^′ →

 



 

R ^′ _x = R sin θ ^′ cos φ ^′ R ^′ _y = R sin θ ^′ sin φ ^′ R ^′ _z = R cos θ ^′

(4)

Całkę (3) wyznaczymy numerycznie stosując złożenie dwóch kwadratur (2 wymiary) i wzór całkowy parabol (Simpsona). W tym celu na powierzchnii sfery wprowadzamy siatkę węzłów

∆θ = π

N (5)

∆φ = 2π

M (6)

θ → θ i = ∆θ · i, i = 0, 1, . . . , N (7)

φ → φ j = ∆φ · j, j = 0, 1, . . . , M (8)

a całkowanie zastępujemy sumowaniem V (⃗ r) = σR ²

4πε 0

∆θ 3

∆φ 3

X N i=0

X M j=0

a _i b _j sin θ ^′ _i

|⃗r − ⃗R

^′

_ij | (9)

Uwaga: ⃗ R

^′

_i,j oznacza że współrzędne R ^′ _x , R ^′ _y , R ^′ _z zależą od pary indeksów (i, j) ↔ (θ _i ^′ , φ ^′ _j ) zgodnie z wzorami (4).

Współczynniki a _i /b _j są określone następująco

a i =

 



 

a _i = 1, i = 0 ∨ i = N

a i = 4, i mod 2 = 1 (i nieparzyste) a i = 2, i mod 2 = 0 (i parzyste)

b j =

 



 

b _j = 1, j = 0 ∨ j = M

b j = 4, j mod 2 = 1 (j nieparzyste) b j = 2, j mod 2 = 0 (j parzyste)

(10) Podobnie możemy wyznaczyć numerycznie wektor pola elektrycznego

E(⃗ ⃗ r) = −∇V (⃗r) = 1 4πε 0

Z π 0

dθ ^′ Z 2π 0

dϕ ^′ R ² σ sin θ ^′ (⃗ r − ⃗R ^′ )

|⃗r − ⃗R ^′ | ³ (11)

po zastąpieniu całki sumą

E(⃗ ⃗ r) = σR ² 4πε 0

∆θ 3

∆φ 3

X N i=0

X M j=0

a i b j

sin θ _i ^′ (⃗ r − ⃗R ij

^′

)

|⃗r − ⃗R

^′

_ij | ³ (12)

1.2 Pole magnetyczne

Na powierzchnii sfery rozłożony jest ładunek, gdy sfera obraca, generuje się prąd, który jest źródłem pola magnetycznego. Indukcję pola magnetycznego także wyznaczymy z zasady superpozycji

B(⃗ ⃗ r) = − µ ₀ 4π

Z (⃗ r − ⃗r ^′ ) ×⃗j(⃗r ^′ )

|⃗r − ⃗r ^′ | ³ d ³ r ^′ (13)

Gęstość prądu deﬁniujemy jako

⃗j(⃗r ^′ ) = ρ(⃗ r ^′ )⃗ v ^′ = σδ(r ^′ − R) ⃗ω × ⃗r ^′ (14) wstawiamy ją do równania (13) a następnie wykonujemy całkowanie po zmiennej radialnej

B(⃗ ⃗ r) = − µ 4π

Z π 0

dθ ^′ Z 2π 0

dϕ ^′ σR ² sin θ ^′ (⃗ r − ⃗R ^′ ) × (⃗ω × ⃗r ^′ )

|⃗r − ⃗R ^′ | ³ (15)

Całkę liczymy numerycznie (analogicznie jak dla pola elektrycznego)

B p (⃗ r) = − σR ² µ 0

4π

∆θ 3

∆φ 3

X N i=0

X M j=0

a i b j

sin θ ^′ _i ^h (⃗ r − ⃗R ^′ ) × (⃗ω × ⃗R ^′ ) ⁱ · e b _p

|⃗r − ⃗R ij

^′

| ³ , p = x, y, z (16)

(3)

2 Zadania do wykonania

1. Przyjmujemy następujące wartości parametrów: µ 0 = ε 0 = 1, R = 1, σ = 1, N = M = 201,

∆θ = π/N , ∆φ = 2π/M , ⃗ ω = [ω _x , ω _y , 0], ω _x = ω · sin α, ω y = ω · cos α 2. Napisać funkcję obliczającą iloczyn wektorowy (jako wyznacznik) ⃗ c × ⃗d = ⃗g:

v o i d i l o c z y n _ w e k t o r o w y ( d o u b l e c [3] , d o u b l e d [3] , d o u b l e g [ 3 ] ) posłużymy się nią licząc składowe pola magnetycznego ⃗ B.

3. Wyznaczyć V (⃗ r), ⃗ E(⃗ r) oraz ⃗ B(⃗ r) dla ⃗ r = (x, y, 0) w zakresie x ∈ [−L, L], y ∈ [−L, L], L = 3.0 z krokiem ∆x = ∆y = 2L/K, K = 41. Sporządzić mapę rozkładu potencjału V (x, y, 0), wykres wektorowy ⃗ E(x, y, 0) oraz wykres wektorowy ⃗ B(x, y, 0) dla α = 0 i α = π/4. Ponadto proszę zrobić wykres V (0, y, 0) i porównać z rozwiązaniem analitycznym.

3 Uwagi

Do wyznaczania wartości współczynników a _i oraz b _j można użyć funkcji int w s p o l c z y n n i k ( int i , int N ){

int k ;

if ( i ==0 || i == N ) k =1;

el s e if ( i % 2 = = 1 ) k =4;

el s e if ( i % 2 = = 0 ) k =2;

r e t u r n k ; }

wówczas a i = wspolczynnik(i, N ) lub b j = wspolczynnik(j, M ).

Natomiast poniższy pseudokod pozwala wygenerować dane potrzebne do utworzenia wykresów

int w s p o l c z y n n i k ( int i , int N ){

int k ;

if ( i ==0 || i == N ) k =1;

el s e if ( i % 2 = = 1 ) k =4;

el s e if ( i % 2 = = 0 ) k =2;

}

i n i c j a l i z a c j a p a r a m e t r ó w :

R, σ , ∆θ , ∆φ, ∆ _x , ∆ y , L, α, ⃗ ω FOR x F R O M - L TO L S T E P ∆ _x DO

FOR y F R O M - L TO L S T E P ∆ y DO z =0

/ / = = = = z e r o w a n i e z m i e n n y c h p r z e d c a ł k o w a n i e m = = = = = = = = = = = = = = = = = = = = = = = = = = = = v =0

ex =0 ey =0 bx =0 by =0

/ / = = = = = = = = = c a ł k o w a n i e = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

for i f r o m 0 to N s t e p 1 do

(4)

for j f r o m 0 to M s t e p 1 do θ _i ^′ = . . .

φ ^′ _j = . . . x ^′ = . . . y ^′ = . . . z ^′ = . . .

|⃗r − ⃗r ^′ | = . . . W = ^σR _4π

²

^∆θ ₃ ^∆φ ₃

v = v + W a i b j sin θ ^′ _i / |⃗r − ⃗r ^′ |

ex = ex + W a _i b _j sin θ ^′ _i (x − x ^′ _i )/ |⃗r − ⃗r ^′ | ³ ey = ey + . . .

i l o c z y n _ w e k t o r o w y (⃗ ω , ⃗ r ^′ , ⃗ g ⁽¹⁾ ) i l o c z y n _ w e k t o r o w y (⃗ r − ⃗r ^′ , ⃗ g ⁽¹⁾ , ⃗ g ⁽²⁾ ) bx = bx + W a i b j sin θ ^′ _i g x ⁽²⁾ / |⃗r − ⃗r ^′ | ³

by = by + W a i b j sin θ _i ^′ g ⁽²⁾ y / |⃗r − ⃗r ^′ | ³ end do

end do

/ / = = = z a c h o w u j e m y o b l i c z o n e w a r t o ś c i = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = z a p i s do p l i k u : x , y , v , ex , ey , bx , by

END DO END DO

Wykres wektorowy w Gnuplocie tworzymy poleceniem pl o t ’ d a n e . dat ’ u A : B : C : D w v e c t o r s

gdzie: A i B to numery kolumn w których zapisane są x i y, a w kolumnach C i D znajdują się składowe x − owa i y − owa pola wektorowego. W takim przypadku wektory (strzałki są zaczepione w węzłach siatki).

Wektory możemy skalować i zmieniać punkt zaczepnia (np. w środku wektora) scale=5

plot ’dane.dat’ u ($A-$Cscale/2):($B-$Dscale/2):($Cscale):($Dscale) w vectors

4 Przykładowe wyniki

(5)

(a)

V(x,y)

-2 -1 0 1 2

x -2

-1 0 1 2

y

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

(b)

-2 -1 0 1 2

y

x E=[Ex,Ey]

Rysunek 2: (a) Rozkład potencjału V (x, y) oraz (b) pola elektrycznego ⃗ E = [E x , E y ] dla z = 0.

(a)

-2 -1 0 1 2

y

x B=[Bx,By], α=0

(b)

-2 -1 0 1 2

y

x B=[Bx,By], α=π/4

Rysunek 3: Rozkład pola magnetycznego ⃗ B = [B _x , B _y ] dla z = 0 oraz: (a) α = 0 i (b) α = π/4

Projekt 6: modelowanie pola elektrycznego i magnetycznego w otoczeniu obracającej się naładowanej sfery

Projekt 6: modelowanie pola elektrycznego i magnetycznego w otoczeniu obracającej się naładowanej sfery

10 czerwca 2021

1 Wstęp

Na powierzchni sfery o promieniu R znajduje się ładunek o gęstości powierzchniowej σ. Sfera obraca się wokół osi leżącej w płaszczyźnie xy (rys.1). Na zajęciach wyznaczymy rozkład potencjału i pola elektrycznego oraz pola magnetycznego w jej otoczeniu.

Rysunek 1: Środek naładowanej sfery znajduje się w początku układu współrzędnych. Sfera obraca się wokół osi leżącej w płaszczyźnie xy.

1.1 Pole elektryczne

Skorzystamy z zasady superpozycji

V (⃗ r) = 1 4πε 0

Z ρ(⃗ r ′ )

|⃗r − ⃗r ′ | d 3 r ′ (1)

Gęstość ładunku deﬁnujemy w postaci ρ(⃗ r ′ ) = σ · δ(r ′ − R) (σ = const) i wprowadzamy współrzędne sferyczne

V (⃗ r) = 1 4πε 0

Z π 0

dθ ′ Z 2π 0

dϕ ′ Z ∞ 0

dr ′ sinθ ′ r ′2 σδ(r ′ − R)

|⃗r − ⃗r ′ | (2)

= 1

4πε 0 Z π 0

dθ ′ Z 2π 0

dϕ ′ sinθ ′ R 2 σ

|⃗r − ⃗R ′ | (3)

gdzie:

⃗ r →

 



 

x = r sin θ cos φ y = r sin θ sin φ z = r cos θ

R ⃗ ′ →

 



 

R ′ x = R sin θ ′ cos φ ′ R ′ y = R sin θ ′ sin φ ′ R ′ z = R cos θ ′

(4)

Całkę (3) wyznaczymy numerycznie stosując złożenie dwóch kwadratur (2 wymiary) i wzór całkowy parabol (Simpsona). W tym celu na powierzchnii sfery wprowadzamy siatkę węzłów

∆θ = π

N (5)

∆φ = 2π

M (6)

θ → θ i = ∆θ · i, i = 0, 1, . . . , N (7)

φ → φ j = ∆φ · j, j = 0, 1, . . . , M (8)

a całkowanie zastępujemy sumowaniem V (⃗ r) = σR 2

4πε 0

∆θ 3

∆φ 3

X N i=0

X M j=0

a i b j sin θ ′ i

|⃗r − ⃗R

ij | (9)

Uwaga: ⃗ R

i,j oznacza że współrzędne R ′ x , R ′ y , R ′ z zależą od pary indeksów (i, j) ↔ (θ i ′ , φ ′ j ) zgodnie z wzorami (4).

Współczynniki a i /b j są określone następująco

a i =

 



 

a i = 1, i = 0 ∨ i = N

a i = 4, i mod 2 = 1 (i nieparzyste) a i = 2, i mod 2 = 0 (i parzyste)

b j =

 



 

b j = 1, j = 0 ∨ j = M

b j = 4, j mod 2 = 1 (j nieparzyste) b j = 2, j mod 2 = 0 (j parzyste)

(10) Podobnie możemy wyznaczyć numerycznie wektor pola elektrycznego

E(⃗ ⃗ r) = −∇V (⃗r) = 1 4πε 0

Z π 0

dθ ′ Z 2π 0

dϕ ′ R 2 σ sin θ ′ (⃗ r − ⃗R ′ )

|⃗r − ⃗R ′ | 3 (11)

po zastąpieniu całki sumą

E(⃗ ⃗ r) = σR 2 4πε 0

∆θ 3

∆φ 3

X N i=0

X M j=0

a i b j

sin θ i ′ (⃗ r − ⃗R ij

V (⃗ r) = 1 4πε ₀

Z ρ(⃗ r ^′ )

|⃗r − ⃗r ^′ | d ³ r ^′ (1)

Gęstość ładunku deﬁnujemy w postaci ρ(⃗ r ^′ ) = σ · δ(r ^′ − R) (σ = const) i wprowadzamy współrzędne sferyczne

dθ ^′ Z 2π 0

dϕ ^′ Z ∞ 0

dr ^′ sinθ ^′ r ^′2 σδ(r ^′ − R)

|⃗r − ⃗r ^′ | (2)

4πε ₀ Z π 0

dθ ^′ Z 2π 0

dϕ ^′ sinθ ^′ R ² σ

|⃗r − ⃗R ^′ | (3)

R ⃗ ^′ →

R ^′ _x = R sin θ ^′ cos φ ^′ R ^′ _y = R sin θ ^′ sin φ ^′ R ^′ _z = R cos θ ^′

a całkowanie zastępujemy sumowaniem V (⃗ r) = σR ²

a _i b _j sin θ ^′ _i

_ij | (9)

_i,j oznacza że współrzędne R ^′ _x , R ^′ _y , R ^′ _z zależą od pary indeksów (i, j) ↔ (θ _i ^′ , φ ^′ _j ) zgodnie z wzorami (4).

Współczynniki a _i /b _j są określone następująco

a _i = 1, i = 0 ∨ i = N

b _j = 1, j = 0 ∨ j = M

dθ ^′ Z 2π 0

dϕ ^′ R ² σ sin θ ^′ (⃗ r − ⃗R ^′ )

|⃗r − ⃗R ^′ | ³ (11)

E(⃗ ⃗ r) = σR ² 4πε 0

sin θ _i ^′ (⃗ r − ⃗R ij

_ij | ³ (12)

B(⃗ ⃗ r) = − µ ₀ 4π

Z (⃗ r − ⃗r ^′ ) ×⃗j(⃗r ^′ )

|⃗r − ⃗r ^′ | ³ d ³ r ^′ (13)

⃗j(⃗r ^′ ) = ρ(⃗ r ^′ )⃗ v ^′ = σδ(r ^′ − R) ⃗ω × ⃗r ^′ (14) wstawiamy ją do równania (13) a następnie wykonujemy całkowanie po zmiennej radialnej

dθ ^′ Z 2π 0

dϕ ^′ σR ² sin θ ^′ (⃗ r − ⃗R ^′ ) × (⃗ω × ⃗r ^′ )

|⃗r − ⃗R ^′ | ³ (15)

B p (⃗ r) = − σR ² µ 0

sin θ ^′ _i ^h (⃗ r − ⃗R ^′ ) × (⃗ω × ⃗R ^′ ) ⁱ · e b _p

| ³ , p = x, y, z (16)

∆θ = π/N , ∆φ = 2π/M , ⃗ ω = [ω _x , ω _y , 0], ω _x = ω · sin α, ω y = ω · cos α 2. Napisać funkcję obliczającą iloczyn wektorowy (jako wyznacznik) ⃗ c × ⃗d = ⃗g:

Do wyznaczania wartości współczynników a _i oraz b _j można użyć funkcji int w s p o l c z y n n i k ( int i , int N ){

R, σ , ∆θ , ∆φ, ∆ _x , ∆ y , L, α, ⃗ ω FOR x F R O M - L TO L S T E P ∆ _x DO

for j f r o m 0 to M s t e p 1 do θ _i ^′ = . . .

φ ^′ _j = . . . x ^′ = . . . y ^′ = . . . z ^′ = . . .

|⃗r − ⃗r ^′ | = . . . W = ^σR _4π

^∆θ ₃ ^∆φ ₃

v = v + W a i b j sin θ ^′ _i / |⃗r − ⃗r ^′ |

ex = ex + W a _i b _j sin θ ^′ _i (x − x ^′ _i )/ |⃗r − ⃗r ^′ | ³ ey = ey + . . .

i l o c z y n _ w e k t o r o w y (⃗ ω , ⃗ r ^′ , ⃗ g ⁽¹⁾ ) i l o c z y n _ w e k t o r o w y (⃗ r − ⃗r ^′ , ⃗ g ⁽¹⁾ , ⃗ g ⁽²⁾ ) bx = bx + W a i b j sin θ ^′ _i g x ⁽²⁾ / |⃗r − ⃗r ^′ | ³

by = by + W a i b j sin θ _i ^′ g ⁽²⁾ y / |⃗r − ⃗r ^′ | ³ end do

plot ’dane.dat’ u ($A-$Cscale/2):($B-$Dscale/2):($Cscale):($Dscale) w vectors