Подпрограмма FCN для вычисления правых частей ДУ. Явная СДУ

SUBROUTINE FCN (N1, Tl, Yl, YPRIME).

USE MSIMSL; USE DESCRIBE IMPLICIT NONE.

INTEGER nl; REAL (8) Tl, Yl (nl), YPRIME (nl) if (dvigen) then.

x=Yl (5); CALL CalculatingSpline! Выч-е Ll, Rl, vzi, dvzdx, LD, SoprD d=Ll*Ld-vzi**2; endif.

il=(Yl (l)*Ld-Yl (2)*vzi)/d; i2=(Yl (2)*Ll-Yl (l)*vzi)/d !расч.токов no (1.32).

! Вычисление правых частей диф. уравнений:

YPRIME (l)=-il*Rl-Y (3); YPRIME (2)=-i2*SoprD; YPRIME (3)=il/emk if (.not.dvigen) RETURN.

if (uch==l) Fmex=FmO+zl*x; if (uch==2) Fmex=Fml+z2*(x-xl).

Fem=il*i2*dvzdx; YPRIME (4)=(Fem-Fmex)/massPodv; YPRIME (5)=Y1(4);

END.

Подпрограммы ParameterJDM, C_serialRL, FCN2. Параметры ИДМ

Вызов: CALL ParameterJDM SUBROUTINE ParameterJDM.

!Пп рассчитывает R, L, M параметры ИДМ для двухконтурной схемы замещения по Заданным размерам катушки и диска и расстоянию между ними. Должны быть заданы! ёмкость и напряжение питающего конденсатора, длина, сечение и материал соединительных проводов. Учитываются сопр. и инд. паразитные USE DESCRIBE; USE INCL_GRF; USE MSIMSL; IMPLICIT NONE EXTERNAL FCN2.

if (i==l) then; hKrach=tWKz;hDrach=hDisk; endif! Нач. высота кат. и диска perechK=.false.;perechD=.false.; nn=0; CRs=0d0; CLs=0d0 WRITE (out," (' hKat f DltK SoprK Lk hDisk DltD SoprD Ld vzi RnarK &.

&ksv'/' мм Гц мм мОм мкГн мм мм mkOm нГн мкГн mm')") RnarK=rKatvn+nWk*(twKr+tlzwwr)-tlzwwr [Наружный радиус катушки по металлу Rprov=roprov*(DIProv+DIProv)*4dO/(pi*DiamPr**2) [Активное сопротивл. проводов LProv=4d-7*DIProv*(LOG ((DiamPr+2dO*TizPr)/(DiamPr/2dO))+0.25dO)!MHfl. проводов CALL C_serialRL (emk, CRs, CLs, nn) [Вычислен.собственных R, L параметров конд-ра [Выдел.памяти масс., кот. хранят нач. и кон. номера площадок:

IF (ALLOCATEDfNACHW)) DEALLOCATE (NACHW, KONW).

ALLOCATE (NACHW (nwK), KONW (nwK)).

kollter=0.

