Дискретно-полевые модели электрических машин Часть 1. Численные метод

£/(!,/) = ц2(/).

Такие задачи также решают конечно-разностными методами, используя при этом безусловно устойчивые неявные схемы [24, 25, 27]. Для аппроксимации дифференциальных операторов возможны два варианта:

1. Линейный вариант, когда коэффициенты уравнения рас­ считываются по ранее найденным значениям искомых функций

к , т -

2. Нелинейный вариант, когда для определения коэффициен­

Линейный вариант реализуется значительно проще, так как значения функций на предыдущем временном слое известны. Нелинейный вариант сложнее в реализации. Значения функций, через которые определяются коэффициенты уравнения, неиз­ вестны. Поэтому приходится использовать итерационный про­ цесс. При этом значения коэффициентов задают приблизитель­ но, ориентируясь на предыдущий опыт в начале счета или ис­ пользуя данные с предыдущего временного интервала. Решают уравнение с принятыми коэффициентами и по найденным зна­ чениям функций уточняют величины коэффициентов. Таким образом производят несколько итераций, пока не будет достиг­ нута необходимая точность.

Если зависимость K(U) имеет монотонный характер, то ите­

рационный процесс сходится быстро, за несколько итераций. Если же эта зависимость имеет экстремум, то процесс сходится очень плохо или даже зацикливается. В этом случае для реше­ ния системы алгебраических уравнений часто используют метод Ньютона и его модификации, а для уточнения значений коэф­ фициентов уравнения применяют алгоритм Вегстейна [30, 31 ].

4. РЕШЕНИЕ МНОГОМЕРНЫХ УРАВНЕНИЙ ПАРАБОЛИЧЕСКОГО ТИПА

4.1.Решение многомерных нестационарных краевых задач

сиспользованием явных и неявных схем

До сих пор рассматривались дифференциальные уравнения с одной пространственной координатой. Многие задачи электро­ динамики многомерны, т.е. исследуемая область является функ­ цией нескольких пространственных координат.

Решение подобных задач традиционными методами связано со значительными временными затратами, так как число неиз­ вестных, а следовательно, число уравнений определяется произ­ ведением числа разбиений исследуемой функции по каждой пространственной координате. Основным требованием к мето­ дам решения многомерных краевых задач является их эконо­ мичность. Под экономичными методами понимаются методы решения краевых задач, при которых число математических операций, затрачиваемых на нахождение одного неизвестно­ го, не зависит от числа неизвестных [24, 27, 33]. Существует ряд способов решения многомерных уравнений параболиче­ ского типа.

Для решения многомерных нестационарных краевых задач могут быть использованы те же численные методы, что и для решения одномерных задач. Рассмотрим для примера решение нестационарной краевой двумерной задачи

ния к{=к2= к , N) = N2 = N , hx = hy = h , то число операций для определения значений функций на каждом временном слое ока­

жется пропорциональным N 2 Явная схема устойчива при вы­ полнении условия [24]

(4.4)

Тогда для достижения конечного времени счета Т потребу­ ется затратить Tjt а дг2 математических операций. Общее ко­ личество операций, затрачиваемых на решение краевой задачи, оказывается пропорциональным N 4 Реализация явных схем отличается предельной простотой и легко программируется.

Пример 4.1. Решение двумерной см еш анной краевой задачи с использованием неявной схемы. В квадрате [0:1,0:1] решим смешанную краевую задачу

при нулевых начальных и граничных условиях

U(x,y,0) = 0; С/(0,у,0 = 1/(1,у,/) = С/(х,0,0 = U(x,\,t) = 0.

Для решения задачи используем явную схему. Пространственные координаты х и у разбиваем на АС = 16

интервалов величиной h = \ / N =0,0625. Исходя из условия ус­ тойчивости решения

х < i f i —1 =0,00097656,

2 \ x h 2)

выбираем величину временного интервала х = 0,0009. Диффе­ ренциальные операторы аппроксимируем конечно-разностными выражениями и получаем систему алгебраических уравнений, записываемых в виде

Рис. 4.1. Результат решения двумерной краевой задачи: а - при выпол­ нении условия устойчивости (т = 0,0009); б - при нарушении условия устойчивости (т = 0,001)

На рис. 4.1 представлено решение смешанной краевой зада­ чи, выполненной по явной схеме для устойчивого варианта, когда временной интервал соответствует условию устойчиво­ сти, и неустойчивого варианта, когда условие устойчивости не выполнено.

4.1.2. Неявная схема

Для неявной схемы дифференциальное уравнение записыва­ ется следующим образом:

для тех же значений i , j . В этом случае на каждом временном

слое необходимо решать систему алгебраических уравне-

ний. Положим для упрощения N x = N 7 = N Тогда число неиз­

вестных будет равно N2. Матрица системы является разряжен­

ной, причем отличные от нуля элементы расположены лишь на пяти диагоналях. Для решения такой системы с учетом слабого

заполнения ленточной матрицы потребуется * №* операций.

Для неявныхсхем временной шаг можно выбирать произвольно. Поэтому для достижения конечного времени счета потребуется выполнить Т/х шагов по времени. В этом случае время решения

Т

краевой задачи будет пропорционально — S * •

Пример 4.2. Решение двумерной см еш анной краевой задача с использованием неявной схемы. В квадрате [0:1,0:1] решим

смешанную краевую задачу

при нулевых начальных и граничных условиях

Щх,у,0) = 0 ; U{0,y,t) = 1/(1,У,О = t/(*,0,/) = U ( x X t ) = 0 .

Длярешения задачи используем неявную схему.

Программа решения двумерной краевой задачи с использо­

ванием неявной схемы:

nl=17;n2=17; dn=16; hx=l./(nl-l); hy=lДп2-1);dt=0.001; nk=500;y(l:nk)=0; z(l:nk)=0; kl=dt/(hx*hx); k2=l.+4.*kl; a(l:225,l:225)=0.; f(l:225)=0.; x0(l:225)=0; f(97:99)=-100.*dt; f(l 12:114)=-100.*dt; Д127:129)=-100.*А; t=0.;n=0; Ю=0.;

il=17; i2=29; while il<212

for i=il :i2

a(i,i-15)=kl;a(M-l)=kl; a(i,i)=-k2; a(i,i+l)=kl;

end

for i=30:15:210

a(i,i-15)=kl; a(i,i-l)=kl; a(i,i)=-k2; a(i,i+15)=kl;

a( 1,1)=-k2; a( 1,2)=k1; a( 1,16)=k1; a(15,15)=-k2; a(15,14)=k1; a( 15,30)=k 1; a(211,212)=k 1; a(211,211 )=-k2; a(211,196)=k1; a(225,225)=-k2;

a(225,224)=kl; a(225,210)=kl; while n<nk

fk( 1:225)=-x0( 1:225)+f( 1:225); x=fk/a; x0(l:225)=x( 1:225); t=tO+dt*n; n=n+l; y(n)=x(l 13);

z(n)=n; disp(n); disp(x(l 13)); end

disp(x);

plot(z,y);

ходящимся (см. рис. 4.1,6).