Численное моделирование и анализ переходных процессов в электрической цепи. 3
Содержание:
- Введение.
- Постановка задачи, исходные
данные.
- Реализация поставленной задачи:
3.1. первый этап решения задачи: численная реализация системы дифференциальных уравнений:
3.1.1. алгоритмически,
построив блок-схему и
написанную на языке программирования Free Pascal,
используя алгоритм метода Рунге-Кутта ;
3.2. второй этап решения задачи:
3.2.1. решение задачи аппроксимации
используя язык
3.3. третий этап решения задачи: вычисление интеграла, определяющего количество теплоты, выделяемое на резисторе R4 за период времени t1<t<t2:
3.3.1. алгоритмически,
построив блок-схему и
написанную на языке программирования Free Pascal,
используя метод трапеций;
- Выводы по проделанной работе.
ВВЕДЕНИЕ
В ходе данной работы будут проанализированы переходные процессы в электрической цепи. Работа подразделена на 3 части:
I часть работы включает в себя решение системы дифференциальных уравнений, описывающих переходный процесс в электрической цепи.
II часть работы представляет из себя процесс решения задачи аппроксимации.
Заключительная III часть работы: численное интегрирование. Для выполнения данной курсовой работы необходимо применение не только знаний по информатике, но также и физике, математике.
Постановка задачи, исходные данные
Дана схема электрической цепи:
Параметры элементов цепи:
E=E0*sin(w*t+φ) – гармонический источник тока
E0=15 В – амплитуда колебаний
w=2П*f – циклическая частота
t– текущий момент времени
f=30 Гц – линейная частота, φ=9П/5 – (параметры заданные по вариантам)
R1=30 Ом, R2=25 Ом , R3=50 Ом , R4=1.88 Ом
R5=15 Ом, R6=50 Ом - резисторы
L=5*57 мГн - катушка индуктивности
С=20 мкФ - конденсатор
В начальный момент t=t0=0 ключ находится в положении 1. При этом цепь разомкнута, напряжение на конденсаторе и ток на катушке равны нулю (U=0, I=0). Происходит первое переключение ключа (ключ мгновенно переводится в положение 2). При этом происходит заряд конденсатора, меняются значения U и I. В момент времени t=t1=0.01c ключ мгновенно переводится в положение 1. Конденсатор разряжается, вновь меняются значения U и I. Анализ схемы заканчивается в момент t=t1=0.02c
Реализация поставленной задачи
Вывод системы дифференциальных уравнений.
В соответствии со схемой запишем соотношения из законов Ома и Кирхгофа
Ключ в положении “2”:
i=i1+i2
ключ в положении “1”:
i=i1+i2
Исключив токи i1 и i2 получили ключ в положении “2”:
Для положения ключа “1”:
т.е. E(t)=0
Для обоих положений ключа можно записать систему дифференциальных уравнений в виде:
Описание методов применяемых в курсовой работе.
Метод Рунге-Кутта:
Метод Рунге—Кутта 4 порядка столь широко распространён, что его часто называют просто методом Рунге—Кутта.
Рассмотрим задачу Коши
Тогда значение в следующей точке вычисляется по формуле:
где h — величина шага сетки по x и вычисление нового значение проходит в четыре этапа:
Этот метод
имеет четвёртый порядок
Метод трапеций (ЛП):
Если функцию на каждом из частичных отрезков аппроксимировать прямой, проходящей через конечные значения, то получим метод трапеций.
Площадь трапеции на каждом отрезке:
Погрешность аппроксимации на каждом отрезке:
Полная формула трапеций в случае деления всего промежутка интегрирования на отрезки одинаковой длины h:
Погрешность формулы трапеций:
Метод Гаусса для решения систем линейных уравнений:
Пусть исходная система выглядит следующим образом
Тогда согласно свойству элементарных преобразований над строками эту систему можно привести к трапециальному виду:
Переменные называются главными переменными. Все остальные называются свободными.
Если , то рассматриваемая система несовместна.
Предположим, что .
Перенесём свободные переменные за знаки равенств и поделим каждое из уравнений системы на свой коэффициент при самом левом (, где — номер строки):
где
Если свободным переменным системы (2) придавать все возможные значения и вычислить через них главные переменные, то мы получим все решения этой СЛАУ. Так как эта система получена путём элементарных преобразований над исходной системой (1), то по теореме об эквивалентности при элементарных преобразованиях полученное нами решение является решением системы (1).
Метод парабол (Симпсона):
Основан на интерполяции подынтегральной функции y=f(x) параболой на паре соседних частичных интервалов (i=1,3), т.е. интерполяционным полиномом второй степени
Тогда оказывается, что частичный интеграл для интерполяционного полинома вычисляется по формуле
Выбирая число частичных интервалов n четным, получим суммированием составную
формулу для приближенного вычисления всего интеграла методом парабол
Аппроксимация полиномом 1-ой степени:
меру близости искомой прямой к графику узловых точек выразим величиной:
Найдем а1 и а2 из условия минимума суммы квадратов отклонений:
Из этого условия
получаем систему линейных
полученная система называется системой нормальных уравнений метода наименьших квадратов и в матричной форме имеет вид:
Аппроксимация полиномом 2-ой степени
В матричной форме решение будет иметь вид:
А=В-1С
В
А С
Составим систему нормальных уравнений в матричной форме:
Решение системы дифференциальных уравнений на языке Free Pascal:
2. Решение для системы дифференциальных уравнений (3):
program cursc3;
var a,b,m,d,f,t,i,u,h,k1,k2,k3,k4,
begin
writeln(' vvod R1,R2,R3,R4,R5,R6');
READLN(R1,R2,R3,R4,R5,R6);
l:=5.57*0.001;
c:=20*0.000001;
A:=(R5+R6)/(R3+R5+R6);
B:=R1*R2/(R1+R2);
M:=1/(R3+R5+R6);
D:=R4+B+A*r3;
F:=R2/(R1+R2);
h:=0.00025;
i:=0;
u:=0;
t:=0;
while (t<=0.01001) do begin
k1:=h*(15*sin(62.8*t+7*pi/5)*
m1:=h*(i*a/c-m*u/c);
k2:=h*(15*sin(62.8*(t+h/2)+7*
m2:=h*((i+k1/2)*a/c-m*(u+m1/2)
k3:=h*(15*sin(62.8*(t+h/2)+7*
m3:=h*((i+k2/2)*a/c-m*(u+m2/2)
k4:=h*(15*sin(62.8*(t+h)+7*pi/
m4:=h*((i+k3)*a/c-m*(u+m3)/c);
i:=i+(k1+2*k2+2*k3+k4)/6;
u:=u+(m1+2*m2+2*m3+m4)/6;
t:=t*1000;
writeln(' pri t=',t:2:2,'*10E-3',' u=',u:2:2,' i=',i:4:4);
t:=t/1000;
t:=t+h;
end;
readln
end.
Результаты работы программы:
Возьмем начальные данные для решения системы дифференциальных уравнений (2) из результата решения системы дифференциальных уравнений (3).
1. Решение для системы дифференциальных уравнений (2):
program cursc;
var a,b,m,d,f,t,i,u,h,k1,k2,k3,k4,
begin
writeln(' vvod R1,R2,R3,R4,R5,R6');
READLN(R1,R2,R3,R4,R5,R6);
l:=5.57*0.001;
c:=20*0.000001;
A:=(R5+R6)/(R3+R5+R6);
B:=R1*R2/(R1+R2);
M:=1/(R3+R5+R6);
D:=R4+B+A*r3;
F:=R2/(R1+R2);
h:=0.00025;
i:=-0.08;
u:=-5.32;
t:=0.01+h;
while (t<=0.02001) do begin
k1:=h*(-d*i/l-a*u/l);
m1:=h*(i*a/c-m*u/c);
k2:=h*(-d*(i+k1/2)/l-a*(u+m1/
m2:=h*((i+k1/2)*a/c-m*(u+m1/2)
k3:=h*(-d*(i+k2/2)/l-a*(u+m2/
m3:=h*((i+k2/2)*a/c-m*(u+m2/2)
k4:=h*(-d*(i+k3)/l-a*(u+m3)/l)
m4:=h*((i+k3)*a/c-m*(u+m3)/c);
i:=i+(k1+2*k2+2*k3+k4)/6;
u:=u+(m1+2*m2+2*m3+m4)/6;
t:=t*1000;
writeln(' pri t=',t:2:2,'*10E-3',' u=',u:2:2,' i=',i:4:4);
t:=t/1000;
t:=t+h;
end;
readln
end.
Результаты работы программы:
В результате решения систем (2) и (3) получим таблицу результатов и графики функций I(t) и U(t)(Гра фики строю в пакете Excel);
t |
U |
I |
0 |
-0,65 |
-0,101 |
0,00025 |
-1,43 |
-0,123 |
0,0005 |
-2,13 |
-0,1232 |
0,00075 |
-2,73 |
-0,1181 |
0,001 |
-3,22 |
-0,1125 |
0,00125 |
-3,62 |
-0,1075 |
0,0015 |
-3,96 |
-0,1033 |
0,00175 |
-4,23 |
-0,1 |
0,002 |
-4,45 |
-0,0972 |
0,00225 |
-4,64 |
-0,095 |
0,0025 |
-4,79 |
-0,0932 |
0,00275 |
-4,91 |
-0,0917 |
0,003 |
-5,02 |
-0,0906 |
0,00325 |
-5,1 |
-0,0896 |
0,0035 |
-5,18 |
-0,0888 |
0,00375 |
-5,23 |
-0,0882 |
0,004 |
-5,28 |
-0,0876 |
0,00425 |
-5,32 |
-0,0871 |
0,0045 |
-5,36 |
-0,0867 |
0,00475 |
-5,39 |
-0,0861 |
0,005 |
-5,41 |
-0,0858 |
0,00525 |
-5,43 |
-0,0858 |
0,0055 |
-5,44 |
-0,0855 |
0,00575 |
-5,45 |
-0,0852 |
0,006 |
-5,46 |
-0,0849 |
0,00625 |
-5,46 |
-0,0846 |
0,0065 |
-5,47 |
-0,0846 |
0,00675 |
-5,47 |
-0,0844 |
0,007 |
-5,47 |
-0,0841 |
0,00725 |
-5,46 |
-0,0838 |
0,0075 |
-5,46 |
-0,0835 |
0,00775 |
-5,45 |
-0,0831 |
0,008 |
-5,44 |
-0,0828 |
0,00825 |
-5,43 |
-0,0824 |
0,0085 |
-5,42 |
-0,082 |
0,00875 |
-5,41 |
-0,0817 |
0,009 |
-5,39 |
-0,0812 |
0,00925 |
-5,38 |
-0,0808 |
0,0095 |
-5,36 |
-0,0804 |
0,00975 |
-5,34 |
-0,08 |
0,01 |
-5,32 |
-0,0794 |
0,01025 |
-4,65 |
0,0208 |
0,0105 |
-3,87 |
0,0422 |
0,01075 |
-3,16 |
0,0417 |
0,011 |
-2,57 |
0,036 |
0,01125 |
-2,09 |
0,0298 |
0,0115 |
-1,69 |
0,0243 |
0,01175 |
-1,37 |
0,0198 |
0,012 |
-1,11 |
0,016 |
0,01225 |
-0,9 |
0,013 |
0,0125 |
-0,73 |
0,0105 |
0,01275 |
-0,59 |
0,0085 |
0,013 |
-0,48 |
0,0069 |
0,01325 |
-0,39 |
0,0056 |
0,0135 |
-0,31 |
0,0045 |
0,01375 |
-0,25 |
0,0037 |
0,014 |
-0,21 |
0,003 |
0,01425 |
-0,17 |
0,0024 |
0,0145 |
-0,13 |
0,0019 |
0,01475 |
-0,11 |
0,0016 |
0,015 |
-0,09 |
0,0013 |
0,01525 |
-0,07 |
0,001 |
0,0155 |
-0,06 |
0,0007 |
0,01575 |
-0,05 |
0,0005 |
0,016 |
-0,04 |
0,0004 |
0,01625 |
-0,03 |
0,0004 |
0,0165 |
-0,02 |
0,0003 |
0,01675 |
-0,02 |
0,0002 |
0,017 |
-0,02 |
0,0002 |
0,01725 |
-0,01 |
0,0001 |
0,0175 |
-0,01 |
0,0001 |
0,01775 |
-0,01 |
0,0001 |
0,018 |
-0,01 |
0,0001 |
0,01825 |
0 |
0,0001 |
0,0185 |
0 |
0,0001 |
0,01875 |
0 |
0,0001 |
0,019 |
0 |
0 |
0,01925 |
0 |
0 |
0,0195 |
0 |
0 |
0,01975 |
0 |
0 |
0,02 |
0 |
0 |
Второй этап решения:
Таблица и график для интервала: 0≤Т≤0,005:
t |
I |
|
0 |
-0,101 |
0,00025 |
-0,123 |
0,0005 |
-0,1232 |
0,00075 |
-0,1181 |
0,001 |
-0,1125 |
0,00125 |
-0,1075 |
0,0015 |
-0,1033 |
0,00175 |
-0,1 |
0,002 |
-0,0972 |
0,00225 |
-0,095 |
0,0025 |
-0,0932 |
0,00275 |
-0,0917 |
0,003 |
-0,0906 |
0,00325 |
-0,0896 |
0,0035 |
-0,0888 |
0,00375 |
-0,0882 |
0,004 |
-0,0876 |
0,00425 |
-0,0871 |
0,0045 |
-0,0867 |
0,00475 |
-0,0861 |
0,005 |
-0,0858 |
Решение задачи аппроксимации и интерполяции на языке Free Pascal:
Разобьем участок на 2:
1: 0≤Т≤0,00075
2: 0,00075≤Т≤0,005
Найдем интерполяцией график квадратичной функции для первого участка, для этого необходимо решить систему линейных уравнений:
02a2+0а1+а0=-0,101,
0.000252a2+0.00025a1+a0=-0.
0.00052a2+0.0005a1+a0=-0.1232;
Здесь и далее решать системы линейных уравнений будем методом гаусса реализованным в Free Pascal(обозначим ее (*)):
Программа:
program laba3;
var c,d,s,t,z:real;
i,j,e,r,n:integer;
a:array [1..1000,1..1000] of real;
b:array [1..1000] of real;
x:array [1..1000] of real;
begin
repeat
writeln('vvod razmer');
readln(n);
for i:=1 to n do
for j:=1 to n do begin
write('vvod a[',i:1,',',j:1,']=');
readln(a[i,j]);
end;
for i:=1 to n do begin
write('vvod b[',i:1,']=');
readln(b[i]);
end;
for i:=1 to n do begin
for j:=1 to n do begin
write(' ',a[i,j]:6:2);
end;
writeln(' ',b[i]:6:2);
end;
for i:=1 to n do begin
d:=a[i,i];
for j:=1 to n do begin
a[i,j]:=a[i,j]/d;
end;
b[i]:=b[i]/d;
for r:=i+1 to n do begin;
c:=-a[r,i];
b[r]:=c*b[i]+b[r];
for e:=1 to n do begin
a[r,e]:=c*a[i,e]+a[r,e];
end; end; end;
writeln; writeln(' new matrica');
for i:=1 to n do begin
for j:=1 to n do begin
write(' ',a[i,j]:6:2);
end;
writeln(' ',b[i]:6:2);
end;
s:=0;
for i:=0 to n-1 do begin
x[n-i]:=b[n-i]+s;
s:=0;
for j:=n-i to n do begin
if (i<n-1) then s:=s-x[j]*a[n-i-1,j] else s:=5;
end; end;
writeln;
for i:=1 to n do begin
writeln(' x[',i:1,']=',x[i]:6:2);
end; writeln;
writeln(' vvod 1 esli vse, 2 - ehe reshat');
readln(z);
t:=2-z;
until(t>0.5);
writeln(' Do svidania!!!');
readln;
end.
Результат работы программы:
Уравнение для первого интервала имеет вид:
I=181600t2-137t-0,101
Аппроксимацией найдем функцию описывающую второй интервал с помощью программ реализованных в Free Pascal:
1)program apro;
var s,st,st2,st3,si,sit,sit2,st4:
k:integer;
i:array [1..19] of real;
t:array [0..19] of real;
begin
t[0]:=0.00025;
for k:=1 to 19 do begin
write(' vvod I[',k:1,']=');
read(i[k]);
t[k]:=t[k-1]+0.00025;
writeln(' t[',k:1,']= ',t[k]:6:5);
end;
s:=0;
st:=0;
st2:=0;
st3:=0;
si:=0;
sit:=0;
sit2:=0;
st4:=0;
for k:=1 to 19 do begin
s:=s+1;
st:=st+t[k];
st2:=st2+t[k]*t[k];
st3:=st3+t[k]*t[k]*t[k];
st4:=st4+t[k]*t[k]*t[k]*t[k];
si:=si+i[k];
sit:=sit+i[k]*t[k];
sit2:=sit2+i[k]*t[k]*t[k];
end;
writeln(' s=',s:3:1,' st=',st:5:4,' st2=',st2:10:10,' st3=',st3:10:10,
' st4',st4:15:15,' si=',si:5:4,' sit=',sit:5:4,' sit2=',sit2:10:10);
readln;
end.
Результат работы программы:
Получили расширенную матрицу:
19 0,0522 0,0001793125 -1,83
0,0522 0,0001793125 0,000000689 -0,0048
0,0001793125 0,000000689 0,00000000282291 -0,000016
С помощью программы (*) получим:
Уравнение, описывающее второй интервал, имеет вид:
I=-94t2+7,09t-0.155
Третий этап решения:
Численное интегрирование
Выполним численное
1) I=181600t2-137t-0,101;
2) I=-94t2+7,09t-0.115;
В результате интегрирования для двух интервалов получим два значения теплоты: Q1 и Q2.
Общее количество теплоты полученное на резисторе R4 получим сложением Q1 и Q2:
Q=Q1+Q2;
Интегрирование проведем методом трапеций, реализованным в Free Pascal:
Для первого интервала:
program integrate1;
var i,i1,s,t,h,tn,tk:REAL;
begin
writeln(' vvod tn,tk');
readln(tn,tk);
t:=tn;
h:=0.00025;
s:=0;
while (t<=tk-h/2) do begin
i:=sqr(181600*t*t-137*t-0.101)
i1:=sqr(181600*(t+h)*(t+h)-
s:=s+(i+i1)*h/2;
t:=t+h;
end;
s:=s*1.88;
writeln(' Q=',s:5:5);
readln;
end.
Результат работы программы:
Для второго интервала:
program integrate2;
var i,i1,s,t,h,tn,tk:REAL;
begin
writeln(' vvod tn,tk');
readln(tn,tk);
t:=tn;
h:=0.00025;
s:=0;
while (t<=tk-h/2) do begin
i:=sqr((-94)*t*t+7.09*t-0.155)
i1:=sqr((-94)*(t+h)*(t+h)+7.
s:=s+(i+i1)*h/2;
t:=t+h;
end;
s:=s*1.88;
writeln(' Q=',s:5:5);
readln;
end.
Результат работы программы:
Найдем суммарную теплоту:
Q=Q1+Q2=0,00011+0,0002=0,
Вывод по проделанной работе:
В данной курсовой работе предложено рассчитать параметры электрической цепи переменного тока. При помощи метода Рунге-Кутта я сначала решил систему дифференциальных уравнений (описывающих процессы, происходящие в цепи, которая была выведена в начале работы) для того, что бы найти зависимость тока (I(t)) и напряжения (U(t)) от времени на заданном временном интервале t1 < t <t2. Затем я произвел аппроксимацию полученных данных средствами Free Pascal. Впоследствии, проинтегрировав полученную при аппроксимации функцию на заданном интервале времени T1 < T <T2, нашел количество теплоты, выделившееся на резисторе R4. При решении интеграла мною был применен метод трапеций. При выполнении курсовой работы были построены соответствующие графики зависимостей, сравнив результаты вычислений представленными методами, были сделаны конструктивные выводы.