1. Ïðèìåð ïðîãðàììû ÌÊÝ
Äàííûé ïðèìåð ÿâëÿåòñÿ ïåðåâîäîì íà ðóññêèé ÿçûê ñîîòâåòñòâóþùåé ãëàâû êíèãè
A. Milenin Podstawy MES. Zagadnienia termomechaniczne, Wydawnictwo AGH, Kraków, 2010.
(À. Ìèëåíèí Îñíîâû ÌÊÝ. Òåðìîìåõàíè÷åñêèå çàäàíèÿ. Èçäàòåëüñòâî ÀÃÕ, Êðàêîâ, 2010)
www
1.1. Ôîðìóëèðîâêà ïðîáëåìû
Ðàññìîòðèì çàäà÷ó íåñòàöèîíàðíîãî òåïëîîáìåíà. Ñèñòåìó óðàâíåíèé ÌÊÝ ïîëó÷àåì ñ ïîìîùüþ ñëåäóþùåãî âûðàæåíèÿ:
, (6.8)
ãäå:
(5.14)
, (5.15)
. (6.3)
 óðàâíåíèÿõ (6.8) {t1} - óçëîâûå çíà÷åíèÿ òåìïåðàòóðû â ìîìåíò âðåìåíè
, {t0} - óçëîâûå
çíà÷åíèÿ òåìïåðàòóðû â ìîìåíò âðåìåíè 
.
1.2. èñõîäíûå äàííûå
Ïðèìåð èñõîäíûõ äàííûõ íà ðèñóíêå 9.1.
100 mTbegin Íà÷àëüíàÿ òåìïåðàòóðà, C
1000 mTime Âðåìÿ ïðîöåññà, s
1 mdTime Øàã ïî âðåìåíè, s
1200 mT_otoczenia Òåìïåðàòóðà îêðóæàþùåé ñðåäû C
300 mAlfa Êîýôôèöèåíò òåïëîîáìåíà W/Cm2
100 mH0 Âûñîòà ñå÷åíèÿ, mm
100 mB0 Øèðèíà ñå÷åíèÿ, mm
25 mNhH ×èñëî óçëîâ ïî âûñîòå
25 mNhH ×èñëî óçëîâ ïî øèðèíå
700 mC Òåïëîåìêîñòü, J/Ckg
25 mK Òåïëîïðîâîäíîñòü
7800 mR Ïëîòíîñòü kg/m3
Rys. 9.1. Ïðèìåð ñîäåðæàíèÿ ôàéëà ñ èñõîäíûìè äàííûìè
1.3. Òåêñò ïðîãðàììû
Ïðîãðàììà íàïèñàíà íà ÿçûêå FORTRAN90. Èñïîëüçîâàí êîìïèëÿòîð Intel. Ïîëíûé òåêñò ìîæíî ñêà÷àòü ñ àäðåñà http://home.agh.edu.pl/~milenin/programs/.
1.3.1. Ãëàâíàÿ ïðîãðàììà TEMP2D
Êîä ãëàâíîé ïðîãðàììû TEMP2D ïîêàçàí íà rys. 9.2. Îáîçíà÷åíèÿ ãëàâíûõ ïåðåìåííûõ ñëåäóþùèå:
GlobData – Íàçâàíèå ìîäóëÿ ñ ãëàâíûìè (ãëîáàëüíûìè) ïåðåìåííûìè;
dTauMax – Ìàêñèìàëüíî äîïóñòèìûé øàã ïî âðåìåíè;
n, Ntau – Íîìåð øàãà ïî âðåìåíè è ÷èñëî ýòõ øàãîâ;
Ãëàâíûå ïðîöåäóðû:
IniEL4 – èíèöèàëèçàöèÿ äàííûõ äëÿ âûáðàííîãî êîíå÷íîãî ýëåìåíòà (âûáðàííûé òèï: 4- óçëîâîé ýëåìåíò);
InpData – ñ÷èòûâàíèå èñõîäíûõ äàííûõ;
GenGrid2d – ãåíåðàöèÿ ñåòêè ÊÝ;
ALLOCATE_Matrix – âûäåëåíèå ïàìÿòè äëÿ ìàññèâîâ;
SOLVER – ðåøåíèå äëÿ òåêóøåãî øàãà ïî âðåìåíè.
program TEMP2D
!*************************************************************************
! PROGRAM: TEMP2D
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
!************************************************************************
use GlobData;
implicit none
real*8 :: Asr,dTauMax,TauP;
integer*4 :: n, Ntau, i, iErr;
call IniEL4();
call InpData();
call GenGrid2d(mH0,mB0,mNhH,mNhB, mGr);
call SetControlPoints();
CALL ALLOCATE_Matrix();
CALL WriteControlPointsBegin();
Asr = mK / (mC * mR);
mdTime = DBLE(mB0/(1e3*mNhB))**2/(0.5*Asr);
CALL WriteControlPoints();
Ntau = int(mTime/mdTime)+1;
mdTime = mTime/DBLE(Ntau);
mTau = 0;
do n=1, Ntau
mTau = mTau + mdTime;
CALL SOLVER( );
CALL WriteControlPoints();
end do;
CLOSE(88); CLOSE(89);
DEALLOCATE ( mA, STAT=iErr );
DEALLOCATE ( mB, STAT=iErr );
DEALLOCATE ( mX, STAT=iErr );
end program TEMP2D
Rys. 9.2. Ãëàâíàÿ ïðîãðàììà TEMP2D
1.3.2. Ñòðóêòóðû äàííûõ è ãëîáàëüíûå ïåðåìåííûå
Ïðàâèëüíî ñïðîåêòèðîâàííûå ñòðóêòóðû äàííûõ ñíèæàþò âåðîÿòíîñòü îøèáêè è óïðîùàþò ðàñ÷åò. Íèæå ïðèâåäåí ïðèìåð òàêèõ ñòðóêòóð.
Íà rys. 9.3 ïîêàçàí êîä ìîäóëÿ My_typ ñ îïðåäåëåíèåì ñòðóêòóð äàííûõ. Îïèñàíû ñëåäóþùèå ñòðóêòóðû:
Node – ñòðóêòóðà äàííûõ, ñîäåðæàùàÿ èíôîðìàöèþ î óçëàõ ÊÝ ñåòêè: x, y – êîîðäèíàòû óçëà, t – òåìïåðàòóðà â óçëå, CR – ñêîðîñòü èçìåðåíèÿ òåìïåðàòóðû;
Element – ñòðóêòóðà äàííûõ, ñîäåðæàùàÿ èíôîðìàöèþ îá ýëåìåíòå ñåòêè ÌÊÝ: nop – íîìåðà óçëîâ, êîòîðûå âõîäÿò â ñîñòàâ ýëåìåíòà, Npov – ÷èñëî ïîâåðõíîñòåé ÊÝ, êîòîðûå êîíòàêòèðóò ñ îêðóæàþùåé ñðåäîé; aPov – ëîêàëüíûå íîìåðà êîíòàêòíûõ ïîâåðõíîñòåé ýëåìåíòà.
gr2d – ñòðóêòóðà, ñîäåðæàùàÿ ñåòêó ÊÝ: nh – ÷èñëî óçëîâ; ne – ÷èñëî ÊÝ; nbn – ÷èñëî óçëîâ â èñïîëüçóåìîì ÊÝ; ncn – ÷èñëî ñòåïåíåé ñâîáîäû â îäíîì ÊÝ, äëÿ òåïëîâîé çàäà÷è ðàâíà nbn; nhPov – ÷èñëî óçëîâ ñåòêè ÊÝ, íàõîäÿùèõñÿ íà ïîâåðõíîñòè êîíòàêòà ñ îêðóæàþùåé ñðåäîé; EL – ìàññèâ äëèíû ne ñîäåðæàùèé äàííûå òèïà typu Element äëÿ êàæäîãî ÊÝ ñåòêè; ND – ìàññèâ äëèíîé nh ñîäåðæàùèé äàííûå òèïà Node äëÿ êàæäîãî óçëà ñåòêè.
cor_L – ñòðóêòóðà äàííûõ, ñîäåðæàùàÿ ëîêàëüíûå êîîðäèíàòû N, E òî÷êè â îáúåìå ýëåìåíòà.
Gran – ñòðóêòóðà äàííûõ, ñîäåðæàùàÿ èíôîðìàöèþ î ïîâåðõíîñòÿõ èñïîëüçóåìîãî êîíå÷íîãî ýëåìåíòà: N_p – ÷èñëî òî÷åê èíòåãðèðîâàíèÿ ïî ïîâåðõíîñòè; P – ìàññèâ ñòðóêòóð òèïà cor_L äëèíîé N_p ñîäåðæàùèé ëîêàëüíûå êîîðäèíàòû òî÷åê èíòåãðèðîâàíèÿ Ãàóññà; W – ìàññèâ äëèíîé N_p ñîäåðæàùèé âåñà òî÷åê èíòåãðèðîâàíèÿ â ñîîòâåòñòâèè; Nf – ìàññèâ äëèíîé nbn ñîäåðæàùèé çíà÷åíèÿ ôóíêöèé ôîðìû, ðàñ÷èòàííûå â òî÷êàõ èíòåãðèðîâàíèÿ íà ïîâåðõíîñòè ýëåìåíòà; N1,N2 – ìàññèâû äëèíîé nbn ñîäåðæàùèå çíà÷åíèÿ ïðîèçâîäíûõ ôóíêöèé ôîðìû ïî ëîêàëüíûì ïåðåìåííûì, ðàñ÷èòàííûå â òî÷êàõ èíòåãðèðîâàíèÿ íà ïîâåðõíîñòè ýëåìåíòà; UZEL – ëîêàëüíûå íîìåðà óçëîâ ïîâåðõíîñòè.
ELEM – ñòðóêòóðà äàííûõ, ñîäåðæàùàÿ äàííûå îá èñïîëüçóåìîì êîíå÷íîì ýëåìåíòå: nbn – ÷èñëî óçëîâ â ýëåìåíòå, nbnp – ÷èñëî óçëîâ, èñïîëüçóåìûõ äëÿ èíòåðïîëÿöèè ñðåäíåãî íàïðÿæåíèÿ (èñïîëüçóåòñÿ â çàäà÷å òåîðèè ïëàñòè÷íîñòè), N_p – ÷èñëî òî÷åê èíòåãðèðîâàíèÿ â îáúåìå ýëåìåíòà; Nf – ìàññèâ ðçìåðîì nbn ñîäåðæàùèé çíà÷åíèÿ ôóíêöèé ôîðìû â òî÷êàõ èíðåãðèðîâàíèÿ â îáúåìå ýëåìåíòà; N1, N2 – ìàññèâû äëèíîé nbn ñîäåðæàùèå çíà÷åíèÿ ïðîèçâîäíûõ ôóíêöèé ôîðìû ïî ëîêàëüíûì ïåðåìåííûì, ðàñ÷èòàííûå â òî÷êàõ èíòåãðèðîâàíèÿ â îáúåìå êîíå÷íîãî ýëåìåíòà; P – ìàññèâ ñòðóêòóð òèïà cor_L äëèíîé N_p ñîäåðæàùèé êîîðäèíàòû òî÷åê èíòåãðèðîâàíèÿ Ãàóññà â îáúåìå ýëåìåíòà â ëîêàëüíîé ñèñòåìå êîîðäèíàò; L – ìàññèâ ñòðóêòóð òèïà cor_L äëèíîé N_p ñîäåðæàùèé ëîêàëüíûå êîîðäèíàòû óçëîâ ýëåìåíòà; W – ìàññèâ äëèíîé N_p ñîäåðæàùèé âåñà òî÷åê èíòåãðèðîâàíèÿ Ãàóññà â îáúåìå ýëåìåíòà; Sf – ìàññèâ äëèíîé, ðàâíîé ìàêñèìàëüíîìó êîëè÷åñòâó ïîâåðõíîñòåé èñïîëüçóåìîãî êîíå÷íîãî ýëåìåíòà (â ïðèìåðå ðàâíà 4), ñîäåðæàùèé ïåðåìåííûå òèïà Gran.
module My_typ;
! ******************************************************************
! Typy danych
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ******************************************************************
implicit none;
type node
real*8 :: x,y,t,CR;
integer*4 :: status;
end type
type element
integer*4 :: nop(4);
integer*4 :: Npov; ! Liczba kontaktowych powierzchni elementu 0-1-2
integer*4 :: aPov(2); ! Lokalne numery kontaktowych powierzchni elementu
end type
type gr2d
integer*4 :: nh, ne, nbn, ncn, nhPov;
type(element) , dimension(:), pointer :: EL;
type(node) , dimension(:), pointer :: ND;
end type
type cor_L;
real*8 N, E;
end type cor_L;
type Gran
integer*4 :: N_p; ! Liczba punktów całkowania po powierzchni grani
type(cor_L), dimension(:), pointer :: P;
real*8, dimension(:), pointer :: W;
real*8, dimension(:,:), pointer :: N1,N2,Nf; ! N1(nbn,N_p)
integer*1, dimension(:),pointer :: UZEL;
! Lokalne numery węzłów grani
end type Gran;
type ELEM;
integer*4 :: nbn, nbnp, N_p;
real*8, dimension(:,:), pointer :: N1,N2,Nf; ! N1(nbn,N_p)
type(cor_L), dimension(:), pointer :: P; ! P(N_p)
type(cor_L), dimension(:), pointer :: L; ! Współrzędne węzłów elementu
real*8, dimension(:), pointer :: W; ! Waga punktów całkowania W(N_p)
type(Gran) :: Sf(4);
end type ELEM;
end module My_typ;
Rys. 9.3. Èñïîëüçóåìûå òèïû äàííûõ
Ìîäóëü GlobData ñîäåðæèò îïðåäåëåíèå ãëîáàëüíûõ ïåðåìåííûõ, äîñòóï ê ýòèì ïåðåìåííûì: use GlobData. Ñîäåðæàíèå ìîäóëÿ íà rys. 9.4.
module GlobData;
! ******************************************************************
! Globalne zmianne
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ******************************************************************
use My_typ;
implicit none;
real*8 :: mTbegin; ! Íà÷àëüíàÿ òåìïåðàòóðà
real*8 :: mTime; ! Âðåìÿ ïðîöåììà
real*8 :: mdTime; ! Øàã ïî âðåìåíè
real*8 :: mTau; ! Òåêóùåå âðåìÿ
real*8 :: mT_otoczenia; ! Òåìïåðàòóðà îêðóæàþùåé ñðåäû
real*8 :: mAlfa; ! Êîýôôèöèåíò òåïëîîáìåíà
real*8 :: mH0; ! Âûñîòà ñå÷åíèÿ
real*8 :: mB0; ! Øèðèíà ñå÷åíèÿ
integer*4 :: mNhH; ! ×èñëî óçëîâ ïî âûñîòå
integer*4 :: mNhB; ! ×èñëî óçëîâ ïî øèðèíå
integer*4 :: mLDA; ! Øèðèíà ëåíòû ìàòðèöû æåñòêîñòè
real*8 :: mC; ! Òåïëîåìêîñòü
real*8 :: mK; ! Êîýôôèöèåíò òåïëîïðîâîäíîñòè
real*8 :: mR; ! ïëîòíîñòü
type(gr2d) :: mGr; ! Ñåòêà ÌÊÝ
type(ELEM) :: mEL4; ! Äàííûå èñïîëüçóåìîãî ÊÝ
integer*4 :: mContrPoints(9); ! Êîíòðîëüíûå òî÷êè
REAL*8 :: mcpX(9),mcpY(9);
real*8, dimension(4,4) :: est; ! Ìàòðèöà æåñòêîñòè òåêóùåãî ýëåìåíòà
real*8, dimension(4) :: r; ! Âåêòîð íàãðóçîê òåêóùåãî ÊÝ
real*8, dimension(:), allocatable :: mB; ! Ãëîáàëüíûé âåêòîð íàãðóçîê
real*8, dimension(:,:), allocatable :: mA; ! Ãëîáàëüíàÿ ìàòðèöà æåñòêîñòè
real*8, dimension(:), allocatable :: mX; ! Ãëîáàëüíûé âåêòîð íåèçâåñòíûõ
end module GlobData;
Rys. 9.4. Ãëîáàëüíûå ïåðåìåííûå
1.3.3. Ïðîöåäóðà èíèöèàëèçàöèè äàííûõ äëÿ ýëåìåíòà
Çàïîëíåíèå ïåðåìåííîé mEL4 òèïà ELEM ïðîèñõîäèò â ïðîãðàììå IniEL4 ïîêàçàííîé íà rys. 9.5.  ïåðâîé ÷àñòè ïðîãðàììû çàäàåòñÿ ÷èñëî óçëîâ ýëåìåíòà nbn ðàâíîå 4. Âûáðàííàÿ ñõåìà èíòåãðèðîâàíèÿ ñîäåðæèò 4 òî÷êè, ïîýòîìó çàäàíî N_p=4. Âûäåëÿåòñÿ ïàìÿòü ìàññèâàì N1, N2, Nf, P, W è L. Äàëåå çàïîëíÿþòñÿ ëîêàëüíûå êîîðäèíàòû òî÷åê èíòåãðèðîâàíèÿ Ãàóññà.  êàæäîé òî÷êå èíòåãðèðîâàíèÿ ðàñ÷èòûâàþòñÿ âåëè÷èíû N1, N2 i Nf.
Àíàëîãè÷íî çàïîëíÿþòñÿ ñòðóêòóðû äëÿ ïîâåðõíîñòåé ýëåìåíòîâ. Ïðè èñïîëüçîâàíèè èíîãî òèïà ÊÝ íåîáõîäèìî ñîîòâåòñòâóþùèì îáðàçîì èçìåíèòü ïåðåìåííóþ mEL4 è ñåòêó ÊÝ.
SUBROUTINE IniEL4();
! ******************************************************************
! Inicjalizacja danych dla elementu
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ******************************************************************
use GlobData;
implicit none;
REAL*8 :: Alfa,L1,L2,e,n,SN;
integer*4 :: iP,i;
Alfa = 1/(sqrt(3.0));
mEL4%nbn=4;
mEL4%N_p=4;
ALLOCATE(mEL4%N1( mEL4%nbn, mEL4%N_p ));
ALLOCATE(mEL4%N2( mEL4%nbn, mEL4%N_p ));
ALLOCATE(mEL4%Nf( mEL4%nbn, mEL4%N_p ));
ALLOCATE(mEL4%P ( mEL4%N_p ));
ALLOCATE(mEL4%W ( mEL4%N_p ));
ALLOCATE(mEL4%L ( mEL4%nbn ));
mEL4%W = 1.0;
mEL4%P(1)%E = -Alfa;
mEL4%P(1)%N = -Alfa;
mEL4%P(2)%E = Alfa;
mEL4%P(2)%N = -Alfa;
mEL4%P(3)%E = Alfa;
mEL4%P(3)%N = Alfa;
mEL4%P(4)%E = -Alfa;
mEL4%P(4)%N = Alfa;
mEL4%L(1)%E = -1;
mEL4%L(1)%N = -1;
mEL4%L(2)%E = 1;
mEL4%L(2)%N = -1;
mEL4%L(3)%E = 1;
mEL4%L(3)%N = 1;
mEL4%L(4)%E = -1;
mEL4%L(4)%N = 1;
DO iP = 1, mEL4%N_p
L1 = mEL4%P(iP)%E
L2 = mEL4%P(iP)%N
mEL4%Nf(1,iP) = 0.25*(1-L1)*(1-L2);
mEL4%N1(1,iP) = -0.25*(1-L2);
mEL4%N2(1,iP) = -0.25*(1-L1);
mEL4%Nf(2,iP) = 0.25*(1+L1)*(1-L2);
mEL4%N1(2,iP) = 0.25*(1-L2);
mEL4%N2(2,iP) = -0.25*(1+L1);
mEL4%Nf(3,iP) = 0.25*(1+L1)*(1+L2);
mEL4%N1(3,iP) = 0.25*(1+L2);
mEL4%N2(3,iP) = 0.25*(1+L1);
mEL4%Nf(4,iP) = 0.25*(1-L1)*(1+L2);
mEL4%N1(4,iP) = -0.25*(1+L2);
mEL4%N2(4,iP) = 0.25*(1-L1);
end do;
DO i = 1, 4
mEL4%Sf(i)%N_p = 2;
ALLOCATE( mEL4%Sf(i)%P( mEL4%Sf(i)%N_p ) );
ALLOCATE( mEL4%Sf(i)%W( mEL4%Sf(i)%N_p ) );
mEL4%Sf(i)%W = 1;
ALLOCATE( mEL4%Sf(i)%N1( mEL4%nbn, mEL4%Sf(i)%N_p));
ALLOCATE( mEL4%Sf(i)%N2( mEL4%nbn, mEL4%Sf(i)%N_p));
ALLOCATE( mEL4%Sf(i)%Nf( mEL4%nbn, mEL4%Sf(i)%N_p));
ALLOCATE( mEL4%Sf(i)%UZEL( 2 ) );
END DO;
mEL4%Sf(1)%UZEL(1) = 4; mEL4%Sf(1)%UZEL(2) = 1;
mEL4%Sf(2)%UZEL(1) = 1; mEL4%Sf(2)%UZEL(2) = 2;
mEL4%Sf(3)%UZEL(1) = 2; mEL4%Sf(3)%UZEL(2) = 3;
mEL4%Sf(4)%UZEL(1) = 3; mEL4%Sf(4)%UZEL(2) = 4;
mEL4%Sf(1)%P(1)%N=Alfa; mEL4%Sf(1)%P(1)%E=-1;
mEL4%Sf(1)%P(2)%N=-Alfa; mEL4%Sf(1)%P(2)%E=-1;
mEL4%Sf(2)%P(1)%N=-1; mEL4%Sf(2)%P(1)%E=-Alfa;
mEL4%Sf(2)%P(2)%N=-1; mEL4%Sf(2)%P(2)%E=Alfa;
mEL4%Sf(3)%P(1)%N=-Alfa; mEL4%Sf(3)%P(1)%E=1;
mEL4%Sf(3)%P(2)%N=Alfa; mEL4%Sf(3)%P(2)%E=1;
mEL4%Sf(4)%P(1)%N=1; mEL4%Sf(4)%P(1)%E=Alfa;
mEL4%Sf(4)%P(2)%N=1; mEL4%Sf(4)%P(2)%E=-Alfa;
do i = 1, 4
DO iP = 1, 2
e = mEL4%Sf(i)%P(ip)%E;
n = mEL4%Sf(i)%P(ip)%N;
mEL4%Sf(i)%Nf(1,ip) = 0.25*(1-e)*(1-n);
mEL4%Sf(i)%Nf(2,ip) = 0.25*(1+e)*(1-n);
mEL4%Sf(i)%Nf(3,ip) = 0.25*(1+e)*(1+n);
mEL4%Sf(i)%Nf(4,ip) = 0.25*(1-e)*(1+n);
mEL4%Sf(i)%N1(1,ip) = -0.25*(1-n);
mEL4%Sf(i)%N1(2,ip) = 0.25*(1-n);
mEL4%Sf(i)%N1(3,ip) = 0.25*(1+n);
mEL4%Sf(i)%N1(4,ip) = -0.25*(1+n);
mEL4%Sf(i)%N2(1,ip) = -0.25*(1-e);
mEL4%Sf(i)%N2(2,ip) = -0.25*(1+e);
mEL4%Sf(i)%N2(3,ip) = 0.25*(1+e);
mEL4%Sf(i)%N2(4,ip) = 0.25*(1-e);
END DO;
end do;
END SUBROUTINE IniEL4; (à)
(á)
Rys. 9.5. Èíèöèàëèçàöèÿ äàííûõ äëÿ ÊÝ (à) è èñïîëüçóåìûé ÊÝ (á)
1.3.4. Ïðîöåäóðà ñ÷èòûâàíèÿ èñõîäíûõ äàííûõ è âûäåëåíèÿ ïàìÿòè ìàññèâàì
Òåêñò ïðîãðàììû íà rys. 9.6. Ïðîãðàììà îòêðûâàåò òåêñòîâûé ôàéë indata.t2d, ñòðóêòóðà êîòîðîãî îïèñàíà âûøå.
Subroutine InpData();
! ******************************************************************
! Wczytywanie danych
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ******************************************************************
use GlobData;
implicit none;
integer nFile,i,j,iost;
nFile = 89
OPEN (nFile, FILE='indata.t2d', action='read')
read(nFile, *);
read(nFile, *);
read(nFile, *);
read(nFile, *);
read(nFile, *) mTbegin;
read(nFile, *) mTime;
read(nFile, *) mdTime;
read(nFile, *) mT_otoczenia;
read(nFile, *) mAlfa;
read(nFile, *) mH0;
read(nFile, *) mB0;
read(nFile, *) mNhH;
read(nFile, *) mNhB;
read(nFile, *) mC;
read(nFile, *) mK;
read(nFile, *) mR;
close (nFile);
end Subroutine InpData
Rys. 9.6. Ïðîãðàììà ñ÷èòûâàíèÿ äàííûõ.
Âûäåëåíèå ïàìÿòè äëÿ ìàòðèöû æåñòêîñòè íà÷èíàåòñÿ ðàñ÷åòîì øèðèíû ïîëîñû íåíóëåâûõ ýëåìåíòîâ mLDA, çàòåì âûäåëåòñÿ ïàìÿòü ìàòðèöå æåñòêîñòè mA, âåêòîðó íåèçâåñòíûõ mX è âåòîðó íàãðóçîê mB.  ìåòîäè÷åñêèõ öåëÿõ ïîêàçàí òàêæå âàðèàíò âûäåëåíèÿ ïàìÿòè äëÿ ïîëíîé ìàòðèöû æåñòêîñòè ðàçìåðîì nh*nh.
SUBROUTINE ALLOCATE_Matrix;
use GlobData;
implicit none;
integer*4 :: NEL, i, j, ii, jj, jB, nk(4),NeMaxB, ierr;
mLDA=0;
do NEL = 1,mGr%ne
do i=1,mEL4%nbn
nk(i) = mGr%EL(NEL)%nop(i);
end do;
do i=1,mEL4%nbn
! * ii - wiersz *
ii=nk(i);
do j=1,mEL4%nbn
! * jj - kolumna *
jj=nk(j);
jB = jj-ii+1;
if(jB>=mLDA) then
mLDA=jB;
NeMaxB = NEL;
end if;
end do;
end do;
end do;
ALLOCATE ( mA(mLDA, mGr%nh), STAT=iErr );
! ALLOCATE ( mA(mGr%nh, mGr%nh), STAT=iErr );
! Pelna matryca MES
ALLOCATE ( mB(mGr%nh), STAT=iErr );
ALLOCATE ( mX(mGr%nh), STAT=iErr );
END SUBROUTINE ALLOCATE_Matrix
Rys. 9.7. Ïðîãðàììà ðàñ÷åòà è âûäåëåíèÿ ïàìÿòè ìàññèâàì.
1.3.5. Ïðîöåäóðà ðåøåíèÿ çàäà÷è íà òåêóùåì øàãå ïî âðåìåíè
Òåêñò ïðîãðàììû ïðåäñòàâëåí íà rys. 9.8.
SUBROUTINE SOLVER( );
! ******************************************************************
! Solver
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ******************************************************************
use GlobData;
use msimsl;
implicit none;
integer*4 :: NEL;
integer*4, dimension(4) :: nk;
integer*4 :: jB,iB,NeMaxB,NCODA,i,j,ii,jj;
integer*4 :: iErr;
mA=0; mB=0; mX=0;
do NEL=1, mGr%ne
do i=1, mEL4%nbn
nk(i) = mGr%EL(NEL)%nop(i);
end do;
CALL FeSM_heat( NEL );
do i=1, mEL4%nbn
ii = nk(i); ! wiersz w macierzy pełnej
do j=1, mEL4%nbn
jj = nk(j); ! Kolumna w macierzy pełnej
iB = mLDA + ii - jj;! Wiersz w macierzy pasmowej
if ((jj>=ii).AND.(iB<=mLDA)) then
mA(iB, jj) = mA(iB, jj) + est(i,j); ! Wypełnienie macierzy pasmowej
end if;
! mA(ii,jj) = mA(ii,jj) + est(i,j); ! Wypełnienie macierzy pełnej
end do;
mB(ii) = mB(ii) + r(i);
end do;
end do;
NCODA = mLDA-1;
CALL DLSAQS (mGr%nh, mA, mLDA, NCODA, mB, mX)
! Program do rozwiązania układu równań z macierzą pasmową
! NCODA=1;
! CALL DLSARG (mGr%nh, mA, mGr%nh, mB, NCODA, mX)
! Program do rozwiązania układu równań z macierzą pełną
do i=1, mGr%nh
mGr%nd(i)%CR = ( mGr%nd(i)%t-mX(i) ) / mdTime;
mGr%nd(i)%t = mX(i);
end do;
END SUBROUTINE SOLVER;
Rys. 9.8. Ïðîãðàììà çàïîëíåíèÿ ãëîáàëüíîé ìàòðèöû æåñòêîñòè è ðåøåíèÿ ñèñòåìû óðàâíåíèé
 äàííîé ïðîãðàììå èñïîëüçîâàíû ñëåäóþùèå îáîçíà÷åíèÿ:
