El modelo anterior de conductividad térmica efectiva uniforme no es bueno porque en el contacto con la FPGA el camino térmico se ensancha, lo que puede contabilizarse también manteniendo el espesor del camino fijo, Lz=LzPCB, y aumentando la k efectiva. El contacto térmico de la FPGA con la PCB es muy bueno porque suele ser con una matriz de bolitas de soldadura (BGA) de un milímetro de separación, i.e. de unos 50·50=2500 contactos para un chip de 50·50 mm².

Si la FPGA cubriese toda la placa, la conductividad efectiva sería:

>	eqkeff_underFPGA:=k=Sum(k[i]*f[i]*delta[i],i)/Sum(delta[i],i);eqkeff_underFPGA:=k=(kC*(1+f)*LzC+kF*LzF+kN*LzN)/Lz;eqkeff_underFPGA_:=k=subs(kC=k,datC,kF=k,datF,kN=k,datN,dat,rhs(%));

pero como la FPGA sólo cubre la mitad de la profundidad en 'y', habrá que tomar la semisuma, y el perfil de conductividades k(x) sería:

>	eqkeff_underFPGA_:=k=(rhs(eqkeff_alongPCB_)+rhs(eqkeff_underFPGA_))/2;k_:=subs(SI0,piecewise(x<0,rhs(eqkeff_alongPCB_),x>0.05,rhs(eqkeff_alongPCB_),rhs(eqkeff_underFPGA_)));plot(k_,x=-0.025..0.075,k=30..200);

por lo que el nuevo perfil de T(x) sería ahora:

>

eqQi:=k1*A1*Diff(Tx1,x)[x=0]=k2*A2*Diff(Tx2,x)[x=0];eqQi:=subs(x=0,Ly=L,k1*L*Lz*diff(T1x,x)=k*L*Lz*diff(T2x,x));Ti_:=subs(k1=k,eqkeff_alongPCB_,k2=k,eqkeff_underFPGA_,dat,solve(%,Ti));'Ti_'=TKC(%);Tm_:=subs(x=LN,Ti=Ti_,Lx=L,Ly=L,eqkeff_underFPGA_,dat,SI0,T2x)*K_;'Tm_'=TKC(%);T123__:=subs(Ti=Ti_,Lx=L,Ly=L,eqkeff_underFPGA_,dat,SI0,piecewise(x<0,T1x,x>0.05,Tm_,T2x)):plot(T123__-To,x=-0.025..0.075,T_C=30..150);

que como se ve es más plano; la Tmax ha bajado de los 150 ºC anteriores a 112 ºC por la mayor conductividad de la FPGA. Como era de esperar, la Ti (comienzo de la FPGA) es la misma, Ti=90 ºC.

e) Volver a resolver el apartado anterior añadiendo el efecto de la radiación térmica por ambas caras con una carcasa que está a 40 ºC, con el modelo de cuerpos negros linealizado.

Como el problema es lineal, todavía admite solución analítica. Integrando el balance energético de un elemento dx con conducción longitudinal y pérdidas radiativas transversales linealizadas para una Textremo conocida (vale para los dos tramos) y para una Qextremo conocida (para el tercer tramo, aunque en él es phi=0):

>	eqBEdx:=0=k*Ly*Lz*diff(T(x),x,x)+phi*Ly*Lz-2*Ly*h*(T(x)-Tinf);eqTo:=T(0)=Ti;solTi:=dsolve({eqBEdx,T(0)=Ti},T(x));eqQo:=Diff(T,x)[x=(L+LN)/2]=0;solQ0:=dsolve({eqBEdx,D(T)(subs(dat,SI0,(L+LN)/2))=0},T(x));

donde h es el coeficiente convectivo correspondiente a las pérdidas radiativas linealizadas, _C2 una constante a obtener por acuerdo en los tramos, y phi la disipación volumétrica (promediada en 'y').

>	eqh:=h*(T-Tinf)=epsilon*sigma*(T^4-Tinf^4);eqh:=h*(T-Tinf)=4*epsilon*sigma*Tmean^3*(T-Tinf);eqh:=h=4*epsilon*sigma*Tmean^3;Tmean=Ti;Tmean_:=Ti_;h_:=subs(Tmean=Ti_,epsilon=1,dat,rhs(eqh));eqphi:=phi=QN/(LN*Ly*Lz);eqphi_:=subs(Ly=L,dat,%);

Ahora hay tres tramos, porque el de la derecha también emite, y hay que imponer las condiciones en los extremos de los tramos. La T(x) en cada tramo es:

