master
John Lauer Add native thermal-pour expansion experiment, finite-difference comparison and recording scripts 33c4ac5 1mo ago
"""First-order 2D copper-sheet conduction plus bottom-surface convection/radiation.
Not a package, whole-board, enclosure, or junction-temperature simulation.
Source copper held at 65 C, ambient/radiative surroundings at 25 C.
"""
import json,sys
import numpy as np
from scipy.sparse import coo_matrix,diags
from scipy.sparse.linalg import spsolve
from scipy.ndimage import label
from shapely import contains_xy
from shapely.geometry import Point
from shapely.ops import unary_union
from analyze import P,regions,copper
SIGMA=5.670374419e-8
BASE=P.parent/'thermal/recording-02/bq25792-astra-thermal-routed.brd'
FINAL=P/'expansion-03/expansion-03-filled.brd'
def ground(file):
 g=copper(file)['PGND'];src=unary_union([Point(x,y).buffer(.225) for x,y in [(2.05,10.3),(2.8,11.05),(35.2,12.3),(35.2,12.9),(35.2,14)]])
 parts=list(g.geoms) if hasattr(g,'geoms') else [g];g=unary_union([p for p in parts if p.intersects(src)])
 return {'geometry':g,'sources':src,'areaMm2':g.area,'components':['IC1','U2']}
def mesh(g,src,step):
 x0,y0,x1,y1=g.bounds
 xx,yy=np.meshgrid(np.arange(x0+step/2,x1,step),np.arange(y0+step/2,y1,step))
 present=contains_xy(g,xx,yy);fixed=contains_xy(src,xx,yy)&present
 labs,n=label(present);keep=np.unique(labs[fixed]);present &= np.isin(labs,keep[keep>0]);fixed &= present
 ids=np.full(present.shape,-1,dtype=int);ids[present]=np.arange(present.sum());N=int(present.sum());bc=fixed[present]
 edges=[]
 for ax in [0,1]:
  a=ids[:-1,:] if ax==0 else ids[:,:-1];b=ids[1:,:] if ax==0 else ids[:,1:];ok=(a>=0)&(b>=0);edges.extend(zip(a[ok],b[ok]))
 a,b=np.array(edges).T
 # k=390 W/m/K, native stackup copper=35 um. Per square sheet conductance k*t.
 conductance=390*35e-6
 rows=np.concatenate([a,b,a,b]);cols=np.concatenate([a,b,b,a]);vals=np.concatenate([np.ones(len(a)),np.ones(len(a)),-np.ones(len(a)),-np.ones(len(a))])*conductance
 L=coo_matrix((vals,(rows,cols)),shape=(N,N)).tocsr();free=~bc
 rhs=-L[free][:,bc].sum(axis=1).A.ravel()*40
 return L,free,rhs,N,step**2*1e-6,ids,present,bc

def solve(m,h,eps,mask=True):
 L,free,rhs,N,A,ids,present,bc=m;T=np.full(N,40.);Ta=298.15
 # Film resistance only for masked areas: provisional 25um / 0.2W/mK.
 film=25e-6/.2 if mask else 0.
 for _ in range(12):
  Ts=Ta+T;hr=eps*SIGMA*(Ts+Ta)*(Ts*Ts+Ta*Ta)
  coeff=1/(1/(h+hr)+film);H=coeff*A
  nxt=spsolve(L[free][:,free]+diags(H[free]),rhs)
  delta=np.max(abs(nxt-T[free]),initial=0.);T[free]=nxt
  if delta<1e-5:break
 Q=float(np.sum(H*T));iso=(1/(1/(h+eps*SIGMA*((Ta+40)+Ta)*((Ta+40)**2+Ta**2))+film))*A*N*40
 # source energy balance, including surface heat lost at fixed-temperature source cells.
 inputQ=float((L@T)[bc].sum()+(H*T)[bc].sum())
 assert abs(inputQ-Q)<1e-6,(inputQ,Q)
 heat=np.full(ids.shape,np.nan);heat[present]=T
 return {'bottomHeatW':Q,'conductanceWPerK':Q/40,'isothermalUpperW':iso,'sheetEfficiency':Q/iso,'energyResidualW':inputQ-Q,'meshAreaMm2':N*A*1e6},heat

def main():
 result={'assumptions':{'sourceCopperC':65,'ambientC':25,'copperUm':35,'copperConductivityWmK':390,'maskUm':25,'maskConductivityWmK':.2,'maskEmissivity':.8,'brightCopperEmissivity':.05,'oxidizedCopperEmissivity':.8,'nominalConvectionWm2K':5,'scope':'Bottom-only copper sheet; no FR4/top/internal conduction, package/via resistance, enclosure, neighboring component heat or actual load power. Existing ground copper is reassigned: do not add these W values to whole-board capability.'},'boards':{}}
 for name,file in [('before',BASE),('expanded',FINAL)]:
  rs=regions(file);rs['PGND']=ground(file);result['boards'][name]={}
  for n,r in rs.items():
   m=mesh(r['geometry'],r['sources'],.1);q,heat=solve(m,5,.8);entry={'areaMm2':r['areaMm2'],'masked':q,'bright':solve(m,5,.05,False)[0],'oxidized':solve(m,5,.8,False)[0],'h3':solve(m,3,.8)[0],'h10':solve(m,10,.8)[0]}
   # Half-area exposed finish approximation via area-weighted emissivity, explicitly not a shaped mask window.
   entry['halfBrightApprox']=solve(m,5,(.8+.05)/2,False)[0]
   entry['coarseGrid']=solve(mesh(r['geometry'],r['sources'],.15),5,.8)[0]
   result['boards'][name][n]=entry
   np.savez_compressed(P/f'{name}-{n}-temperature.npz',temperatureRiseC=heat,extent=r['geometry'].bounds)
   print(name,n,round(r['areaMm2'],2),round(q['bottomHeatW'],4),'W',flush=True)
 (P/'thermal-results.json').write_text(json.dumps(result,indent=2))
if __name__=='__main__':main()