Cho and Munro: Kilovision: Kilovoltage x-ray source …. 10
Appendix A: THEORY FOR FINITE ELEMENT MODELING
There are two general approaches to thermal calculations when modeling a structure containing both solids and fluids. Conjugate heat-fluid transfer calculations estimate both energy and mass transfer in the solid and the fluid. This form of calculation requires that the user have substantial computer resources, since a large number of finite elements and calculation iterations are required for the calculations to converge. The alternative method separates fluid flow from heat transfer. The fluid flow simulation calculates the flow rates and flow distributions in the cooling channels, which is the most important information needed to estimate the convection heat transfer at the solid/fluid boundary as a function of wall temperature. The pre-calculated estimates of the convection heat transfer are used as inputs into the thermal model. This approach has the advantage that we can account for the wide range of the convection heat transfer due to the phase change of the cooling fluid, such as nucleate boiling, fully developed boiling and critical heat flux. In addition, the approach substantially reduces the calculations times and computer resources needed for the calculations. The only disadvantage is that this type of modeling simplifies the geometry of the object. In reality heat transfer depends upon flow rate and cooling tube geometry, which can vary in complex ways from location to location in the object. Assigning the heat transfer characteristics of the solid/fluid boundary tends to simplify the variations in convection heat transfer over the surfaces of the cooling tubes. Because of limitations in our computer resources and in the number of nodes that could be used in our academic version of the finite element modeling software, it was necessary that we use the second calculation approach for our calculations.
A. Heat Transfer Mechanisms and Thermal Modeling
There are three heat transfer mechanisms: conduction, convection and radiation. In our thermal models, conduction was modeled, convection was accounted for using the techniques described below and radiative heat transfer was ignored. Our preliminary calculations suggested that the maximum temperature reached in the model would not exceed 1600 K (750W for megavoltage target and 1200W for kilovoltage target). The radiative heat flux is less than 7.7 W/cm2 at this temperature{Yih & Wang 2001 YIH2001 /id} and can thus be safely ignored.
The convective heat transfer is a very complex phenomenon. As the wall temperature increases, the processes that govern heat transfer pass through a number of distinct stages: forced convection heat transfer (qfc), nucleate boiling, fully developed boiling (qb), and critical heat flux (qchf). When the wall temperature of the cooling tube is below the boiling point of the fluid, the heat transfer is dominated by forced convection heat transfer, which is governed by the thickness of a stationary layer of fluid close to the wall, which in turn is due to the viscosity of the fluid. This thin layer of fluid reaches an elevated temperature and acts as an insulating layer, reducing the heat transfer across the solid/fluid boundary. Forced convection heat transfer is a function of thermal properties and velocity of the coolant, the shape of cooling path and the temperature difference between cooling wall and bulk of the fluid. We have used the equation for forced convection suggested by Bjorge et. al.,{Bjorge, Hall, et al. 1982 BJORGE1982 /id} which applies to fully developed turbulent flow in sub-cooled and regions of low thermal quality [The thermal quality is defined by the vapor content of a liquid. When all of the liquid is changed into vapor at the boiling temperature, the quality is 1.0. When the liquid and the vapor are mixed, the thermal quality is between 0 to 1. When the temperature of liquid is lower than its boiling temperature, the thermal quality has a negative value.]:
qfc = 0.023(GD/μ)0.8 (μCp /k)0.3 k/D (T-Tb) (1)
where G, D, μ, Cp, k, T, and Tb are mass flux, tube diameter, dynamic viscosity, specific heat, thermal conductivity, temperature of the wall and temperature of the bulk of the coolant, respectively. The dynamic viscosity was evaluated at the film temperature, which is the average temperature of T and Tb .
As nucleate boiling starts, the thin layer of fluid is converted into a gas (boiling bubble) that is much more efficient at transferring heat into the bulk of the fluid - due to the agitation created by bubble motion. Fully developed boiling depends mainly on the thermal properties and the applied pressure of the coolant. We have used the equation derived by Rohsenow to predict the heat flux as a function of wall temperature.{Bjorge, Hall, et al. 1982 BJORGE1982 /id}
qb = A ·ΔTsat3 (2)
where ΔTsat is the temperature difference with respect to saturation temperature. The “A” is a complex coefficient dependent upon the properties of water.{Bjorge, Hall, et al. 1982 BJORGE1982 /id} These properties (density, thermal conductivity, specific heat, enthalpy of vaporization, surface tension, viscosity, and saturation temperature) have been obtained from the literature.{Gray 1972 GRAY1972 /id}
When the bubble generation becomes too great, it may generate an insulating vapor blanket and prevent the cooling liquid from reaching the hot surface. The heat flux where the bubbles start to aggregate and unstable cooling conditions starts to occur is called the critical heat flux. It depends on many parameters such as the geometry of cooling path, surface finish, coolant velocity, pressure, and the size of the heat source. We have used an equation derived by Inasaka et.al. to predict the critical heat flux.{Inasaka & Nariai 1997 INASAKA1997 /id}
qchf = C Hfg G 0.4 µ 0.6 D -0.6 (3)
where: Hfg is latent heat of evaporation (J/kg), G is mass velocity (kg/m2s), µ is viscosity of liquid, D is the tube inside diameter, and C is the correction factor proposed by Inasaka as a function of pressure, mass velocity, and exit equilibrium quality of the fluid. Equation (3) applies to low-pressure (0.1-7.0MPa), medium flow-rate (5.5-30Mg/m2s) and long heating length (10cm) conditions, which closely match our physical situation (0.3MPa, ~8Mg/m2s).{Nariai, Inasaka, et al. 1987 NARIAI1987 /id}
The total heat flux curve is formed by combining the three curves from each heat transfer stage: forced convection heat transfer, fully developed boiling heat transfer, and critical heat flux. We have used the mathematical relationship suggested by Bjorge et.al.{Bjorge, Hall, et al. 1982 BJORGE1982 /id} to combine the heat transfer curves for forced convection and fully developed boiling.
qr = (q2fc + q2b[1-w3])1/2 (4)
w=(T-Tib)/(T-Tsat)
where T, Tib and Tsat are the temperature of solid surface, incipient boiling temperature, and saturation temperature of liquid, respectively. To extend the temperature range to the region of critical heat flux, the following equation has been used in this paper:
qtot = (q-2r + q-2chf [1-w3])–1/2 (5)
This generates a smoothly changing curve of heat flux as a function of temperature of the solid/fluid boundary. Equation (5) can be used to generate a table of values for the total heat flux versus wall temperature, which can then be used in the finite element calculations to represent the convection heat transfer coefficients.
B. Thermal-Mechanical Modelling and Thermal Damage
Once the thermal calculations have been completed, a temperature distribution throughout the model is known. From this temperature distribution, the stress and strain in the model can be calculated if the mechanical properties of the materials forming the model are known. These quantities (coefficient of thermal expansion, Young’s modulus, Poisson’s ratio, strength coefficient, and hardening exponent) have been obtained from the literature.{White & Minges 1997 WHITE1997 /id}{Yih & Wang 2001 YIH2001 /id}{Upthegrove & Burghoff 1953 UPTHEGROVE1953 /id} Typical finite element calculations yield the material behavior for monotonic loading conditions (e.g., linear stress/strain quantities). Unfortunately, the material behavior under cyclic loading conditions (e.g., repeated thermal loading and unloading) is the quantity needed to estimate the life-time of our target assemblies. Thus, we need to convert from monotomic to cyclic material properties. One method of converting between linear stress/strain (monotonic material properties) and nonlinear cyclic stress/strain is Neuber’s rule.{Neuber 1961 NEUBER1961 /id} We have used Neuber’s rule to convert linear monotonic stress/strain into cyclic stress/strain using the equations below:
ε = σ/E +(σ/K’)1/n’ (6)
ε•σ = (Kσ S) • (Kε e) (7)
where: K’ is the cyclic strength coefficient and n’ is the cyclic hardening exponent. Those can be calculated from the material properties such as Young’s modulus (E), ductility (εf or true fracture strain), and ultimate strength(σu).{Muralidharan & Manson 1988 MURALIDHARAN1988 /id} Kσ is elastic stress concentration factor, Kε is the elastic strain concentration factor, S is the engineering or nominal stress and e is the engineering or nominal strain. The (Kσ S) and (Kε e) represent the true stress and strain of interest, respectively and these quantities are output by the linear finite element analysis. Equations (6) and (7) can be solved for cyclic stress (σ), or cyclic strain (ε) by iteration or by various numerical techniques.
There are a number of approaches to convert linear into cyclic stress/strain which are not as conservative as Neuber’s rule, such as the use of a fatigue notch factor instead of elastic concentration factor,{Topper, Wetzel, et al. 1969 TOPPER1969 /id} use of a modified relationship between linear and cyclic stress-strain, or use of a strain energy density rule. However, we have used Neuber’s rule since it represents the most conservative technique of converting between linear and cyclic stress/strain and thus likely to lead to safe estimates of the maximum thermal load.
When a stress load is removed, the material behavior differs from that under stress loading. This is known as the Bauschinger effect. Materials for which the hysteresis loop (between loading and unloading stress) can be described by magnifying the cyclic stress/strain curve by a factor of 2 are said to exhibit “Masing-type” behavior. This type of behavior is typical of many common metals.{Stephens 2001 STEPHENS2001 /id} All the materials in our model are assumed to exhibit Masing-type behavior. The following equations are used to model the material deformation when heat loading is removed.
Dε = Dσ/E +2• (Dσ/2K’)1/n’ (8)
Dε•Dσ = (Kσ DS) • (Kε De) (9)
where the symbols are the same as in equations (6) and (7). Knowing the strain range Dε, we can obtain the strain amplitude, εa = Dε/2, for the fatigue analysis.
Once the strain amplitude is known, the modified universal slope method{Muralidharan & Manson 1988 MURALIDHARAN1988 /id} can be used to estimate the number of loadings to fatigue failure:
εa = 0.623(σu/E)0.832(2Nf)-0.09 + 0.0196(εf)0.155(σu/E)-0.53(2Nf)-0.56 (10)
where σu is ultimate strength, εf is true fracture strain (or ductility), and Nf is the number of loading before fatigue failure. This equation can be solved for the number of loading to fatigue failure, Nf, by an iterative or numerical method, given the strain amplitude, εa.
In most situations variable amplitude loading occurs, where the object is exposed to different stress loads (i.e., different heat loading) for different operating modes. When variable amplitude loading occurs, determining the fatigue life of a component requires a method that determines how each load condition contributes to the overall fatigue of the object. One such method is Palmgren-Miner’s damage rule, where accumulation of damage is assumed to be linear.{Miner 1945 MINER1945 /id}
(11)
where ni is the number of cycles of damage and Ni is the life at an i-th level of amplitude. The fatigue failure is expected to occur when the sum of the damage becomes 1.0. The reciprocal of Ni provides a measure of the damage at an i-th level of loading amplitude.
Appendix B: SUPPLEMENTAL MATERIALS AND METHODS
A. Heat Flux Calculations
The fluid flow in the cooling channels can be quite complex. Therefore, to get a good estimate of G, the mass flow rate, (see Eqs. 1 and 3), we first modeled fluid flow in the cooling tubes using ANSYS 5.6. Due to license limitations our software could only model the cooling tubes and insufficient nodes were available to model the rest of the megavoltage target assembly. The fluid velocity across the surface of the inlet of the cooling tubes was assigned to the finite elements assuming that the bulk fluid flow rate is 60 ml/s (1 gallon/minute{Varian Medical Systems 1999 VARIANMEDICALSY1999 /id}). This is a minimum value and so should represent a conservative estimate of the flow. The inlet flow was assumed to be uniform across the opening except for the elements right at the wall of the tube where the flow was assumed to be zero. Since, fully developed turbulent flow occurs quite rapidly in the model,{Rohsenow 1973 ROHSENOW1973 /id} and our flow model has a long introductory inlet, the boundary conditions have negligible effects on the final results. The result of the calculations is the flow rate at all locations in the cooling tube. The average flow rate (flow velocity) times the density of the fluid gives G, the mass flow rate. This numerical value has been used to calculate the heat flux for that region of the cooling tube. In one of our models, the flow rate on one side of the cross-over tube is different from the flow rate on the other side. For this situation the mass flow rate for each half of the tube has been calculated independently.