>	T1x:=subs(phi=0,k=k1,_C2=_C21,rhs(solTi));T1x_:=subs(k1=k,eqkeff_alongPCB_,h=h_,dat,SI0,%);T2x:=subs(k=k2,_C2=_C22,rhs(solTi));T2x_:=subs(k2=k,eqkeff_underFPGA_,h=h_,eqphi_,dat,SI0,%);T3x:=subs(phi=0,k=k1,_C2=_C23,rhs(solQ0));T3x_:=evalf(subs(k1=k,eqkeff_alongPCB_,h=h_,dat,SI0,%));

y las 4 incógnitas (Ti,C21,C22,C23) se obtienen imponiento las 4 ecuaciones: T1x(-0,025)=Tb, k1·dT1x(0)=k2·dT2x(0), T2x(LN)=T3x(LN), y k2·dT2x(LN)=k3·dT3x(LN), ya que las de T1x(0)=T2x(0)=Ti y dT3x(0,075) ya se han impuesto en la integración.

>

eq1_:=evalf(subs(x=-0.025,dat,SI0,Tb=T1x_));eq2:=k1*A*Diff('T1x',x)[x=0]=k2*A*Diff('T2x',x)[x=0];eq2_:=evalf(subs(x=0,k1=k,eqkeff_alongPCB_,k2=k,eqkeff_underFPGA_,dat,SI0,k1*diff(T1x_,x)=k2*diff(T2x_,x)));eq3:='T2x'[x=LN]='T3x'[x=LN];eq3_:=evalf(subs(x=LN,dat,SI0,T2x_=T3x_));eq4:=k2*A*Diff('T2x',x)[x=LN]=k1*A*Diff('T3x',x)[x=LN];eq4_:=evalf(subs(x=LN,k1=k,eqkeff_alongPCB_,k2=k,eqkeff_underFPGA_,dat,SI0,k2*diff(T2x_,x)=k1*diff(T3x_,x)));sol_:=solve({eq1_,eq2_,eq3_,eq4_},{Ti,_C21,_C22,_C23});T_:=subs(sol_,piecewise(x<0,T1x_,x>0.05,T3x_,T2x_));plot(T_-To,x=-0.025..0.075,T_C=30..70);Ti_:=subs(sol_,Ti)*K_;'Ti_'=TKC(%);Ti2:='T2x'[x=LN];Ti2_:=evalf(subs(sol_,x=LN,dat,SI0,T2x_))*K_;'Ti2_'=TKC(%);Tend:='T3x'[x=(L+LN)/2];Tend_:=evalf(subs(sol_,x=(L+LN)/2,dat,SI0,T3x_))*K_;'Tend_'=TKC(%);N:=40:Tmax:=max(seq(evalf(subs(x=L*i/N,dat,SI0,T_)),i=1..N))*K_;'Tmax'=TKC(%);

i.e. contando la radiación (linealizada) vemos que la Tmax baja de los 112 ºC antes calulado a unos 70 ºC. También se ve que el tramo 3 emite calor y por eso se va enfriando desde la fuente (el IC) hacia el borde adiabático. Las pendientes en los empalmes no son continuas, pero el flujo de calor sí (Q=kAdT/dx).

f) Resolver numéricamente el apartado anterior.

Lo vamos a hacer en régimen transitorio desde un estado inicial con T(x)=Tb, aunque este no sea de equilibrio si Tinf<>Tb, pero es irrelevante para llegar al estacionario..

Solución numérica. Sea N el número de tramos en que dividimos la chapa completa (N+1 puntos): . Desde i=1 hasta i=N/3+1 está soleado y el resto no.

El balance energético de un elemento discreto es:

>

eqBE:=rho[i]*c[i]*A*Dx*(T[j+1,i]-T[j,i])/Dt=k[i+1]*A*(T[j,i+1]-T[j,i])/Dx-k[i-1]*A*(T[j,i]-T[j,i-1])/Dx+phi[i]*A*Dx-p*h*Dx*(T[j,i]-Tinf);eqTi:=T[j+1,i]=T[j,i]+(Dt/(rho[i]*c[i]))*((k[i+1]*(T[j,i+1]-T[j,i])-k[i-1]*(T[j,i]-T[j,i-1]))/Dx^2+phi[i]-(p*h/A)*(T[j,i]-Tinf));eqT1:=T[j+1,1]=Tb;eqTN1:=T[j+1,N+1]=T[j,N];

Como hay discontinuidades (en phi, en k, y en rho·c), si no se suaviza habría que tomar una discretización fina.

>

