Từ trường tạo bởi cuộn solenoid

Ta sẽ tính từ trường \(\vec{B}\) tại các điểm trong không gian giới hạn bởi \(\mathrm{xmin}\), \(\mathrm{xmax}\), \(\mathrm{ymin}\), \(\mathrm{ymax}\), \(\mathrm{zmin}\) và \(\mathrm{zmax}\):

1
2
3
4
5
6
7
8
9
10
11

Nx = 20;
Ny = 10;
Nz = 10;
xmin = -1e-2; xmax = 11e-2;
ymin = -1e-2; ymax = 1e-2;
zmin = -1e-2; zmax = 1e-2;

x = linspace(xmin,xmax,Nx);
y = linspace(ymin,ymax,Ny);
z = linspace(zmin,zmax,Nz);
[X,Y,Z] = meshgrid(x,y,z);

Cuộn solenoid có chiều dài \(L=10\,\mathrm{cm}\), bán kính \(R=0.5\,\mathrm{cm}\) mang dòng điện \(I=1\mathrm{A}\) cần có thể mô tả trong Matlab qua hàm số của một đường xoắn ốc. Ở đây ta khai báo một đường xoắn ốc gồm \(\mathrm{N\_turn}=30\) vòng xoắn:

1
2
3
4
5
6
7
8
9
10
11

I = 1; % Ampere
R = 0.5e-2; % Solenoid's radius
L = 10e-2; % Solenoid's length
N_turn = 30; % Number of turns
N_phi = 60;
N1 = N_phi*N_turn;
phi_max = 2*pi*N_turn;
phi = linspace(0,phi_max,N1);
x_wire = linspace(0,L,N1);
y_wire = R*cos(phi);
z_wire = R*sin(phi);

Từ trường \(\vec{B}\) tại mỗi điểm trong không gian là sự chồng chập từ trường tạo ra bởi từng đoạn dây riêng rẽ, theo định luật Biot-Savart:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

Bx = zeros(size(X));
By = zeros(size(Y));
Bz = zeros(size(Z));

kB = mu0/4/pi*I;
cond = [];
for i_wire = 1:N
dx_wire = x_wire(i_wire+1)-x_wire(i_wire);
dy_wire = y_wire(i_wire+1)-y_wire(i_wire);
dz_wire = z_wire(i_wire+1)-z_wire(i_wire);

Delta_x = X-x_wire(i_wire);
Delta_y = Y-y_wire(i_wire);
Delta_z = Z-z_wire(i_wire);
r = sqrt(Delta_x.^2+Delta_y.^2+Delta_z.^2);

Bx = Bx+kB*(dy_wire*Delta_z-dz_wire*Delta_y)./r.^3; % x-component of vector field
By = By+kB*(dz_wire*Delta_x-dx_wire*Delta_z)./r.^3; % y-component of vector field
Bz = Bz+kB*(dx_wire*Delta_y-dy_wire*Delta_x)./r.^3; % z-component of vector field

cond = [cond;find(r<0.1e-2)];
end

Bx(cond) = 0;
By(cond) = 0;
Bz(cond) = 0;

Hình dưới biểu diễn kết quả mô phỏng cho từ trường tạo bởi cuộn solenoid.

Kết quả tính toán từ trường sinh ra bởi cuộn solenoid

Code chương trình Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

function M_field_vec_Solenoid
% Author: Tran Hai Cat
% Lecturer in Physics, HCM University of Technology and Education
% - Dai hoc Su pham Ky thuat Tp. Ho Chi Minh
% Created: 2019.05.02
clc;
clear variables;
close all;

%% CONSTANTS
mu0 = 4*pi*1e-7;

%% PARAMETERS
I = 1; % Ampere
R = 0.5e-2; % Solenoid's radius
L = 10e-2; % Solenoid's length
N_turn = 30; % Number of turns
N_phi = 60;
N1 = N_phi*N_turn;
phi_max = 2*pi*N_turn;
phi = linspace(0,phi_max,N1);
x_wire = linspace(0,L,N1);
y_wire = R*cos(phi);
z_wire = R*sin(phi);

streak_arrow = 1; % 1-MUI TEN CONG, 0-MUI TEN THANG
Nx = 20;
Ny = 10;
Nz = 10;
xmin = -1e-2; xmax = 11e-2;
ymin = -1e-2; ymax = 1e-2;
zmin = -1e-2; zmax = 1e-2;
%% CALCULATION
N = length(x_wire)-1;
x = linspace(xmin,xmax,Nx);
y = linspace(ymin,ymax,Ny);
z = linspace(zmin,zmax,Nz);
[X,Y,Z] = meshgrid(x,y,z);

