Giáo trình Matlab - Chương 8: Phương trình vi phân đạo hàm riêng (Phần 2)

ây nó 
có trị số là 1. Hình dạng hình học của vùng xét cho thấy vec tơ A đối xứng so 
với trục y và đối xứng đổi dấu so với trục x. Do đó ta chỉ cần tính trường trong 
vùng x ≥ 0, y ≥ 0 với điều kiện biên Neumann  ⎟⎠
⎞⎜⎝
⎛ ∇µ A
1nr trên trục x và điều kiện 
biên Dirichlet A = 0 trên trục y. Trường bên ngoài động cơ có thể bỏ qua nên 
trên biên bên ngoài ta dùng điều liện biên Dirichlet A = 0. 
Hình dạng động cơ rất phức tạp, gồm 5 cung tròn và hai hình chữ nhật. 
Ta dùng pdetool GUI, đặt trục x có giới hạn từ [‐1.5 1.5] và trục y có giới hạn [‐
1 1]. Đặt kiểu  ứng dụng  là Magnetostatics và dùng grid  spacing  là 0.1. Mô 
hình  là  tổ hợp của hình  tròn và hình chữ nhật. Ta khởi động pdetool vẽ mô 
hình bằng các lệnh sau(lưu trong file ct8_17.m): 
pdecirc(0,0,1,ʹC1ʹ) 
pdecirc(0,0,0.8,ʹC2ʹ) 
pdecirc(0,0,0.6,ʹC3ʹ) 
pdecirc(0,0,0.5,ʹC4ʹ) 
pdecirc(0,0,0.4,ʹC5ʹ) 
pderect([–0.2 0.2 0.2 0.9],ʹR1ʹ) 
pderect([–0.1 0.1 0.2 0.9],ʹR2ʹ) 
pderect([0 1 0 1],ʹSQ1ʹ) 
Nhập lệnh sau (C1 + C2 + C3 + C4 + C5 + R1 + R2)*SQ1 rồi bấm icon ∂Ω để đưa 
mô hình về góc phần tư thứ nhất. Ta cần bỏ một số biên của vùng con. Dùng 
shift‐click,  chọn  biên  và  bỏ  nó  bằng  cách  dùng Remove  Subdomain  Border 
trong menu Boundary cho đến khi mô hình có 4 vùng: stator, rotor, khe hở và 
cuộn dây. Trước khi chuyển sang PDE mode, chọn biên dọc theo trục x và đặt 
điều kiện biên Neumann với g = 0 và q = 0. 
  Trong  PDE mode  chọn  Show  Subdomain  Labels. Nhấp  đúp  lên  từng 
vùng để xác định các thông số của PDE. 
  • trong cuộn dây cả µ và J đều bằng 1, như vậy không cần thay đổi giá trị  
177
mặc định. 
  • trong stator và rotor µ là phi tuyến và xác định bằng phương trình trên. 
Nhập trị của µ: 5000./(1+0.05*(ux.^2+uy.^2))+200 và J = 0. 
  • trong không khí µ = 1 và J = 0 
Ta khởi gán lưới và tiếp tục bằng cách mở hộp thoại Solve Parameters và 
chọn cách giải phi tuyến. Giải bài toán và vẽ B bằng mũi tên và đường đẳng 
thế bằng contour. Bài toán lưu trong ct8_18.m. 
  e. Trường  điện  từ  của nguồn ac: Bài  toán về  trường  điện  từ của nguồn 
xoay chiều xuất hiện khi ta nghiên cứu trường của các động cơ, m.b.a, vật dẫn 
có dòng điện xoay chiều. Chúng ta sẽ khảo sát trường trong môi trường điện 
môi  đồng  nhất  có  các  hệ  số  ε  &  µ  và  không  có  điện  tích  trong  toàn 
miền.Trường thoả mãn hệ phương trình Maxwell: 
JEH
HE
rrr
rr
+∂
∂ε=×∇
∂
∂µ−=×∇
t
t  
Khi không có dòng điện, hệ trường thoả mãn phương trình sóng với vận 
tốc truyền sóng là  εµ : 
0
t
0
t
2
2
2
2
=∂
∂µε−∆
=∂
∂µε−∆
EH
HE
rr
rr
Chúng ta nghiên cứu một điện môi đồng nhất không có điện tích, có các 
hệ số điện môi là ε, độ từ thẩm là µ và độ dẫn điện là σ. Mật độ dòng điện là: 
 và các sóng bị tắt dần do điện trở. Phương trình đối với E là: EJ
rr σ=
0
tt 2
2
=∂
∂µε−∂
∂µσ−∆ EEE
rrr
Phương trình đối với H cũng có dạng tương tự. 
Trường hợp trường biến thiên điều hoà ta dùng dạng phức, thay E bằng  tjce
ωE . 
Trường hợp bài toán phẳng ta có Ec = (0,0,Ec) và J = (0,0,Jejωt) và từ trường là: 
  cyx j
