nano HUB .org online simulations and more Modeling and Simulation of Sub-Micron Thermal Transport Jayathi Y. Murthy School of Mechanical Engineering Purdue University
nano HUB .org online simulations and more Outline • Review of phonon transport • Phonon Boltzmann Transport Equation (BTE) • Finite volume method • Improved BTE models • Future directions and closure
nano HUB .org online simulations and more Sub-micron Heat Conduction L L Λ <<L Λ >>L (a) (b) Two regimes, (a) Fourier’s law is valid, (b) Fourier’s law is invalid, Λ is the mean free path of the heat carrier (phonons or electrons) Phonons are quanta of lattice vibrations. They are the main carriers of energy in semiconductors and dielectrics.
nano HUB .org online simulations and more Heat Transport in Solids m 1 m 2 One-dimensional spring-mass system x n y n x n+1 2 d x m n g y ( y 2 x ) = + � 1 n n 1 n 2 + dt 2 d y n m g x ( x 2 y ) = + � 2 n 1 n n 2 + dt Dispersion relation for crystal vibrations, Majumdar (1998)
Boltzmann Transport nano HUB .org online simulations and more Equation Boltzmann transport equation for phonons: f f � � � � ( v f s ) + �� = � � t t � � � � scat . Ballistic term Phonons are characterized by ( r , t, K ) and polarization Frequency vs. reduced wave number in (100) direction for silicon
nano HUB .org online simulations and more Scattering Processes • Phonons scatter on impurities, grain boundaries, isotopes, physical boundaries, other phonons and carriers • Three-phonon interactions are major contributors to thermal resistance at room temperature in silicon • Three-phonon interactions must satisfy momentum and energy conservation rules:
nano HUB .org online simulations and more Three-Phonon Interactions General scattering term for 3-phonon processes very complex: Valid only for phonons satisfying conservation rules
nano HUB .org online simulations and more Relaxation Time Approximation f f � � � � ( v f s ) + � = � � t t � � � � scat . Planck distribution function f f ( K ) f ( r,K ) � � � � 0 = � � t ( K ) � � � � scat . Relaxation time
nano HUB .org online simulations and more Gray Phonon BTE 0 e e e • Energy form: �� �� � � ( ) s v e �� + �� = g t � � eff • e ” is energy per unit volume per unit solid angle • There are as many pde’s as there are s directions. In each direction, e ” varies in space and time • Different directions are related to each other because of e 0 in the scattering term: 1 1 0 e ( , ) r t e d e sin d d �� �� = � � = �� � � � 4 4 � � 0 e e d • Notice that �� � 0 � = � No net � eff 4 � energy source
nano HUB .org online simulations and more When is a Gray Model Good? • BTE in relaxation time approximation incapable of transferring energy across frequencies and polarizations – Not useful in MOSFET modeling • Too expensive at present to solve full BTE, though work is in progress • Gray approximations to BTE are good for problems with small departures from equilibrium provided modeling is done well Courtesy IBM
Overview of Finite Volume nano HUB .org online simulations and more Method • Divide spatial domain into control volumes of extent Δ x Δ y • Divide angular domain into control angles of extent ΔΩ s • Divide time into steps of Δ t - but will only do steady here for simplicity • Integrate gray BTE in each direction over control volume and control angle. Get discrete equation set. • Solve each direction sequentially and iteratively • Back out “temperature” from e 0 upon convergence using 0 4 e � T T = + ref C
nano HUB .org online simulations and more Discretization • Divide domain into rectangular control volumes of extent Δ x and Δ y. Assume 2D, so that depth into page (z) is one unit. z y s x y x • Divide angular domain of 4 π in N θ xN φ control angles per octant. Centroid of each control angle is ( θ i , φ i ), extents are ( Δθ , Δφ ). For each control angle i: s sin sin i + sin cos j+ cos k = � � � � � i i i i i i • Important: The directions s are 3D even though we are considering 2D
Angular Discretization nano HUB .org online simulations and more Nomenclature • Control angle extent is /2 /2 � �� � � + � � i i /2 sin d d 2sin sin 0.5 ( ) �� = � � � � � = � � � � � i /2 � �� � � �� � i i • In 2D, only directions in the “front” hemisphere are necessary. Thus θ ranges from 0- π /2 and φ =0-2 π • Thus, increase control angle extent to: /2 /2 � �� � � + � � i i 2* sin d d �� = � � � � � = /2 /2 � �� � � �� � i i Weight _ Factor *2sin sin 0.5 ( ) � � � � � i Weight _ Factor 2 for 2D = • Define for future use: /2 /2 � �� � � + � � = � S i i s d d � � � /2 /2 � �� � � �� � i i
nano HUB .org online simulations and more Formula for S sin sin 0.5 cos 2 sin i ( ) ( ) � � � � � � � � � � � � � � � i i � � S ( ) Weight _ Factor + cos sin 0.5 cos 2 sin j ( ) � � � � = � � � � � � � � � i � j � � � + 0.5 sin 2 sin k � � � � � � � i � �
nano HUB .org online simulations and more Spatial Discretization N n W E P e” stored w e at cell Δ y centroids s Δ x S
nano HUB .org online simulations and more Control Volume Balance • Integrate governing equation over control volume and control angle: 0 " e e s � ( ) s v e dVd i i dVd �� �� �� � = �� � i g n � eff V , V , � �� � �� Look at LHS. Apply Divergence Theorem: P � � � � � � w e " " v e s n dA d v e n dA s d � � � � = � � � � � � � � � � g i i g i i � � � � � � A A �� � � �� � �� s S " = v e n dA � � � faces f i g i � � A � � �� S=S S=S " " v e n dA v e A n � � � � � ii ii g i g if f f � � faces A � SSSS SSSS " " " " = v e y v e y v e x v e x � � � + � � � g ie xi g iw xi g in yi g is yi
Control Volume Balance nano HUB .org online simulations and more (cont’d) • Now look at RHS 0 " 0 " e e e e � � ( 2 ) i i dVd iP iP V O x �� � = � �� + � � � eff eff V , � �� • Collecting terms: SSSS SSSS " " " " v e y v e y v e x v e x � � � + � � � g ie xi g iw xi g in yi g is yi 0 " e e � i P , i P , V = � �� � eff • Control volume balance says that net rate of energy entering the CV in direction s i must be balanced by net in-scattering to the direction i in the CV
Upwinding nano HUB .org online simulations and more • e” is stored at cell centroids, but we need it on the CV faces • Need to interpolate from cell centroid to face • Can use a variety of schemes to perform interpolation – Central difference scheme 0.5 on uniform mesh ( ) � = � + � e P E • Second-order accurate, but wiggles in spatial solution – Upwind difference scheme S if >= 0 � = � E e P x P W S if < 0 = � e E x • Upwind scheme is only first order accurate • Higher order schemes are available
nano HUB .org online simulations and more Discrete Equation • Using upwinding and collecting terms, we obtain an algebraic equation: " " a e a e b � = + i P i P , , i nb i nb , , i P , nb Here, nb are the spatial neighbors E,W,N,S. • We obtain one such equation for each grid point P for each direction i • The b term contains e 0iP • Once we have boundary conditions discretized, we can solve the set
nano HUB .org online simulations and more Discussion • Prefer to solve iteratively and if possible, sequentially to keep memory requirements low • For upwind scheme, diagonal dominance is guaranteed, making it possible to use iterative schemes • Conservation of energy is guaranteed regardless of spatial and angular discretization – Confirm that sum of all scattering source terms at a point is zero regardless of discretization • Any linear solver can be used
Overall Solution nano HUB .org online simulations and more Algorithm 1. Initialize all e ” i values for all cell centroids and directions 2. Find e 0 P for each point P from current e ” values. 3. Start with direction i=1 4. For direction i: – Find discretized equations for direction i, assuming e 0 temporarily known – Solve for e ” i at all grid points using a linear solver 5. If (i.le. total number of directions) go to 4 6. If (all directions are done) check for convergence. If converged, stop. Else, go to 2.
Recommend
More recommend