3 CALL RAZBIVKA (hKrach, twKr, shrzbK, shWKz, shWKr, nrWKz, nrWKr, nrWK)!Pa36.витка Кат. nrK=nrWK*nWK [Число разбиений Катушки (полное число элементарн. витков в кат.).

CALL RAZBIVKA (hDrach, rDnar-rDvn, shrzbD, shDz, shDr, nrDz, nrDr, nrD)!pa6nBKa Диска ndimrz=nrD+nrK; NACHW (1)=1; KONW (l)=nrWK [размерность массивов RAD и Z rWvn=rKatvn; rWnar=rWvn+tWKr [радиусы внутр. и наружный для Витка Катушки SoprK=0d0 [Обнуление суммы для расчёта в цикле сопрот. кат-ки; заполнен. масс-ов do k=l, nwK; SoprK=SoprK+ldO/LOG (rWnar/rWvn); rWvn=rWnar+tlzWWr; rWnar=rWvn+tWKr if (k>l) then; NACHW (k)=NACHW (k-l)+nrWK; KONW (k)=KONW (k-l)+nrWK; endif enddo; SoprK=SoprK*piu2*udSoprK/hKrach [Сопротивление кат-ки RdPer=udSoprD*piu2/(hDrach*LOG (rDnar/rDvn)) [Сопротивление диска IF (ALLOCATED (RAD)) DEALLOCATE (RAD); IF (ALLOCATED (Z)) DEALLOCATE (Z);

ALLOCATE (RAD (ndimrz), Z (ndimrz)) !Выделение памяти для массивов RAD и Z RAD (l)=rKatvn+shWKr/2dO;Z (l)=shWKz/2dO^{г-коорд.центра} элемент.площ.в кат. NACHWD (l)=nrK+l; KONWD (l)=ndimrz !нач.и кон. номера площадок в диске! Расчёт координат центров площадок в кат. и заполн. ими массивов RAD и Z:

CALL KRD_PLOSK_KAT (RAD, Z, ndimrz, l, nrWKr, nrWKz, shWKz, shWKr, tlzwwr, otl, nrK) RAD (nrk+l)=rDvn+shDr/2dO;Z (nrk+l)=hKrach+sO+shDz/2dO+x !R, Z-KOop.l-ti площ. в диске CALL KRD_PLOSK_KAT (RAD, Z, ndimrz, nrK+l, nrDr, nrDz, shDz, shDr, OdO, otl, nrD) CALLIndKat_plosk_q (RAD_;Z, ndimrz, nwK, NACHW, KONW, nrWKr_;nrWKz, shWKr, shWKz, LkPer) CALL lndMWtk_qq (RAD, Z, ndimrz, nrDr, nrDz, nrK+l, ndimrz, shDr, shDz, LdPer).

CALL VZI_KK (RAD, Z, ndimrz, nwK, l, NACHW, KONW_;NACHWD, KONWD, vziPer); ksv=vziPer/SQRT (LkPer*LdPer) !коэф.магнитной связи катушки и диска RlPer=SoprK+Rprov+RTiristor !Активное сопротивление 1-го контура LlPer=LkPer+LTiristor+LProv !Индуктивность 1-го контура! Присвоение начальных значений переменным для решения СДУ (1.31):

Y=OdO; Y (3)=UcO; t=OdO; tend=OdO; k=l; XARR (1)=0.; YARR (3,l)=UcO; ido=l dPer=LlPer*LdPer-vziPer*vziPer IdPerДетерминант САУ (1.19), см.(1.32).

!Бесконечный цикл по времени для решения системы диф. ур-й на каждом шаге:

DO; k=k+l !к-Счётчик шагов для заполн. массивов XARR, YARR для постр. графиков tend=tend+ht !Изменение текущего времени CALL DIVPRK (ido, 3, FCN2,t, tend, tol, PARAM, Y) !решение СДУ (1.31).

il=(Y (l)*LdPer-Y (2)*vziPer)/dPer; i2=(Y (2)*LlPer-Y (l)*vziPer)/dPer !Реш.(1.32) относит. токов !Запоминание переменных на каждом шаге по времени для построения графиков: XARR (k)=t; YARR (l, k)=il; YARR (2,k)=i2; uc=Y (3); YARR (3,k)=uc !Условие для вылавливания первого перехода через ноль тока il катушки: if (il>OdO.AND.ipred<=OdO) then.

tO=(ipred*t-il*tpred)/(ipred-il); Tperiod=2d0*t0; f=0.5d0/t0; omga=piu2*f exit endif.

ipred=il; tpred=t !Запоминание тока и времени на каждом шаге для использования их в условии опред. перехода тока il через ноль if (t>=tkon) exit! Условие выхода из цикла по времени ENDDO;

ido=3; CALL DIVPRKfido, 3, FCN2,t, tend, tol, PARAM, Y) Остановка работы пп DIVPRK kollter=kollter+l; if (kollter>5)goto 5.

DeltK=SQRT (2dO*udSoprK/(omga*rmu));DeltD=SQRT (2dO*udSoprD/(omga*rmu))!3KB.my6. WRITE (out, 110) hKrach, fjDeltl^Soprl^LkPe^hDracl^DeltD^dPer, LdPer, vziPe^Rnarl^ksv 110 FORMAT (3pF5.2,0pF6.0,3pF5.2,F7.2,6pF6.1,3pF6.2,F5.2,6pF6.2,9pF6.2,6pF5.2 &, 3pF5.1,0pF5.2).

!Условие пересчёта R, L параметров диска по эквивал.глуб.проникн.поля: if (ABS (DeltD-hDrach)/hDrach>epsDel.AND.DeltD.

else; if (DeltD>=hDisk)then; perechD=.false.; hDrach=hDisk; endif endif.

!Условие пересчёта R, L параметров кат-ки по эквивал.глуб.проникн.поля: if (ABS (DeltK-hKrach)/hKrach>epsDel.AND.DeltK.

else; if (DeltK>=tWKz)then;perechK=.false.; hKrach=tWKz; endif endif.

!Условие пересчёта параметров катушки или диска: if (perechD.OR.perechK) then; perechD=.false.; perechK=.false.

!условие однократного уменьшения вдвое шагов разбразбиений катушки и диска:

! if (j==0)then; shrzbK=shrzbK/2dO; shrzbD=shrzbD/2dO; j=l;endif goto 3 {Освобождение массивов и переход к новому расчёту endif.

• 5 WRITE (out, 115) RlPer, LlPer, RdPer, LdPer, vziPer, x
• 115 FORMAT ('R, L, M параметры для решения ур-й (1.31),(1.32):'/'Rl=', 3pF8.3,^,MOM&

& И='_/6рРб.1_/'мкГн; SoprD=', 6pF7.3,'MKOM; Ld=', 9pF6.2,'Hl" H; vzi='_/6pF5.2 &.

'мкГн x='_/3pf4.2_/^,MM') !FDATA (1:)=L1; (2:)=R1; (3:)=M; (4:)=LD; (5:)=RD END SUBROUTINE ParameterJDM.

SUBROUTINE C_serialRL (emk, CRs, CLs, n).

!Пп вычисляет собственное активное сопротивление CRs и собственную индуктив- !ность CLs конденсатора по его ёмкости emk.

{Данные подпрограммы зависят от типа конденсатора, выбранного пользователем. {Здесь CRs аппроксимируется экспоненциальной зависимостью по справочным данным {конденсаторов типа DC 85-С, a CLs — линейной.

INTEGER: п.

REAL (8):emk, CRs, CLs, a, RO, b &.

Cl=280d-6,C2=3300d-6,Rl=4.1d-3,R2=l.ld-3,Ll=50d-9,L2=95d-9 !Из справочника SAVE a, RO, b if (n==0) then.

a=log (Rl/R2)/(Cl-C2); RO=Rl/EXP (a*Cl); b=(L2-Ll)/(C2-Cl); n=l endif.

CRs=RO*EXP (a*emk); CLs=Ll+b*(emk-Cl).

END SUBROUTINE C_serialRL.

SUBROUTINE FCN2 (nl, tl, Yl, YPR) {Подпрограмма расчёта правых частей СДУ (1.31).

USE DESCRIBE; IMPLICIT NONE; INTEGER nl REAL (8): Yl (nl), YPR (nl), tl il=(Yl (l)*LdPer-Yl (2)*vziPer)/dPer.

i2=(Yl (2)*LlPer-Yl (l)*vziPer)/dPer {расчёт токов по ур-ям (1.32).

{Вычисление правых частей системы диф. уравнений (1.31):

YPR (l)=-il*RlPer-Y (3); YPR (2)=-i2*RdPer; YPR (3)=il/emk END SUBROUTINE FCN2.

Подпрограмма CalculatingSpline

Вызов: CALL CalculatingSpline SUBROUTINE CalculatingSpline USE DESCRIBE; USE MSIMSL.

Ll=DCSVAL (x, nlntervX, BREAK (l:l:), CSCOEF (l:l:)) !Инд-ть катушки.

Rl=DCSVAL (x, nlntervX, BREAK (2:2:), CSCOEF (2:2:)) {Акт.сопр.катушки.

vzi=DCSVAL (x, nlntervX, BREAK (3:3:), CSCOEF (3:3:)) !Вз.инд.кат-ки и диска.

dvzdx=DCSDER (l, x, nlntervX, BREAK (3:3:), CSCOEF (3:3:)) {Производная dvzi/dx LD=DCSVAL (x, nlntervX, BREAK (4:4:), CSCOEF (4:4:)) !Инд-ть диска.

SoprD=DCSVAL (x, nlntervX, BREAK (5:5:), CSCOEF (5:5:)) {Акт.сопр.диска END.