1)0,H,H( EH
rr ×∇ωµ−==   
Phương trình vô hướng đối với Ec trở thành: 
  0E)j(E1 c
2
c =εω−ωσ+⎟⎠
⎞⎜⎝
⎛ ∇µ∇−  
178
  Phương trình này được dùng trong PDE Toolbox với dạng ứng dụng AC 
Power  Electromagnetics. Nó  là  phương  trình Helmholz  phức mô  tả  sự  lan 
truyền  của một  sóng  điện  từ phẳng  trong môi  trường  điện môi không hoàn 
hảo  và  dẫn  điện  tốt(σ  >> ωε). Hệ  số  điện môi  phức  εc  định  định  nghĩa  là: 
ωσ−ε=ε jc . Điều kiện biên ở bề mặt phân cách hai môi trường là điều kiện tự 
nhiên đối với công thức biến phân nên ta không cần chú ý. Các thông số PDE 
cần nhập vào hộp  thoại PDE Specification  là  tần số góc ω, độ  từ  thẩm µ, độ 
dẫn điện σ và hệ số điện môi ε. Điều kiện biên kết hợp với mode này là điều 
kiện biên Dirichlet mô  tả giá  trị của Ec  trên biên và điều kiện biên Neumann 
mô tả đạo hàm của Ec theo hướng pháp tuyến. Điều này tương đương với việc 
cho thành phần tiếp tuyến của vec tơ cường độ từ trường H: 
  ⎟⎠
⎞⎜⎝
⎛ ∇µω= ct E
1njH  
Các đại lượng có thể tính từ nghiệm là: 
  EB
rr ×∇ω=
j  
và  
J = σE. 
  Các đại lượng E, J, B, H và nhiệt lượng Q = Ec*Ec/σ có thể vẽ ra. Vec tơ 
B,H vẽ được nhờ dùng mũi tên. 
Ta xác định hiệu ứng mặt ngoài khi dòng điện xoay chiều chạy trong vật 
dẫn bằng đồng có tiết diện tròn. Độ dẫn điện của đồng là σ = 57.106 và độ từ 
thẩm là 1, nghĩa là µ = 4π.10‐7. Tại tần số f = 50Hz, ω2ε ≈ 0 và có thể bỏ qua. Do 
hiện tượng cảm ứng nên mật độ dòng điện bên trong thanh dẫn nhỏ hơn mặt 
ngoài là nơi có JS = 1 và điều kiện biên của trường là Ec = 1/σ. Trong trường hợp 
này có thể tìm nghiệm giải tích dạng: 
)kR(J
)kr(JJJ
0
0
S=  
Trong đó  ωµε= jk , và R  là bán kính của dây, r  là khoảng cách đến  tâm và 
J0(x) là hàm Bessel loại 1 bậc zero. 
Khởi động pdetool GUI và dùng mode AC Power Electromagnetics. Vẽ 
hình  tròn  bán  kính  0.5  để  biểu  diễn  mặt  cắt  ngang  của  dây  dẫn  (file 
ct8_19.m)và xác định điều kiện biên. Dùng tuỳ chọn Select All để chọn tất cả 
biên và nhập 1/57E6 và cho r trong hộp thoại Boundary Condition để xác định 
điều kiện biên Dirichlet(E = J/σ). Mở hộp thoại PDE Specification và nhập các 
thông số. Tần số góc ω = 2*pi*50. Đánh dấu ô Adaptive mode. Đặt số tam giác 
179
là Inf và maximum numbers of refinements là 1. Khởi gán lưới và giải phương 
trình. Vẽ mật độ dòng điện là 3‐D. Bài toán lưu trong ct8_19.m. 
  f. Môi trường dẫn dc: Khi với tính toán quá trình điện phân và điện trở 
nối đất ta gặp một môi trường dẫn có độ dẫn σ và một dòng điện ổn định. Mật 
độ dòng điện J liên quan với cường độ điện trường E bằng phương trình J=σE. 
Kết hợp  tính  liên  tục  của phương  trình ∇J = Q  (Q  là nguồn dòng) với  định 
nghĩa điện thế U cho ta phương trình Poisson: 
   0)U( =∇σ∇−
Ta chỉ có 2 thông số PDE  là độ dẫn σ và nguồn dòng điện Q. Điều kiện biên 
Dirichlet gán các giá trị của điện thế U vào biên, thường là vật dẫn bằng kim 
loại. Điều kiện biên Neumann đòi hỏi giá trị pháp tuyến của mật độ dòng điện 
(nσ(∇(U))) đã cho. Cũng có thể mô tả điều kiện biên Neumann tổng quát hoá 
được xác định bằng nσ(∇(U)) + qU = g, trong đó q có thể coi là một lớp mỏng 
dẫn điện. Điện thế U, cường độ điện trường E và mật độ dòng điện J có thể vẽ 
ra. Ta muốn thấy các đường dòng điện(trường vec tơ  của J) và các đường thế 
của U. Các đường thế trực giao với đường dòng điện khi σ đẳng hướng. 
Ta xét hai thanh dẫn kim loại hình tròn được đặt trong một mặt phẳng, 
dẫn điện mỏng như tờ giấy thấm thấm đẫm nước biển. Mô hình vật lí của bài 
toán gồm phương trình Laplace đối với điện thế U: 
0)U( =∇σ∇−  
và các điều kiện biên: 
  • U = 1 trên thanh dẫn tròn bên trái 
  • U = ‐1 trên thanh dẫn tròn bên phải 
  • điều kiện biên Neumann tự nhiên trên biên ngoài ∂U/∂n  = 0 
Độ dẫn σ = 1. 
Trong pdetool GUI dùng Conductive Media mode. Vẽ  tờ giấy  thấm  là 
hình chữ nhật R1 có các góc (‐1.2, ‐0.6),(1.2, ‐0.6),(1.2, 0.6) và (‐1.2, 0.6). Vẽ thêm 
2 hình tròn C1 và C2 có bán kính 0.2 và tâm tại (‐0.6, 0) và (0.6, 0). Bài toán 2‐D 
(file ct8_20.m)được biểu diễn trong miền R1 ‐ (C1 + C2). 
  Chọn toàn bộ biên ngoài và đặt điều kiện biên Neumann vào hộp thoại 
Boundary Condition. Đối với  thanh dẫn  tròn bên  trái  ta nhập điều kiện biên 
Dirichlet U = 1. Thanh dẫn tròn bên phải có điều kiện biên Dirichlet u = ‐1. 
  Tiếp  theo mở hộp  thoại PDE Specification và nhập q  =  0. Giá  trị mặc 
định của σ =1 không cần thay đổi. 
  Khởi gán  lưới và  tinh chỉnh 2  lần. Giải phương  trình bằng nhấn nút =. 
Xem  J bằng cách vẽ giá  trị  tuyệt đối và dùng contour plot và  trường vec  tơ 
bằng cách dùng mũi tên. 
180
  g. Truyền nhiệt: Phương trình truyền nhiệt có dạng: 
  )TT(hQ)Tk(
t
TC ext −+=∇∇−∂
∂ρ  
Đây là  phương trình parabolic có các thông số: 
  • mật độ ρ 
  • nhiệt dung C 
  • hệ số dẫn nhiệt k 
  • nguồn nhiệt Q 
  • hệ số truyền nhiệt bằng đối lưu h 
  • nhiệt độ bên ngoài Text
Các  điều kiện biên  có  thể  là  điều kiện biên Dirichlet(cho nhiệt  độ  trên biên) 
hay  điều  kiện  biên  Neumann(cho  dòng  nhiệt  n.(k∇(T)).  Điều  kiện  biên 
Neumann tổng quát n.(k∇(T) + qt = g (q là hệ số truyền  nhiệt)cũng có thể được 
dùng. Ta có thể xem nhiệt độ, gradient nhiệt độ và dòng nhiệt k∇(T). 
Ta xét bài toán truyền nhiệt với các vật liệu khác nhau. Bài toán 2‐D gồm 
một hình vuông quay 450(hình thoi). Vùng hình vuông khác bằng vật  liệu có 
hệ số dẫn nhiệt 10 và mật độ 2. Vùng hình thoi có nguồn nhiệt phân bố đều 
với trị số 4, có hệ số dẫn nhiệt 2 và mật độ 1. Cả 2 vùng có nhiệt dung 0.1. 
Khởi động pdetool ở kiểu Heat Transfert. Đặt giới hạn của trục x và y là 
[‐0.5 3.5] và chọn Axis Equal từ menu Option(file ct8_21.m). 
  Vẽ hình vuông có các góc (0, 0), (3, 0), (0, 3) và (3, 3) và hình thoi có các 
góc (1.5 0.5), (2.5 1.5), (1.5 2.5) và (0.5 1.5). 
  Nhiệt độ bằng 0 trên tất các các biên ngoài và như vậy ta không cần thay 
đổi điều kiện biên mặc định. 
  Nhấp  đúp  vào  từng  vùng  để  nhập  các  thông  số  PDE.  Ta muốn  giải 
phương trình truyền nhiệt parabolic vì vậy cần chọn cách giải parabolic. 
  Trong vùng hình vuông nhập mật độ 1, nhiệt dung 0.1, hệ số dẫn nhiệt 
10. Do không có nguồn nhiệt nên nhập 0. 
  Trong vùng hình thoi nhập mật độ 1, nhiệt dung 0.1, hệ số dẫn nhiệt 2, 
nguồn nhiệt 4. Số hạng (Text‐T) không dùng nên nhập h = 0. 
  Do  ta giải phương  trình PDE  động nên  cần  cho  điều kiện  đầu và  thời 
gian  tính nghiệm. Vì vậy mở hộp  thoại Solve Parameter. Bài  toán xảy ra rất 
nhanh, đạt giá  trị xác  lập  trong 0.1 đơn vị  thời gian. Để bắt được phần quan 
trọng của đặc tính động, nhập time [0:0.01:0.1] như là vec tơ thời gian để giải 
phương trình nhiệt. Đặt giá trị đầu của nhiệt độ là 0. Giải bài toán. Mặc định, 
nhiệt độ cuối quá trình được vẽ. Cách tốt nhất để nhìn đặc tính động là hoạt 
181
hình hoá nghiệm. Khi hoạt hình hoá,  đánh dấu Height và  chọn Plot  in  x‐y 
grid. 
182