DISTRIBUTIA TEMPERATURII INTR-O PLACA METALICA DREPTUNGHIULARA FARA SURSA DE CALDURA

DISTRIBUTIA TEMPERATURII INTR-O PLACA METALICA DREPTUNGHIULARA FARA SURSA DE CALDURA

Se considera o placa dreptungiulara metalica de dimensiuni 0,5X0,4[m] avand grosimea de 0,01m.Conductivitatea termica a materialului placii este λ=1000[W/m/K.Frontiera de V a placii primeste un flux q=500W/m² iar frontierele S si E sunt mentinute la temperatura constanta de 100C

Sa se calculeze distributia stationara de temperatura din placa folosind reteaua de dicretizare din figura de mai jos.

Rezolvare :

Fig.1

Reteaua si conditiile la limita pentru o problema de conductie termica 2D stationara

Ecuatia conductiei termice 2D stationara pentru conditiile enuntate este:

Ecuatia discretizata pentru un nod interior este:

unde:

Sistemul de ecuatii care rezulta este:

40T₈=10T₂+10 T₁₄+10T₇+10T₉

40T₉=10T₃+10 T₁₅+10T₈+10T₁₀

40T₁₀=10T₄+10 T₁₆+10T₉+10T₁₁

40T₁₁=10T₅+10 T₁₇+10T₁₀+10T₁₂

40T₁₄=10T₈+10 T₂₀+10T₁₅

40T₁₅=10T₉+10 T₂₁+10T₁₄+10T₁₆

40T₁₆=10T₁₀+10 T₂₂+10T₁₅+10T₁₇

40T₁₇=10T₁₁+10 T₂₃+10T₁₆+10T₁₈

40T₂₀=10T₁₄+10 T₂₆+10T₁₉+10T₂₁

40T₂₁=10T₁₅+10 T₂₇+10T₂₀+10T₂₂

40T₂₂=10T₁₆+10 T₂₈+10T₂₁+10T₂₃

40T₂₃=10T₁₇+10 T₂₉+10T₂₂+10T₂₄

Pentru nodurile de pe frontiera V (2, 3, 4, 5) ecuatia discretizata se obtine integrand ecuatia pe jumatatea de volum hasurata din figura :

Ecuatia discretizata pentru un nod pe frontiera Vest este:

b=qAp=500000(0.1·0,01)=500

Sistemul de ecuatii care rezulta este :

20T₂=10T₈+5 T₁+5T₃+500

20T₃=10T₉+5 T₂+5T₄+500

20T₄=10T₁₀+5 T₃+5T₅+500

20T₅=10T₁₁+5 T₄+5T₆+500

Pentru nodurile de pe frontiera E (26, 27, 28, 29) ecuatia discretizata se obtine integrand ecuatia pe jumatatea de volum hasurata din figura :

Ecuatia pentru un nod pe frontiera Est este:

Sistemul de ecuatii care rezulta este :

20T₂₆=10T₂₀+5 T₂₅+5T₂₇

20T₂₇=10T₂₁+5 T₂₆+5T₂₈

20T₂₈=10T₂₂+5 T₂₇+5T₂₉

20T₂₉=10T₂₃+5 T₂₈+5T₃₀

Pentru nodurile de pe frontiera S (7, 13, 19) ecuatia discretizata se obtine integrand ecuatia pe jumatatea de volum hasurata din figura :

Ecuatia pentru nodul situat pe frontiera Sud este:

Sistemul de ecuatii care rezulta este :

20T₇=10T₁+5 T₁₃+5T₈

20T₁₃=10T₇+5 T₁₉+5T₁₄

20T₁₉=10T₁₃+5 T₂₅+5T₂₀

Ne raman de tratat nodurile de colt.

Ecuatia discretizata pentru nodul (1) se obtine integrand ecuatia pe sfertul de volum hasurat din figura :

Ecuatia pentru nodul "1" este:

b=Apq=250

Ecuatia ce trebuie rezolvata pentru nodul (1) devine :

10T₁=5T₇+5 T₂+250

Ecuatia discretizata pentru nodul (25) se obtine integrand ecuatia pe sfertul de volum hasurat din figura :

Ecuatia pentru nodul "25" este:

Ecuatia devine :

10T₂₅=5T₁₉+5 T₂₆

Sistemul de ecuatii este urmatorul:

Programul in Fortran :

c------------ecuatia rezolvata:(divk(gradt))+s=0-------- ----- ------ ---

c------------variabilele de intrare-------- ----- ------ ----- ----- ----------

c nnx : numarul de noduri pe directia x

c nny : numarul de noduri pe directia y

c l : lungimea domeniului de calcul

c h : inaltimea domeniului de calcul

c tb : conditia Dirichlet

c------------conditiile la limita-------- ----- ------ ----- ----- ----- ----- ------

c Dirichlet pe toate frontierele(W,E,S,N)

c t(0,y)=0,t(w,y)=0,t(x,0)=0,t(x,h)=0

c------------variabilele de intrare-------- ----- ------ ----- ----- -----------

c temp : temperatura t(x,y)

c-------- ----- ------ -------- ----- ------ -------- ----- ------ ------

parameter(nnx=21, nny=26, nit=1500)

double precision temp(nnx,nny),ax(nnx,nny),ay(nnx,nny)

double precision a(nnx,nny),b(nnx,nny),c(nnx,nny),wk(nnx,nny)

double precision dif(nnx,nny),temp2(nnx,nny)

double precision a1(nny),b1(nny),c1(nny),wk1(nny),temp1(nny)

double precision dx,dy,tb,l,h,la,gz,s,err

data tb/0.0/,l/0.04/,h/0.05/,la/4.0/,gz/1.0/,err/0.001/

datas/40000000.0/

c------------calculul pasului spatiului(dx)-------- ----- ------ ---

dx=l/(nnx-1)

dy=h/(nny-1)

write(*,*)'dx=',dx

write(*,*)'dy=',dy

c------------initializarea temperaturii-------- ----- ------ -------

do i=1,nnx

do j=1,nny

temp(i,j)=0.01

enddo

enddo

c-----------initializarea conditiilor la limita de tip Dirichlet----- ----- -----

c----- ----- -------------la frontiera NORD-------- ----- ------ -----

do i=1,nnx

temp(i,nny)=tb

enddo

c----- ----- -------------la frontiera SUD-------- ----- ------ -----

do i=1,nnx

temp(i,1)=tb

enddo

c----- ----- -------------la frontiera WEST-------- ----- ------ ----

do j=1,nny

temp(1,j)=tb

enddo

c----- ----- -------------la frontiera EST-------- ----- ------ -----

do j=1,nny

temp(nnx,j)=tb

enddo

c----- ----- ------------bucla de iteratie-------- ----- ------ ------

do k=1,nit

c------------formarea vectorilor diagonalelor-------- ----- ------ -

do i=2,nnx-2

if(i.eq.2)then

c------------diagonala inferioara-------- ----- ------ -------------

do j=2,nny-1

if(j.eq.2)then

a(i,j)=0.0

a1(j-1)=a(i,j)

else

a(i,j)=-la*(dx*gz)/dy

a1(j-1)=a(i,j)

endif

enddo

c------------diagonala superioara-------- ----- ------ ----- ----- -----

do j=2,nny-2

if(j.eq.2)then

c(i,j)=-la*(dx*gz)/dy

c1(j-1)=c(i,j)

else

c(i,j)=-la*(dx*gz)/dy

c1(j-1)=c(i,j)

endif

enddo

c(i,nny-1)=0.0

c1(nny-2)=c(i,nny-1)

c------------diagonala principala-------- ----- ------ ----- ----- -----

do j=2,nny-1

if(j.eq.2)then

b(i,j)=2.0*la*(dy*gz)/dx+2.0*la*(dx*gz)/dy

b1(j-1)=b(i,j)

else

b(i,j)=2.0*la*(dy*gz)/dx+2.0*la*(dx*gz)/dy

b1(j-1)=b(i,j)

endif

enddo

c------------formarea initiala a termenului liber wk----- ----- -----------

do j=2,nny-2

if(j.eq.2)then

wk(i,j)=(la*(dx*gz)/dy)*temp(i,j-1)+

(la*(dy*gz)/dx)*temp(i-1,j)+(la*(dy*gz)/dx)*temp(i+1,j)+ s*dx*dy*1.0

wk1(j-1)=wk(i,j)

else

wk(i,j)=(la*(dy*gz)/dx)*temp(i-1,j)+

* (la*(dy*gz)/dx)*temp(i+1,j)+s*dx*dy*1.0

wk1(j-1)=wk(i,j)

endif

enddo

wk(i,nny-1)=(la*(dy*gz)/dx)*temp(i-1,nny-1)+

* (la*(dy*gz)/dx)*temp(i+1,nny-1)+s*dx*dy*1.0+

* (la*(dx*gz)/dy)*temp(i,nny)

wk1(nny-2)=wk(i,nny-1)

c-----------rezolvarea sistemului-------- ----- ------ ----- ----- -----

CALL TRIDAG(a1,b1,c1,wk1,temp1,nny-2)

c-------formarea solutiei pe toata linia verticala(linia aplicarii TDMA)----- ----- -----

do j=2,nny-1

temp(i,j)=temp1(j-1)

enddo

else

c------------diagonala inferioara-------- ----- ------ -------------

do j=2,nny-1

if(j.eq.2)then

a(i,j)=0.0

a1(j-1)=a(i,j)

else

a(i,j)=-la*(dx*gz)/dy

a1(j-1)=a(i,j)

endif

enddo

c------------diagonala superioara-------- ----- ------ ----- ----- -----

do j=2,nny-2

if(j.eq.2)then

c(i,j)=-la*(dx*gz)/dy

c1(j-1)=c(i,j)

else

c(i,j)=-la*(dx*gz)/dy

c1(j-1)=c(i,j)

endif

enddo

c(i,nny-1)=0.0

c1(nny-2)=c(i,nny-1)

c------------diagonala principala-------- ----- ------ ----- ----- -----

do j=2,nny-1

if(j.eq.2)then

b(i,j)=2.0*la*(dy*gz)/dx+2.0*la*(dx*gz)/dy

b1(j-1)=b(i,j)

else

b(i,j)=2.0*la*(dy*gz)/dx+2.0*la*(dx*gz)/dy

b1(j-1)=b(i,j)

endif

enddo

c------------formarea initiala a termenului liber wk----- ----- --------- ----- -----

do j=2,nny-2

if(j.eq.2)then

wk(i,j)=(la*(dx*gz)/dy)*temp(i,j-1)+

(la*(dy*gz)/dx)*temp(i-1,j)+(la*(dy*gz)/dx)*temp(i+1,j)+

* s*dx*dy*1.0

wk1(j-1)=wk(i,j)

else

wk(i,j)=(la*(dy*gz)/dx)*temp(i-1,j)+

* (la*(dy*gz)/dx)*temp(i+1,j)+(s*dx*dy*1.0)

wk1(j-1)=wk(i,j)

endif

enddo

wk(i,nny-1)=(la*(dy*gz)/dx)*temp(i-1,nny-1)+

* (la*(dy*gz)/dx)*temp(i+1,nny-1)+(la*(dx*gz)/dy)*temp(i,nny)+

* s*dx*dy*1.0

wk1(nny-2)=wk(i,nny-1)

c-----------rezolvarea sisitemului-------- ----- ------ ----- ----- -----

CALL TRIDAG(a1,b1,c1,wk1,temp1,nny-2)

c----- ----- ----formarea solutiei pe toata linia verticala----- ----- ----- ----- -----

do j=2,nny-1

temp(i,j)=temp1(j-1)

enddo

endif

enddo

c------------diagonala inferioara-------- ----- ------ ----- ----- -----

do j=2,nny-1

if(j.eq.2)then

a(nnx-1,j)=0.0

a1(j-1)=a(nnx-1,j)

else

a(nnx-1,j)=-la*(dx*gz)/dy

a1(j-1)=a(nnx-1,j)

endif

enddo

c------------diagonala superioara-------- ----- ------ ----- ----- -----

do j=2,nny-2

if(j.eq.2)then

c(nnx-1,j)=-la*(dx*gz)/dy

c1(j-1)=c(nnx-1,j)

else

c(nnx-1,j)=-la*(dx*gz)/dy

c1(j-1)=c(nnx-1,j)

endif

enddo

c(nnx-1,nny-1)=0.0

c1(nny-2)=c(nnx-1,nny-1)

c------------diagonala principala-------- ----- ------ ----- ----- -----

do j=2,nny-1

if(j.eq.2)then

b(nnx-1,j)=2.0*la*(dy*gz)/dx+2.0*la*(dx*gz)/dy

b1(j-1)=b(nnx-1,j)

else

b(nnx-1,j)=2.0*la*(dy*gz)/dx+2.0*la*(dx*gz)/dy

b1(j-1)=b(nnx-1,j)

endif

enddo

c------------formarea initiala a termenului liber wk----- ----- --------- ----- -----

do j=2,nny-2

if(j.eq.2)then

wk(nnx-1,j)=(la*(dx*gz)/dy)*temp(nnx-1,j-1)+

* (la*(dy*gz)/dx)*temp(nnx-2,j)+(la*(dy*gz)/dx)*temp(nnx,j)+

* s*dx*dy*1.0

wk1(j-1)=wk(nnx-1,j)

else

wk(nnx-1,j)=(la*(dy*gz)/dx)*temp(nnx-2,j)+

* (la*(dy*gz)/dx)*temp(nnx,j)+s*dx*dy*1.0

wk1(j-1)=wk(nnx-1,j)

endif

enddo

wk(nnx-1,nny-1)=(la*(dy*gz)/dx)*temp(nnx-2,nny-1)+

* (la*(dy*gz)/dx)*temp(nnx,nny-1)+

* (la*(dx*gz)/dy)*temp(nnx-1,nny)+s*dx*dy*1.0

wk1(nny-2)=wk(nnx-1,nny-1)

c-----------rezolvarea sistemului-------- ----- ------ ----- ----- -----

CALL TRIDAG(a1,b1,c1,wk1,temp1,nny-2)

c----- ----- ----formarea solutiei pe toata linia verticala----- ----- ----- ----- -----

do j=2,nny-1

temp(i,j)=temp1(j-1)

enddo

c-----------verificarea criteriului de convergenta----- ----- --------- ----- --------

do i=2,nnx-1

do j=2,nny-1

dif(i,j)=abs(100.0*(temp(i,j)-temp2(i,j))/temp(i,j))

enddo

enddo

difmax=dif(2,2)

do i=2,nnx-1

do j=2,nny-1

if(dif(i,j).gt.difmax)then

difmax=dif(i,j)

else

endif

enddo

enddo

write(*,*)'iter,difmax,err',k,difmax,err

write(*,*)'iter,difmax,err',difmax,err

if(difmax.le.err)then

write(*,*)'sortarea dupa criteriul de convergenta,k=',k

goto 100

else

endif

c------------initializarea temperaturii pasului precedent ----- ----- ----- ----- ----

do i=2,nnx-1

do j=2,nny-1

temp2(i,j)=temp(i,j)

enddo

enddo

c-----------sfarsitul buclei de iteratie-------- ----- ------ --------

enddo

100 continue

c-----------scrierea solutiei-------- ----- ------ ----- ----- ----------

open(20,file='laborator8.prn')

do i=1,nnx

write(20,101)(temp(i,j),j=1,nny)

enddo

close(20)

101 format(26(1x,f7.2))

stop

end

subroutine tridag(a,b,c,r,u,n)

parameter (nmax=200000)

integer j

double precision a(n),b(n),c(n),r(n),u(n)

double precision bet, gam(nmax)

c----------a=diagonala inferioara-------- ----- ------ ----- ----- ----

c----------b=diagonala principala-------- ----- ------ ----- ----- ----

c----------c=diagonala superioara-------- ----- ------ ----- ----- ----

c----------r=termenul din partea dreapta-------- ----- ------ -------

c----------u=termenul necunoscut-------- ----- ------ ----- ----- -----

c----------n=dimensiunea sistemului-------- ----- ------ ------------

c----------nmax=numarul maxim al ecuatiei-------- ----- ------ ------

if (b(1).eq.0) pause 'tridag:rewrite equations !'

bet=b(1)

u(1)=r(1)/bet

do j=2,n

gam(j)=c(j-1)/bet

bet=b(j)-a(j)*gam(j)

if(bet.eq.0) pause 'tridag failed'

u(j)=(r(j)-a(j)*u(j-1))/bet

enddo

do j=n-1,1,-1

u(j)=u(j)-gam(j+1)*u(j+1)

enddo

return

end

Distributia temperaturii(in Quickfield) :