Bx = zeros(size(X));
By = zeros(size(Y));
Bz = zeros(size(Z));

kB = mu0/4/pi*I;
cond = [];
for i_wire = 1:N
dx_wire = x_wire(i_wire+1)-x_wire(i_wire);
dy_wire = y_wire(i_wire+1)-y_wire(i_wire);
dz_wire = z_wire(i_wire+1)-z_wire(i_wire);

Delta_x = X-x_wire(i_wire);
Delta_y = Y-y_wire(i_wire);
Delta_z = Z-z_wire(i_wire);
r = sqrt(Delta_x.^2+Delta_y.^2+Delta_z.^2);

Bx = Bx+kB*(dy_wire*Delta_z-dz_wire*Delta_y)./r.^3; % x-component of vector field
By = By+kB*(dz_wire*Delta_x-dx_wire*Delta_z)./r.^3; % y-component of vector field
Bz = Bz+kB*(dx_wire*Delta_y-dy_wire*Delta_x)./r.^3; % z-component of vector field

cond = [cond;find(r<0.1e-2)];
end

Bx(cond) = 0;
By(cond) = 0;
Bz(cond) = 0;

%% FIGURES
figure('name','Electric Field Rod','color','k','numbertitle','off');
if streak_arrow
streakarrow3d(X,Y,Z,Bx,By,Bz,2,1);
else
quiver3(X,Y,Z,Bx,By,Bz,'AutoScaleFactor',3);
end
hold on

% Draw the solenoid:
plot3(x_wire,y_wire,z_wire,'linewidth',3);

%% Axis properties
axis equal
axis([xmin xmax ymin ymax zmin zmax]);
rotate3d on

xlabel('X, m');
ylabel('Y, m');
zlabel('Z, m');
view(-20,30);
box off
grid off
set(gca,'color','k','xcolor','w','ycolor','w','zcolor','w')

function hh=streakarrow3d(X,Y,Z,U,V,W,np,AC)
%% Bertrand Dano 03-09
% Copyright 1984-2009 The MathWorks, Inc.
DX=abs(X(1,1,1)-X(1,2,1)); DY=abs(Y(1,1,1)-Y(2,1,1)); DZ=abs(Z(1,1,1)-Z(1,1,2));
DD=min([DX DY DZ]);
ks=DD/100; % Size of the "dot" for the tuft graphs
np=np*10;
alpha = 3; % Size of arrow head relative to the length of the vector
beta = .15; % Width of the base of the arrow head relative to the length
XYZ=stream3(X,Y,Z,U,V,W,X,Y,Z);
Vmag=sqrt(U.^2+V.^2+W.^2);
Vmin=min(Vmag(:));Vmax=max(Vmag(:));
Vmag=Vmag(:);
cmap=colormap;
for k=1:length(XYZ)
F=XYZ(k); [L M]=size(F{1});
if L<np
F0{1}=F{1}(1:L,:);
if L==1
F1{1}=F{1}(L,:);
else
F1{1}=F{1}(L-1:L,:);
end

else
F0{1}=F{1}(1:np,:);
F1{1}=F{1}(np-1:np,:);
end
P=F1{1};
if AC==1
vcol=floor((Vmag(k)-Vmin)./(Vmax-Vmin)*64); if vcol==0; vcol=1; end
COL=[cmap(vcol,1) cmap(vcol,2) cmap(vcol,3)];
else
COL='k';
end
hh=streamline(F0);
set(hh,'color',COL,'linewidth',.5);

if L>1
x1=P(1,1); y1=P(1,2); z1=P(1,3);
x2=P(2,1); y2=P(2,2); z2=P(2,3);
u=x2-x1; v=y2-y1; w=z2-z1; uv=sqrt(u.*u + v.*v);

xa1=x2+u-alpha*(u+beta*(v+eps)); xa2=x2+u-alpha*(u-beta*(v+eps)); xa3=x2+u-alpha*u;
ya1=y2+v-alpha*(v-beta.*(u+eps)); ya2=y2+v-alpha*(v+beta.*(u+eps)); ya3=y2+v-alpha*v;
za1=z2+w-alpha*w; za2=z2+w-alpha*(w+beta.*(uv+eps)); za3=z2+w-alpha*(w-beta.*(uv+eps));

plot3([x2 xa1 xa2 x2 xa3 xa3 x2],[y2 ya1 ya2 y2 ya3 ya3 y2],[z2 za1 za1 z2 za2 za3 z2],'color',COL); hold on
end

end
axis image