N:=20;Dx:=evalf(subs(dat,SI0,L/N));phi:='phi':phi__:=subs(SI0,rhs(eqphi_)):for i from 1 to N+1 do phi[i]:=0;od:for i from N/4+2 to 3*N/4 do phi[i]:=phi__:od:phi[N/4+1]:=phi__/2:phi[3*N/4+1]:=phi__/2:plot([seq([(i-1)*Dx,phi[i]],i=1..N+1)],x=0..0.1,'phi'=0..phi__);k1__:=subs(SI0,rhs(eqkeff_alongPCB_));k2__:=subs(SI0,rhs(eqkeff_underFPGA_));for i from 1 to N+1 do k[i]:=k1__;od:for i from N/4+2 to 3*N/4 do k[i]:=k2__:od:k[N/4+1]:=(k1__+k2__)/2:k[3*N/4+1]:=%:plot([seq([(i-1)*Dx,k[i]],i=1..N+1)],x=0..0.1,'k'=0..k2__);i:='i':eqrc:=rho*c=Sum(rho[i]*c[i]*f[i]*delta[i],i)/Sum(delta[i],i);eqrc1:=rho*c=(rhoC*cC*LzC+rhoF*cF*LzF+rhoC*cC*f*LzC)/Lz;rhoc1:=subs(dat,SI0,(subs(datC,rho*c*LzC*(1+f))+subs(datF,rho*c*LzF))/Lz);eqrc2:=rho*c=(rhoC*cC*LzC+rhoF*cF*LzF+rhoC*cC*f*LzC+rhoN*cN*LzN)/Lz;rhoc2:=subs(dat,SI0,(subs(datC,rho*c*LzC*(1+f))+subs(datF,rho*c*LzF)+subs(datN,rho*c*LzN))/Lz);'rhoc2_=(rhoc1+rhoc2)/2';rhoc2_:=(rhoc1+rhoc2)/2;for i from 1 to N+1 do rc[i]:=rhoc1;od:for i from N/4+2 to 3*N/4 do rc[i]:=rhoc2_:od:rc[N/4+1]:=(rhoc1+rhoc2_)/2:rc[3*N/4+1]:=%:plot([seq([(i-1)*Dx,rc[i]],i=1..N+1)],x=0..0.1,'rc'=0..rhoc2);

donde ya se ha tenido en cuenta para el producto rho·c (rc=rho·c), al igual que para la conductividad k, que la FPGA sólo cubre la mitad de la profundidad en y,

La discretización temporal es (Note: to simulate case d (h=0) use tsim=750 s and M=3500).

>	T:='T':tsim:=300;M:=1500;Dt:=evalf(tsim/M);Fo_:=k1__*Dt/(rhoc1*Dx^2);Lz_:=subs(dat,SI0,Lz);h_:=subs(SI0,h_);Tinf_:=subs(dat,SI0,Tinf);T0_:=subs(dat,SI0,Tb);T:=Matrix(1..M+1,1..N+1,T0_):

Y la simulación espacio-temporal:

>

for j from 1 to M do T[j+1,1]:=T0_;T[j+1,N+1]:=T[j,N];for i from 2 to N do T[j+1,i]:=T[j,i]+(Dt/(rc[i]))*(((k[i+1]+k[i])*(T[j,i+1]-T[j,i])-(k[i]+k[i-1])*(T[j,i]-T[j,i-1]))/(2*Dx^2)+phi[i]-(2/Lz_)*h_*(T[j,i]-Tinf_));od:od:i:='i':j:='j':Mjump:=50:sx:=seq([seq([i*Dx-Dx,T[1+j_*Mjump,i]-To],i=1..N+1)],j_=0..M/Mjump):plot([sx],x_m=0..0.1,T_C=30..70);st:=seq([seq([j*Dt,T[j,i]-To],j=1..M+1)],i=1..N+1):plot([st],t_s=0..tsim,T_C=30..70);Tmax:=max(max(T))*K_;'Tmax'=TKC(%);

Vemos que la T(x) final coincide con la del modelo analítico, y que se tarda 2 o 3 minutos en alcanzar prácticamente el régimen estacionario.

Puede repetirse la simulación numérica sin linealizar la radiación, i.e. substituyendo h·(T-Tinf) por sigma·(T^4-Tinf^4), obteniéndose un perfil T(x) similar (la Tmax aumentaría de 70 ºC a 74 ºC).

g) Resolver el caso bidimensional estacionario.

Se ha resuelto en Matlab y el resultado es:

Como era de esperar, al estar más concentrada la disipación, la Tmax aumenta bastante, de 70 ºC a 87 ºC, por lo que el modelo 1D no es muy bueno.

Si nos fijamos en el perfil de temperaturas en los bordes aislados, e.g. en T(x,y=0), podría parecer erróneo que el modelo 1D no sirva de aproximación al modelo 2D (la Tmax era 70 ºC en el 1D y vale 80 ºC en el 2D), i.e. la T en el borde es mayor cuanto más alejada del borde está la disipación. La explicación es que, manteniendo fija la disipación total (10 W), cuanto más concentrada esté más subirá la T en la FPGA, y como hacia los lados laterales fluye poco calor (los bordes son adiabáticos), toda la tarjeta PCB a una cierta cota x aumenta su temperatura, no solo la parte bajo la FPGA.

>