Giáo trình Matlab - Chương 8: Phương trình vi phân đạo hàm riêng (Phần 2)
menu con Select All trong menu Edit và chọn điều kiện biên Neumann cho tất
cả các biên. Sau đó sửa lại điều kiện biên cho hai phía của thanh. Phía trái chọn
điều kiện biên Dirichlet với r = 100. Phía bên phải chọn điều kiện biên
Neumann với g = ‐10. Bước tiếp theo là mở hộp thoại PDE Specification và
nhập vào các hệ số của PDE. PDE parabolic tổng quát mà PDE Toolbox xử lí
có dạng:
∂u
∂t
d
−∇.(∇u) + au = f
với điều kiện đầu u0 = u(t0) và thời gian tính nghiệm mô tả trong mảng tlist.
Như vậy trong trường hợp này ta có d = 1, c = 1, a = 0 và f = 0. Khởi gán các
lưới và làm tinh lại. Điều kiện đầu u0 = 0 và khoảng thời gian được nhập vào
là [0:0.5:5]. Ta nhập chúng vào hộp thoại Solve Parameters từ menu Solve. Bây
giờ ta có thể giải bài toán. Để thấy được quá trình truyền nhiệt ta đánh dấu
vào ô Animation trong hộp thoại Plot selection. Nên chọn màu là colormap
hot. Chú ý là nhiệt độ của khối tăng rất nhanh. Bài toán này được lưu trong
ct8_9.m.
b. Phân bố nhiệt trong thanh phóng xạ: Bài toán phân bố nhiệt này là
một ví dụ về bài toán 3‐D PDE parabolic được biến đổi thành bài toán 2‐D nhờ
dùng toạ độ trụ. Ta khảo sát một thanh phóng xạ hình trụ. Tại cuối bên trái
của thanh nhiệt được gia tăng liên tục. Đầu cuối bên phải có nhiệt độ không
đổi. Tại biên bên ngoài, nhiệt được trao đổi với mô trường bằng truyền nhiệt.
Tại một thời điểm,nhiệt độ được tạo ra không đồng đều trong toàn bộ thanh
do quá trình phóng xạ. Giả sử ban đầu nhiệt độ bằng 0. Điều này đưa tới bài
toán sau:
∂u
∂t
ρC
−∇.(k∇u) = f
Trong đó ρ là mật độ, C là nhiệt dung riêng của thanh, k là hệ số dẫn nhiệt và f
là nguồn nhiệt phóng xạ. Mật độ của kim loại là 7800kg/m3, nhiệt dung riêng
là 500Ws/kg0C, độ dẫn nhiệt là 40W/m0C. Nguồn nhiệt là 20000W/m3. Nhiệt độ
ở một đầu thanh là 1000C. Nhiệt độ môi trường bên ngoài là 1000C và hệ số
truyền nhiệt là 50W/m20C. Dòng nhiệt ở cuối bên trái là 5000 W/m2. Nhưng đây
là bài toán hình trụ, như vậy ta cần biến đổi phương trình, dùng các toạ độ trụ
r , z và θ. Do tính đối xứng, nghiệm không phụ thuộc θ. Như vậy phương trình
đã biến đổi là:
∂u ∂
∂t ∂r
∂u
∂r
∂
∂z
∂u
∂z
⎛
⎜
⎝
⎞
⎟
⎠
⎛
⎜
⎝
⎞
⎟
⎠
rρC
−
kr
−
kr
= fr
Điều kiện biên là:
168
r
•
n.(k∇u) = 5000 ở đầu cuối bên trái của thanh(điều kiện biên Neumann).
r
Do điều kiện Neumann tổng quát hoá trong PDE Toolbox là n.(c∇u)+qu = g và
c phụ thuộc vào r trong bài toán này( c= kr),điều kiện biên này được biểu diễn
r
bằng biểu thức n.(c∇u) = 5000r.
•
u = 100 tại đầu cuối bên phải của thanh(điều kiện biên Dirichlet)
n.(k∇u)= 50(100‐u) tại biên bên ngoài(điều kiện biên Neumann tổng
r
•
quát hoá). Trong PDE Toolbox nó được biểu diễn bằng:
r
n.(c∇u)+ 50r.u = 50r.100.
•
trục của hình trụ r = 0 không phải là biên trong bài toán gốc nhưng khi
r
biến đổi thành 2‐D thì lại là biên.Ta phải cho một điều kiện biên n.(c∇u) = 0 tại
đây. Giá trị đầu là u(t0) = 0
Mô hình thanh là hình chữ nhật dọc theo trục x và trục y hướng r. Ta vẽ
hình chữ nhật với các góc (‐1.5,0), (1.5,0), (1.5,0.2), (‐1.5,0.2), nghĩa là cần nhập
các số [‐1.5 0.0 3 0.2] vào Object Dialog của phần tử R1. Nhập điều kiện biên
Neumann cho đầu cuối bên trái với q = 0 và g = 5000*y. Nhập điều kiện biên
Dirichlet cho đầu cuối bên phải với h = 1và r = 100. Đối với biên ngoài dùng
điều kiện biên Neumann với q = 50*y và g = 50*y*100. Trên trục ta dùng điều
kiện biên Neumann với q = 0 và g = 0. Các hệ số của phương trình c = 40*y, a =
0, d = 7800*500*y và f =20000*y.
3. Các ví dụ về bài toán hyperbolic:
a. Phương trình sóng: Ta khảo sát sóng tạo ra từ dao động của một màng
hình vuông có các góc (‐1,‐1),(‐1,1),(1,‐1) và (1,1). Phương trình dao động có
dạng:
∂2u
− ∆u = 0
∂t2
∂u
∂n
⎛
⎜
⎝
⎞
⎠
Màng được cố định(u = 0) tại cạnh phải và cạnh trái và tự do
= 0 ở
⎟
∂u(t )
⎛
⎜
⎝
⎞
⎟
⎠
0
cạnh trên và cạnh dưới. Ngoài ra, ta cần giá trị đầu u(t0) và
. Giá trị
∂t
đầu phải khớp với điều kiện biên. Nếu ta bắt đầu tại t =
π
2
⎛
⎜
⎝
⎞
⎟
⎠
sin
y
π
2
∂u(0)
∂t
⎛
⎝
⎞
⎠
0,thìu(0) = arctan cos x và
= 3sin(πx)e
là các giá trị đầu thoả mãn
⎜
⎟
điều kiện biên.
Ta dùng PDE Toolbox với mode Generic Scalar. Vẽ hình chữ nhật với
các góc như trên, nghĩa là ta phải điền vào Object Dialog các số: [ ‐1 ‐1 2 2].
Sau đó ta xác định điều kiện biên và khởi gán lưới. Mở hộp thoại PDE
169
Specification để nhập các giá trị của hệ số của phương trình. Do phương trình
tổng quát có dạng:
∂2u
d
− ∇.(c∇u) + au = f
∂t2
nên với phương trình sóng ta có c = 1, a = 0, f = 0 và d = 1.
Trước khi giải phương trình chọn Parameters... từ menu Solve để mở
hộp thoại Solve Parametes. Trong mục times nhập linspace(0, 5, 31), giá trị
đầu của u ta nhập atan(cos(pi/2*x)), và đạo hàm của u bằng:
3*sin(pi*x).*exp(sin(pi/2*y)). Cuối cùng nhấn nút = để giải phương trình. Bài
toán này được lưu trong ct8_11.m.
4. Các ví dụ về các bài toán giá trị riêng:
a. Các giá trị riêng và hàm riêng của màng dạng L: Bài toán tìm các giá trị
riêng và các hàm riêng tương ứng của một màng dạng L hấp dẫn người dùng
MATLAB, vì vẽ hàm riêng đầu tiên là logo của MathWorks. Thực tế, có thể so
sánh các giá trị riêng và hàm riêng tính bằng PDE Toolbox và các giá trị riêng
và hàm riêng tạo bởi hàm membrane của MATLAB.
Bài toán được giải với tất cả các mode riêng và có các giá trị riêng nhỏ
hơn 100 đối với bài toán PDE eigenmode:
‐∆u = λu
trên hình dạng của màng L, u = 0 trên biên(điều kiện biên Dirichlet).
Kích hoạt pdetool và kiểm tra xem ta đã ở Generic Scalar chưa. Sau đó
vẽ hình L có các góc (0,0), (‐1,0), (‐1,‐1), (1,‐1),(1,1), và (0,1) bằng cách dùng nút
polygon. Không cần xác định điều kiện biên đối với bài toán này vì điều kiện
biên mặc định u = 0 trên biên là phù hợp. Do đó ta có thể thực hiện bước tiếp
theo là tạo lưới. Sau đó ta tinh chỉnh lại lưới hai lần. Xác định bài toán giá trị
riêng PDE cũng dễ. Ta mở hộp thoại PDE Specification và chọn Eigenmodes.
Các giá trị mặc định của các hệ số của phương trình c =1, d = 1 và a = 0 khớp
với bài toán. Như vậy ta có thể thoát khỏi hộp thoại PDE Specification bằng
cách nhấn nút OK. Ta mở hộp thoại Solve Parameters và xác nhận phạm vi
[0,100] là đúng.Cuối cùng giải bài toán bằng nhấn nút =. Giá trị riêng đầu
tiên(nhỏ nhất) được hiển thị. Ta tìm số lượng giá trị riêng trên đường thông tin
ở cuối GUI. Ta có thể mở hộp thoại Plot Selection và chọn giá trị riêng cần vẽ
bằng cách chọn menu có giá trị riêng tương ứng. Bài toán này được lưu trong
ct8_12.m.
b. Các giá trị riêng và mode riêng của màng hình vuông: Ta nghiên cứu
các giá trị riêng và mode riêng của một hình vuông. Các góc của hình vuông là
170
(‐1,‐1), (‐1,1), (1,1), and (1,‐1). Điều kiện biên như sau:
•
trên biên bên trái, điều kiện biên Dirichlet u = 0
∂u
∂n
•
trên biên trên và dưới, điều kiện biên Neumann
= 0
∂u 3
− u = 0.
∂n 4
•
trên biên bên phải, điều kiện biên Neumann tổng quát hoá
Bài toán PDE giá trị riêng là:
‐∆u = λu
Chúng ta quan tâm đến các giá trị riêng nhỏ hơn 10 và mode riêng tương ứng
sao cho phạm vi nghiên cứu là [ ‐Inf 10]. Dấu của điều kiện biên Neumann
tổng quát hoá sao cho ta có giá trị riêng âm.
Ta dùng pdetool GUI theo kiểu Generic Scalar. Vẽ hình vuông dùng
menu Draw. Sau đó ta xác định điều kiện biên trên biên bên phải là g = 0 và q =
‐3/4(điều kiện biên Neumann); điều kiện biên bên trái là u = 0(điều kiện biên
Dirichlet); điều kiện biên trên và dưới là g = 0 và q = 0(điều kiện biên
Neumann). Khởi gán lưới và tinh chỉnh lại một lần. Các hệ số của phương
trình là: c = 1, a = 0 và d = 1. Khoảng giải bài toán trong Solve Parameters là [‐
Inf 10]. Cuối cùng nhấn = để giải bài toán. Bài toán được lưu trong ct8_13.m.
5. Các dạng ứng dụng:
a. Các kiểu ứng dụng và GUI: PDE Toolbox có thể áp dụng cho một số rất
lớn bài toán khoa học và kĩ thuật. Các kiểu ứng dụng có thể có của nó là:
• Generic scalar (kiểu mặc định)
• Generic system
• Structural Mechanics ‐ Plane Stress
• Structural Mechanics ‐ Plane Strain
• Electrostatics
• Magnetostatics
• AC Power Electromagnetics
• Conductive Media DC
• Heat Transfer
• Diffusion
Chú ý là nếu dùng GUI, bài toán bị giới hạn ở hàm 2 biến. Sử dụng các hàm
dòng lệnh sẽ không gặp giới hạn này. Lựa chọn mode được thực hiện trên
menu và khi đó các hệ số và điều kiện biên mặc định sẽ được tự động thay đổi
cho phù hợp.
b. Cơ học kết cấu ‐ ứng suất bề mặt: Trong cơ học kết cấu, các phương
171
trình quan hệ giữa ứng suất và lực kéo phát sinh từ sự cân bằng lực trong vật
thể. Ứng suất bề mặt là một điều kiện phổ biến trong một đĩa phẳng trong mặt
phẳng xy chịu tải chỉ trong mặt phẳng không liên quan đến hướng z. Quan hệ
lực ‐ ứng suất có thể được viết trong điều kiện đẳng hướng và đẳng nhiệt:
⎛
⎞
⎛
⎞
⎛
⎞
σx
1 ν
ν 1
0 0
0
εx
⎜
⎟
⎜
⎟
⎜
⎟
E
⎜
⎜
⎟
⎟
⎟
σ
=
0
ε
⎜
⎜
⎝
⎟
⎟
⎠
⎜
⎜
⎝
⎟
⎟
⎠
y
y
1− ν2
1− ν
τxy
γxy
⎜
⎝
2 ⎠
Trong đó σx và σy là ứng suất theo hướng x và y; τxy là ứng suất cắt. Tính chất
của vật liệu được biểu diễn như là kết hợp của E‐mô đun đàn hồi hay mô đun
Young và ν‐tỉ số Poisson. Sự biến dạng của vật liệu được mô tả bằng sự dịch
chuyển theo hướng x và hướng y là u và v và được xác định như sau:
∂u
∂x
∂v
∂y
∂u ∂v
∂y ∂x
εx =
;εy =
;γxy=
+
Phương trình cân bằng lực là:
∂τxy
∂σx
∂x
−
−
= Kx
∂y
∂σy ∂τxy
−
−
= Ky
∂y
∂x
Trong đó Kx và Ky là lực khối.
Kết hợp các quan hệ trên, chúng ta có phương trình dịch chuyển :
− ∇.(c ⊗ ∇u) = k
Trong đó c là hạng của tenxơ bậc 4 mà ta có thể viết như 4 ma trận 2*2:
2G + µ 0
⎛
⎞
c11 = ⎜
⎟
⎟
⎜
0
G
⎝
⎛
⎠
0 µ
G 0
0 G
µ 0
G
⎞
c12 = ⎜
⎟
⎟
⎜
⎝
⎛
⎠
⎞
c21 = ⎜
⎟
⎟
⎜
⎝
⎛
⎠
0
⎞
c22 = ⎜
⎟
⎟
⎜
0 2G + µ
⎝
⎠
Trong đó G là mo đun cắt, xác định bằng:
E
G =
2(1+ ν)
và µ được xác định bằng:
172
ν
1− ν
µ = 2G
và
K
Ky
⎛
⎞
x
⎜
⎜
⎟
⎟
k =
⎝
⎠
là lực khối.
Đây là phương trình PDE elliptic của một hệ(u là 2 chiều) nhưng ta chỉ
cần để lựa chọn kiểu ứng dụng Structural Mechanics, Plan Stress và nhập các
thông số phụ thuộc vật liệu E và ν và lực khối vào hộp thoại PDE Specification.
Trong kiểu này chúng ta cũng giải bài toán giá trị riêng được mô tả bằng:
−∇.(c ⊗∇u) = λdu
ρ 0
0 ρ
⎛
⎞
d = ⎜
⎟
⎟
⎜
⎝
⎠
với ρ là mật độ và cũng được nhập bằng hộp thoại PDE Specification. Trong
hộp thoại Plot Selection, dịch chuyển theo trục x và trục y là u và v và giá trị
tuyệt đối của vec tơ dịch chuyển (u,v) có thể nhìn thấy bằng cách dùng màu,
đường đồng mức hay độ cao z và trường vec tơ dịch chuyển (u,v) có thể vẽ
bằng cách dùng mũi tên hay lưới biến đổi. Hơn nữa,ta có thể chọn từ 15 biểu
thức ten xơ vô hướng:
∂u
∂x
∂u
∂y
∂v
•
•
ux =
uy =
•
•
vx =
vy =
∂x
∂v
∂y
•
•
•
•
•
•
•
•
•
•
exx, lực căng theo hướng x(εx)
eyy, lực căng theo hướng y (εy)
exy, lực cắt(γxy)
sxx, ứng suất theo hướng x(σx)
syy, ứng suất theo hướng y(σy)
sxy, ứng suất cắt (τxy)
e1, lực chính thứ nhất (ε1)
e2, lực chính thứ 2 (ε2)
s1, ứng suất chính thứ nhất(σ1)
s2, ứng suất cơ bản thứ 2(σ2)
173
•
ứng suất hiệu dụng theo Mises ( σ12 + σ22 − σ1σ2 )
Ta khảo sát một tấm thép bị nén bởi một lực vuông góc với cạnh bên
dưới và kéo ở đoạn cắt tròn bên trên. Các cạnh khác tự do. Tấm thép có các đặc
tính sau: kích thước 1m*1m; dày 1mm; ghép 1/3‐1/3m; vết cắt tròn từ (2/3,1)
đến (1,2/3). Mo đun Young là 196.103MN/m2; tỉ số Poisson 0.31. Biên cong chịu
một lực kéo ra ngoài với trị số 500N/m. Ta cần mô tả lực kéo trên bề mặt nên ta
chia cho chiều dày 1mm và như vậy có lực 0.MN/m2. Ta cần tính lực và ứng
suất theo hướng x và hướng y.
Với bài toán này ta chọn mode Structural Mechanics, Plane Stress. Mô
hình CSG có thể thực hiện bằng cách vẽ đa giác có các góc với x = [0 2/3 1 1 1/3
1/3 0] và y = [1 1 2/3 0 0 1/3 1/3] và hình tròn có tâm tại x = 2/3, y = 2/3 và bán
kính1/3. Mô hình là P1+E1. Tiếp đó chọn Boundary Mode để mô tả điều kiện
biên. Trước hết bỏ hết các biên của vùng con bằng cách chọn Remove All
Subdomain Borders từ menu Boundary. Hai biên ở góc dưới phải bị nén, nghĩa
là điều kiện biên Dirichlet với 0 bị bỏ. Góc cắt tròn có điều kiện biên Neumann
với q = 0 và g1 = 0.5*nx và g2 = 0.5*ny. Các biên còn lại có điều kiện biên
Neumann với q = 0 và g = 0. Bước tiếp theo là mở hộp thoại PDE Selection và
nhập các thông số của PDE. Do không có lực nên Kx và Ky bằng 0 và ρ không
dùng trong mode này. Vật liệu đồng nhất nên E và ν áp dụng cho toàn miền.
Khởi gán lưới và tinh chỉnh lại nó. Giải bài toán bằng cách nhấn nút = .
Bài ttoán được lưu trong ct8_14.m.
c. Cơ học kết cấu‐lực căng bề mặt: Trạng thái biến dạng trong đó không
có dịch chuyển theo hướng z và dịch chuyển theo hướng x và y là hàm của x
và y mà không phải của z được gọi là lực căng bề mặt. Ta có thể giải bài toán
về lực căng bề mặt với PDE Toolbox bằng chọn kiểu ứng dụng Structural
Mechanics, Plane Strain. Quan hệ lực ‐ ứng suất chỉ khác ít so với trường hợp
ứng suất phẳng và các chọn lựa như trên được dùng. Giao diện ứng dụng là
đồng nhất với 2 kiểu ứng dụng cơ học kết cấu. Sự biến dạng trong bài toán lực
căng bề mặt khác với phương trình ứng suất phẳng ở chỗ:
•
thông số m trong ten xơ c được xác định là:
ν
µ = 2G
1− 2ν
•
ứng suất hiệu dụng Mises được tính bằng:
(
σ12 + σ22 )(ν2 − ν + 1
)
− σ1σ2
2ν2 − 2ν − 1
)
Bài toán tương tự như trên được lưu trong ct8_15.m.
174
c. Điện trường tĩnh: Các ứng dụng liên quan đến điện trường tĩnh bao
gồm các thiết bị điện cao áp, các dụng cụ điện tử và các tụ điện. Chữ “tĩnh”
bao hàm ý là sự thay đổi của các đại lượng theo t rất chậm và bước sóng rất
lớn so với kích thước của vùng đang xét. Trong các bài toán tĩnh điện, điện thế
vô hướng V quan hệ với cường độ điện trường E qua biểu thức:
r
E = −∇U
Dùng phương trình Maxwell :
r
∇D = ρ
r
r
và quan hệ D = εE ta có phương trình Poisson:
− ∇.(ε∇U) = ρ
Trong đó ε là hệ số điện mối và ρ là mật độ điện tích không gian.
Sử dụng các ứng dụng PDE Toolbox kiểu Electrostatics, ta có thể giải các
bài toán tĩnh điện được mô hình hoá bằng phương trình trên. Hộp thoại PDE
Specification chứa mục vào cho ε và ρ. Điều kiện biên đối với các bài toán tĩnh
điện có thể là điều kiện biên Dirichlet hay Neumann. Với điều kiện biên
Dirichlet, điện thế tĩnh U được mô tả trên biên. Với điều kiện biên Neumann,
r
điện tích bề mặt n(ε∇U) được mô tả trên biên. Để xem các nghiệm của bài toán
điện tĩnh phần vẽ bao gồm điện thế U, cường độ điện trường E, và vec tơ dịch
chuyển điện D.
Ta khảo sát bài toán xác định điện thế tĩnh trong một khung hình vuông
chứa đầy không khí với biên trong có cạnh dài 0.2 và biên ngoài có cạnh 0.5.
Tại biên trong điện thế là 1000V. Tại biên ngoài điện thế là 0V. Không có điện
tích trong vùng khảo sát. Như vậy ta có phương trình Laplace:
∆U = 0
với điều kiện biên Dirichlet U = 1000 ở biên trong và U = 0 tại biên ngoài.
Sau khi chọn kiểu ứng dụng là Electrostatics ta vẽ vùng xác định của bài
toán là 2 hình vuông có cạnh 0.5 gọi là R1 và 0.2 gọi là R2(file ct8_16.m). Vùng
2‐D ta xét là R2 ‐ R1. Tiếp đó ta xác định điều kiện biên. Chọn tất cả biên trong
bằng cách dùng Shift‐click và đặt điều kiện biên Dirichlet 1000. Điều kiện biên
bên ngoài là 0. Mở hộp thoại PDE Specification và đặt rho = 0. Hệ số điện môi
có thể đặt bằng 1 vì nó là hằng nên không ảnh hưởng đến kết quả. Khởi gán
lưới rồi nhấn nút = để giải phương trình. Để xem các đường đẳng thế, chọn
Contour plot từ hộp thoại Plot Selection.
d. Từ trường tĩnh: Nam châm, động cơ điện và máy biến biến áp là
những lĩnh vực có từ trường tĩnh. Phương trình Maxwell đối với trường tĩnh
là:
175
r r
∇× B = J
r
∇.B = 0
r
r
và quan hệ B = µH
r
Do ∇.B = 0 nên tồn tại từ thế vec tơ A sao cho B = ∇×A và:
r
r
⎛ 1
⎞
∇× ∇× A = J
⎜
⎟
µ
⎝
⎠
Do ∇.B = 0 nên tồn tại từ thế vec tơ A sao cho:
B = ∇A
r
r
⎛ 1
⎞
Và: ∇× ∇× A = J
⎜
⎟
µ
⎝
⎠
Trong trường hợp bài toán 2‐D ta coi dòng điện chạy theo hướng z và như vậy
vec tơ A chỉ có thành phần theo trục z:
A = (0,0,A),J = (0,0,J),
Và phương trình trên được đơn giản hoá thành phương trình ellipptic:
⎛ 1
⎞
− ∇. ∇A = J
⎜
⎟
µ
⎝
⎠
Trong đó : J = J(x,y).
Với bài toán 2‐D ta tính từ cảm B theo:
r
B =
⎛
⎞
∂A ∂A
∂y ∂x
⎜
⎜
⎟
⎟
,−
,0
⎝
⎠
và H được tính theo:
r
r
1
µ
H = B
Điều kiện biên trên mặt phân cách hai môi trường khác nhau là H×n phải liên
1 ∂A
tục,nghĩa là
liên tục. Điều này không đòi hỏi sự xử lí đặc biệt vì ta dùng
µ ∂n
công thức biến phân của bài toán PDE.
Trong vật liệu sắt từ, µ thường phụ thuộc vào B và do đó ta phải dùng
phương pháp giải phi tuyến. Điều kiện biên Dirichlet mô tả giá trị của từ thế A
trên biên. Điều kiện biên Neumann mô tả giá trị của thành phần pháp tuyến
r
⎛ 1
⎞
n
∇A . Điều này tương đương với việc mô tả giá trị tiếp tuyến của H trên
⎜
⎟
µ
⎝
⎠
biên. Ta khảo sát từ trường tĩnh tạo bởi cuộn dây stator của một động cơ điện
một chiều hai cực. Ta coi động cơ rất dài nên có thể bỏ qua hiệu ứng biên và
mô hình tính 2‐D là đủ. Ta tính trường trong 4 vùng:
176
•
•
•
hai vùng sắt từ là stator và rotor
khe hở không khí giữa stator và rotor
cuộn dây phần ứng có điên một chiều.
Độ từ thẩm µ = 1 trong cuộn dây và trong khe hở không khí. Trong stator và
µmax
rotor ta có: µ =
+ µmin
1+ c ∇(A) 2
µmax = 5000, µmin = 200 và c = 0.05 là các giá trị thường dùng đối với thép làm lõi
máy điện. Mật độ dòng điện J = 0 tại mọi nơi trừ cuộn dây. Trong cuộn dâ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
r
⎛ 1
⎞
vùng x ≥ 0, y ≥ 0 với điều kiện biên Neumann n ∇A 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:
r
r
∂H
∇×E = −µ
∂t
r
r
r
∂E
∂t
∇×H = ε
+ J
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à εµ :
r
∂2H
r
∆E − µε
= 0
∂t2
r
∂2E
r
∆H − µε
= 0
∂t2
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à:
r
r
J = σE và các sóng bị tắt dần do điện trở. Phương trình đối với E là:
r
r
∂2E
r
∂E
∂t
∆E − µσ − µε
= 0
∂t2
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 Ecejωt .
ω
t
Trường hợp bài toán phẳng ta có Ec = (0,0,Ec) và J = (0,0,Jej ) và từ trường là:
r
r
1
jωµ
H = (Hx ,Hy ,0) = −
∇×Ec
Phương trình vô hướng đối với Ec trở thành:
⎛ 1
µ
⎞
− ∇ ∇E + (jωσ − ω2ε)E = 0
⎜
⎝
⎟
⎠
c
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à:
εc = ε − jσ ω. Đ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:
j ⎛ 1
⎞
H = n ∇Ec
⎜
⎟
t
ω
µ
⎝
⎠
Các đại lượng có thể tính từ nghiệm là:
r
r
j
ω
B = ∇×E
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:
J0(kr)
J0(kR)
J = JS
Trong đó k = jωµε , 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:
− ∇(σ∇U) = 0
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:
− ∇(σ∇U) = 0
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:
∂T
ρC
− ∇(k∇T) = Q + h(Text − T)
∂t
Đâ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
15 trang
21560
Free
Bạn đang xem tài liệu "Giáo trình Matlab - Chương 8: Phương trình vi phân đạo hàm riêng (Phần 2)", để tải tài liệu gốc về máy hãy click vào nút Download ở trên
File đính kèm:
- giao_trinh_matlab_chuong_8_phuong_trinh_vi_phan_dao_ham_rien.pdf