NEL – òåêóùèé íîìåð êîíå÷íîãî ýëåìåíòà;
Nk – ìàññèâ, ñîäåðæàùèé íîìåðà óçëîâ, îáðàçóþùèé äàííûé ýëåìåíò;
FeSM_heat – ïðîãðàììà ðàñ÷åòà ìàòðèöû æåñòêîñòè è âåêòîðà íàãðóçîê òåêóùåãî ýëåìåíòà NEL;
DLSAQS – ïðîãðàììà ðåøåíèÿ ñèñòåáó óðàâíåíèé ñ ëåíòî÷íîé ìàòðèöåé;
DLSARG – ïðîãðàììà ðåøåíèÿ ñèñòåìû óðàâíåíèé ñ ïîëíîé ìàòðèöåé.
1.3.6. Ïðîãðàììà ðàñ÷åòà ìàòðèöû æåñòêîñòè è âåêòîðà íàãðóçîê òåêóùåãî ýëåìåíòà
Ïðîãðàììà FeSM_heat (rys. 9.9) ñîñòîèò èç äâóõ ÷àñòåé:
PRE_heat_mat – ðàñ÷åò ìàòðèöû è âåêòîðà íàãðóçîê;
PRE_heat_pov_mat – ó÷åò ãðàíè÷íûõ óñëîâèé.
SUBROUTINE FeSM_heat( NEL );
! ******************************************************************
! Generacja macierzy elementu NEL
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ******************************************************************
use GlobData;
implicit none;
integer*4 :: NEL; ! Numer bieżącego Elementu
est = 0;
r = 0;
call PRE_heat_mat( NEL );
call PRE_heat_pov_mat(NEL );
END SUBROUTINE FeSM_heat;
Rys. 9.9. Ðàñ÷åò ìàòðèöû æåñòêîñòè òåêóùåãî ÊÝ (NEL)
Òåêñò ïðîãðàììû PRE_heat_mat ïîêàçàíî íà rys. 9.10. Ïðîãðàììà Jacob_2d ðàñ÷èòûâàåò ìàòðèöó ßêîáè, îáðàòíóþ åé ìàòðèöó è îïðåäåëèòåëü ìàòðèöû ÿêîáè â êàæäîé òî÷êå Ãàóññà.
subroutine PRE_heat_mat( NEL );
use GlobData;
implicit none;
integer*4 :: NEL;
INTEGER*4 :: I,j,N,P,Id;
REAL*8, DIMENSION(2,2) :: J_,J_inv; ! Macierz Jakobiego i macierz odwrotna
REAL*8 :: DetJ,Ni,Nn,Hin,Cin;
real*8 :: T0p;
REAL*8, DIMENSION(4) :: Ndx,Ndy; ! Pochodne funkcji kształtu względem X i Y
REAL*8, DIMENSION(4) :: X,Y,Temp_0;
DO I=1,mEL4%NBN
Id = ABS( mGr%EL(NEL)%NOP(I) );
X(I) = mGr%nd(Id)%x /1000.0; ! [m]
Y(I) = mGr%nd(Id)%y /1000.0; ! [m]
Temp_0(i)= mGr%nd(Id)%t;
END DO
DO P=1, mEL4%N_p
CALL Jacob_2d(J_, J_inv, P, mEL4%N_p, mEL4%NBN, mEL4%N1,mEL4%N2, X, Y, DetJ);
T0p=0;
DO i=1, mEL4%NBN
Ndx(i)=mEL4%N1(i,p)*J_inv(1,1) + mEL4%N2(i,p)*J_inv(1,2) ! [m-1]
Ndy(i)=mEL4%N1(i,p)*J_inv(2,1) + mEL4%N2(i,p)*J_inv(2,2) ! [m-1]
Ni = mEL4%Nf(i,p);
T0p = T0p + Temp_0(i)*Ni;
END DO;
detJ = abs(detJ)*mEL4%W(p); ! [m2]
DO N=1,mEL4%NBN
DO I=1,mEL4%NBN
Ni = mEL4%Nf(i,p);
Nn = mEL4%Nf(n,p);
Hin = mK*(Ndx(n)*Ndx(i) + Ndy(n)*Ndy(i))*DetJ; ! Wzór (4.14)
Cin = mC*mR*Nn*Ni*detJ; ! Wzór (5.3)
est(n,i) = est(n,i) + Hin + Cin/mdTime; ! Wzór (5.8)
r(n) = r(n) + (Cin/mdTime)*T0p;
END DO
! Q = 0.8*Tp*Hp
! r(n)=r(n) + Nn*Q*DetJ;
END DO
END DO
end subroutine PRE_heat_mat;
Rys. 9.10. Ïðîãðàììà ðàñ÷åòà êîìïîíåíòîâ ìàòðèöû æåñòêîñòè è âåêòîðà íàãðóçîê ÊÝ NEL.
1.3.7. Ïðîãðàììà ó÷åòà ãðàíè÷íûõ óñëîâèé
Ïðîãðàììà ó÷åòà ãðàíè÷íûõ óñëîâèé â ìàòðèöå æåñòêîñòè è âåêòîðå íàãðóçîê ïîêàçàíà íà rys. 9.11. Äëÿ òåêóùåãî ÊÝ âûïîëíÿåòñÿ öèêë ïî âñåì êîíòàêíòûì ïîâåðõíîñòÿì (èõ ÷èñëî ðàâíî mGr%EL(NEL)%nPov). Çàòåì âûïîëíÿåòñÿ èíòåãðèðîâàíèå ïî êàæäîé òàêîé ïîâåðõíîñòè.
subroutine PRE_heat_pov_mat( NEL );
! ****************************************************
! WARUNKI BRZEGOWE
! Milenin Andre, AGH, 2008, milenin@agh.edu.pl
! ****************************************************
use GlobData;
implicit none;
integer*4 :: NEL;
INTEGER*4 :: I,j,N,P,Id,iPov;
REAL*8 :: DetJ,Ni,Nn,Pn;
REAL*8, DIMENSION(4) :: X,Y;
DO I=1,4
Id = ABS(mGr%EL(NEL)%NOP(I));
X(I) = mGr%nd(Id)%x /1000.0; ! [m]
Y(I) = mGr%nd(Id)%y /1000.0; ! [m]
END DO
DO iPov=1, mGr%EL(NEL)%nPov
id = mGr%EL(NEL)%aPov(iPov); ! Lokalny numer powierzchni elementu
SELECT CASE (id)
case(1)
DetJ=sqrt( (X(4)-X(1))**2 + (Y(4)-Y(1))**2 );
case(2)
DetJ=sqrt( (X(1)-X(2))**2 + (Y(1)-Y(2))**2 );
case(3)
DetJ=sqrt( (X(2)-X(3))**2 + (Y(2)-Y(3))**2 );
case(4)
DetJ=sqrt( (X(3)-X(4))**2 + (Y(3)-Y(4))**2 );
END SELECT
DO P=1, 2
DO n=1, 4
! Pętla po grupach równań równych wierszom macierzy
DO i=1, 4
! Pętla po kolumnach macierzy
Ni = mEL4%Sf(id)%Nf(i,p);
Nn = mEL4%Sf(id)%Nf(n,p);
est(n,i) = est(n,i) + mAlfa*Nn*Ni*DetJ; ! Wzór (4.14)
END DO
Pn = mAlfa*mT_otoczenia*Nn*DetJ; ! Wzór (4.15)
R(n) = R(n) + Pn;
END DO
END DO ! P
END DO; ! iPov
end subroutine PRE_heat_pov_mat;
Rys. 9.11. Ïðîãðàììà ó÷åòà ãðàíè÷íûõ óñëîâèé äëÿ òåêóùåãî ÊÝ.
1.4. Ïðèìåð ðàñ÷åòà
Ðåçóëüòàòû ðàñ÷åòà äëÿ äàííûéõ, ïðèâåäåííûõ íà rys. 9.1 ïîêàçàíû íà rys. 9.12 äëÿ òî÷êè â öåíòðå ñå÷åíèÿ (êðèâàÿ 1) è íà åãî ïîâåðõíîñòè (êðèâàÿ 2).
Rys. 9.12. Ïðèìåð ðåçóëüòàòîâ ðàñ÷åòà äëÿ òî÷êè â öåíòðå ñå÷åíèÿ (êðèâàÿ 1) è íà åãî ïîâåðõíîñòè (êðèâàÿ 2).