Mixing Two Solutions of Different Initial Temperature for a Continuous Variations Experiment

Linear System Applications

William Ford , in Numerical Linear Algebra with Applications, 2015

12.2 Finite Difference Approximations

Differential equations are at the heart of many engineering problems. They appear in fluid mechanics, thermodynamics, vibration analysis, and many other areas. For most of the problems that appear in practice, we cannot find an analytical solution and must rely on numerical techniques. Applications of linear algebra appear particularly in what are termed initial-boundary value problems. In problems like these, an initial condition is known at t = 0, and the value of the solution is known on the boundary of an object. We must use these known values to approximate the solution in the interior of the object.

12.2.1 Steady-State Heat and Diffusion

Suppose a thin rod is given an initial temperature distribution, then insulated on the sides. The ends of the rod are kept at the same fixed temperature; e.g., suppose at the start of the experiment, both ends are immediately plunged into ice water. We are interested in how the temperature along the rod varies with time. Suppose that the rod has a length L (in meters), and we establish a coordinate system along the rod as illustrated in Figure 12.3. Let u (x, t) represent the temperature at the point x meters along the rod at time t (in seconds). We start with an initial temperature distribution u (x, 0) = f (x). This problem is modeled using a partial differential equation called the heat equation:

Figure 12.3. The heat equation: a thin rod insulated on its sides.

$\begin{array}{ll}\frac{\partial u}{\partial t}& =c\frac{{\partial }^{2}u}{\partial x^2},0\le x\le L,0\le t\le T,\\ u(0,t)& =u(L,t)=0,\\ u(x,0)& =f(x).\end{array}$

The constant c is the thermal diffusivity. For most functions f (x), we cannot obtain an exact solution to the problem, so we must resort to numerical methods. So where does linear algebra come in? We divide the space interval 0 ≤ x ≤ L into small subintervals of length h = L/m, and the time interval into small subintervals of length k = T/n, where m and n are integers (Figure 12.4). We wish to approximate u (x, t) at the grid points $\left(x_i,t_j\right)$ , where (x ₁, x ₂, x ₃, …,x _m , x _m+1) = (0, h, 2h, …, L − h, L) and (t ₁, t ₂, t ₃, …, t _n , t _n+1) = (0, k, 2k, …, T − k, T). We denote this approximate solution by

Figure 12.4. Numerical solution of the heat equation: subdivisions of the x and t axes.

$u_{ij}\approx u(x_i,t_j).$

If m and n are large, the values of h and k are small, and we can approximate each derivative at a point $\left(x_i,t_j\right)$ using finite difference equations as follows:

(12.4) $\frac{\partial u}{\partial t}\left(x_i,t_{j+1}\right)\approx \frac{u_{i,j+1}-u_{ij}}{k},1\le j\le n,$

(12.5) $\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{j+1}\right)\approx \frac{1}{h^2}\left(u_{i-1,j+1}-2u_{i,j+1}+u_{i+1,j+1}\right),2\le i\le m.$

These approximations are developed using Taylor series for a function of two variables, and the interested reader can refer to Ref. [33, Chapter 6], for an explanation. By equating these approximations, we have

$\frac{u_{i,j+1}-u_{ij}}{k}=\frac{c}{h^2}\left(u_{i-1,j+1}-2u_{i,j+1}+u_{i+1,j+1}\right).$

All but one term in the finite difference formula involves j + 1, so isolate it to obtain

(12.6) $u_{i,j}=-r{u}_{i-1,j+1}+(1+2r)u_{i,j+1}-r{u}_{i+1,j+1},2\le i\le m,1\le j\le n,$

where r = ck/h ². Notice that Equations 12.4 and 12.5 specify relationships between points in the grid pattern of Figure 12.5. Figure 12.6 provides a view of the entire grid. In Figure 12.6, the black circles represent boundary and initial values and the open circles represent a general pattern of the four points used in each difference equation.

Figure 12.5. Numerical solution of the heat equation:locally related points in the grid.

Figure 12.6. Grid for the numerical solution of the heat equation.

Using matrix notation, Equation 12.6 becomes

(12.7) $u_j=B{u}_{j+1}-r{b}_{j+1},$

where

$B=\left[\begin{array}{ccccc}1+2r& -r& 0& \cdots & 0\\ -r& 1+2r& -r& \cdots & 0\\ 0& \ddots & \ddots & \ddots & ⋮\\ ⋮& ⋮& -r& 1+2r& -r\\ 0& 0& 0& -r& 1+2r\end{array}\right]$

and b _j+1 accounts for the boundary or initial conditions. Now, we can write Equation 12.7 as

(12.8) $B{u}_{j+1}=u_j+r{b}_{j+1}.$

This is a linear algebraic system to solve for u _j+1. Most of the elements of B are 0, so it is a sparse matrix; furthermore, it is tridiagonal. The Thomas algorithm presented in Section 9.4 takes advantage of the tridiagonal structure and solves system 12.7 in O (n) flops as opposed to $O\left(n^3\right)$ . Think of Equation 12.8 this way. Values of u _j and rb _j+1 give us values for u _j+1, and we work our way up the grid. This is not a "chicken and an egg" problem, since the right-hand side of Equation 12.8 is known at time t = 0. We now give an example of our finite difference equations in action by approximating and graphing the solution to the heat flow problem

(12.9) $\begin{array}{ll}\frac{\partial u}{\partial t}& =0.875\frac{{\partial }^{2}u}{\partial x^2},0\le x\le 5.0,0\le t\le 1.0,\\ u(0,t)& =u(5,t)=0,\\ u(x,0)& =e^{x}sin(\pi x).\end{array}$

The required MATLAB code, heateq.m, is located in the software distribution. We present a graph of the approximate solution obtained in Figure 12.7.

Figure 12.7. Graph of the solution for the heat equation problem.

We will deal with more problems like this in Chapters 20 and 21, where we discuss the Poisson and biharmonic equations.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123944351000120

Differential Equations

Bastian E. Rapp , in Microfluidics: Modelling, Mechanics and Mathematics, 2017

8.3.3.11 Solution

Combining Eq. 8.80 with Eq. 8.78 allows us to find the general solution to the one-dimensional heat equation that is given by

(Eq. 8.81) $f_{\lambda <0}\left(x,t\right)={\displaystyle \sum _{n=1}^{\infty }\frac{2}{L}}\left({\int }_0^L{f}_0\left(x\right)sin\left(\frac{n\pi }{x}L\right)dx\right){\mathrm{e}}^{-{\left(\frac{n\pi c}{L}\right)}^{2_t}}sin\left(\frac{n\pi }{x}L\right)$

Solution for the Special Case f ₀ (x)   = T ₀

In our original description of the heat equation we assumed the initial temperature distribution to be a constant T ₀. We will now simplify Eq. 8.81 for this special case. First, we solve the integral

$\begin{array}{lll}{\int }_0^L{T}_{0}sin\left(\frac{n\pi }{x}L\right)dx& =& T_0\frac{L}{n\pi }{\left[-cos\left(\frac{n\pi }{x}L\right)\right]}_0^L\\ & =& T_0\frac{L}{n\pi }\left(-cos\left(n\pi \right)+1\right)\end{array}$

As you can see, for all n  =   0, 2, 4, … the value for cos (n π)   =   1 in which case the bracket term is 0 and thus the whole function will be zero. For all n  =   1, 3, 5, …, cos (n π)   =   −   1 and the bracket term will amount to 2. In this case the integral evaluates to

${\int }_0^L{T}_{0}sin\left(\frac{n\pi }{x}L\right)dx=2T_0\frac{L}{n\pi }$

This allows us to simplify Eq. 8.81 to

(Eq. 8.82) $\begin{array}{lll}T\left(x,t\right)& =& {\displaystyle \sum _{n=1}^{\infty }\frac{2}{L}}2T_0\frac{L}{\left(2n+1\right)\pi }{\mathrm{e}}^{-{\left(\frac{\left(2n+1\right)\pi c}{L}\right)}^{2}t}sin\left(\frac{\left(2n+1\right)\pi }{x}L\right)\\ & =& \frac{4T_0}{\pi }{\displaystyle \sum _{n=0}^{\infty }\frac{1}{\left(2n+1\right)}}{\mathrm{e}}^{-{\left(\frac{\left(2n+1\right)\pi c}{L}\right)}^{2}t}sin\left(\frac{\left(2n+1\right)\pi }{x}L\right)\end{array}$

where we have substituted 2n  +   1 for n in order to only account for odd values of n. As you can see Eq. 8.82 consists of two terms. The sine term will give the local distribution of the temperature at a given time point t as it is not dependent on the time. The second term is an exponential function that depends on the time t and on the constant c, which is (if you reconsider the equation given by Eq. 8.3) the thermal diffusivity, i.e., a measure of how quickly the material can transport heat. The exponential term only depends on the physical parameters, e.g., the length of the rod L but not on the independent variable x, i.e., the location. The exponential function makes the temperature profile decay over time while the sine sum will create the distribution.

Fig. 8.6 shows plotted results for the simplified one-dimensional heat equation, Eq. 8.82, for two different materials: stainless steel and glass. In both cases, rods of length 10   cm were assumed that were heated to the uniform temperature T ₀. Obviously, the two materials differ in their property to conduct heat, i.e., in their thermal diffusivity α (see Eq. 6.82). For stainless steel we find λ  =   14   W   m^{−   1} K^{−   1}, a density of around ρ  =   7.9   g   cm^{−   3}, and a specific heat capacity of c_p   =   0.49   J   g^{−   1} K^{−   1} from Tab. 6.8. In this case we find the thermal diffusivity α of 3.62   mm²  s^{−   1} according to Eq. 6.82. For glass we find λ  =   0.095   W   m^{−   1} K^{−   1}, a density of around ρ  =   2.5   g   cm^{−   3}, and a specific heat capacity of c_p   =   0.67   J   g^{−   1} K^{−   1} from Tab. 6.8. In this case we find the thermal diffusivity α of 0.057   mm²  s^{−   1}.

Fig. 8.6. Heat conduction in two rods made of stainless steel (a) and glass (b). The rods are both 10   cm long and were at the temperature T₀ at the start of the experiment. At t₀  =   0 they were clamped between two fixtures both of which are at the temperature 0. The temperature profiles for (a) are at 0   s, 1   s, 10   s, 100   s, and 500   s. For (b) the profiles are drawn for 0   s, 10   s, 100   s, 1000   s, and 10 000   s. As the glass rod conducts heat significantly worse than stainless steel it will keep its temperature longer.

As can be seen, the temperature profile for the stainless steel rod is (more or less) uniform at the start of the experiment. The artifacts caused by the Fourier series expansion of the constant are still obvious. Due to the high thermal diffusivity, the temperature profile of the stainless steel rod drops quickly whereas the glass rod hardly changes temperature even after very long-time periods.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781455731411500083

Results of computational experiments

Dmitriy Tarkhov , Alexander Vasilyev , in Semi-Empirical Neural Network Modeling and Digital Twins Development, 2020

4.4.3 Problem for heat conduction equation with time reversal

This section discusses the classic incorrect problem, set in Section 1.4.3 . The problem is to restore the initial temperature distribution of a rod of unit length from the temperature distribution, measured at a specified time interval T, which was chosen equal to 1. The function sin(πx)   exp(−  π ² t) was chosen as the model solution. A stable approximation in such a formulation cannot be obtained without invoking additional information about the solution. We present the results of two ways to attract such information.

The first method consists of adding one or several terms corresponding to the initial condition to the functional J ₃. The computational experiment showed that it is enough to use one value, i.e., it is assumed that we know the initial condition at one point.

For network training, we used the method of growing networks with the addition of neurons one at a time and rejection of the unsuccessfully added one (an error after learning the entire network with an added neuron on the test set is more than an error without it). As a nonlinear optimization algorithm, a combination of the cloud method and RProp was used (see Chapter 3).

Let us give some results of calculations for the case when the number of test points at which the differential operator is calculated N  =   200, the number of test points in time at which the boundary conditions are checked N _b   =   50, the number of points at which the condition is checked at the final time point N _d   =   50. In the first experiment, the number of attempts to add a neuron 10, the number of neurons is 11. In this case, the maximum error was 8%.

In the second experiment, the number of attempts to add a neuron is 20, the number of neurons is 21. In this case, the maximum error was 2.5%.

In the third experiment, the number of attempts to add a neuron is 50, the number of neurons is 48. In this case, the maximum error was 1%.

The second method, instead of a known point in the initial condition, applies one or several random points inside the region (notable measurements at intermediate times). Below five such points were used, the number of attempts to add a neuron is 1559, the number of neurons is 51. In this case, the maximum error was 1%.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128156513000043

Thermal Effects in Friction Brakes

Andrew Day , in Braking of Road Vehicles, 2014

Brake Temperature, Stress and Deformation Prediction Using Computational Methods

Brake thermal analysis using computational methods such as the FE method enables the prediction of bulk component temperature, stress, and deformation profiles and distributions to a much higher accuracy and realism than by using analytical methods. The advantages of FE analysis include the realistic modelling of component geometry, material properties and boundary conditions. Even so, brake friction interface contact and pressure distribution and frictional heat generation remain complex thermophysical phenomena and simplifications are still required in the modelling. Thermomechanical properties and boundary conditions must be accurately specified, and results must be validated and verified. FE analysis can also be used to model 'macroscopic' effects resulting from variation in contact and pressure across the friction interface, arising from local deformation and distortion of the rubbing surfaces (see Chapter 2). Such deformation and distortion is usually transient in nature, and is associated with the modelling of thermoelastic instability (TEI). Since the work done and therefore the frictional power generated at any point in the brake friction interface is proportional to the local friction coefficient (μ), contact pressure (p) and linear sliding speed (v) as indicated in Equation (7.21), as a result of local μ variation with temperature and other factors, interface pressure variations with position and time, and sliding speed variation in the radial direction across the rubbing path, different regions of the friction surfaces will generate heat at different rates. In the regions where the heat generation is higher, the temperatures will be higher and these regions will thus experience more thermal expansion, causing an increase in local pressure and a further increase in local heat generation. The pressure distribution will therefore change across the friction surface, further increasing pressure in high-pressure areas and reducing interface pressure in other areas, and the increased local heat generation may result in the formation of hot bands and/or hot spots. However, the increased temperature in localised areas will eventually cause a reduction of local friction coefficient and increased wear of the friction material (see Chapter 1). The combined effect will be a reduction of local heat generation and, consequently, interface pressure in those areas, transferring contact and pressure increases to previously less loaded areas. Ultimately the hot areas may substantially reduce local interface pressures and even lose contact in a cyclic process that is broadly known as 'thermoelastic instability' (Barber, 1969).

(7.21) $q=\mu pv$

In hot spots and hot bands, friction surface temperatures can be much higher than those predicted using any assumption of uniform heat flux distribution across the friction interface. These high temperatures can cause significant changes in the local material properties at the friction surfaces and can also affect bulk material properties such as elastic and plastic modulus, Poisson's ratio, thermal conductivity, specific heat, thermal expansion coefficient, density, wear rate and friction coefficient, and cause degradation of the friction material in the friction surface (Day, 1983; Koetniyom et al., 2000; Tirovic and Sarwar, 2004), all of which must be adequately represented in FE modelling.

Thermal modelling of the friction pair enables the problem of heat partition to be avoided. More recent research has shown that the previously mentioned assumption that the friction surface temperatures are the same for the rotor and the stator is not realistic for a brake friction pair. Experimental evidence indicates that the friction surface temperature of a resin-bonded composite friction material is higher than the friction surface temperature of a cast-iron rotor, and there is a 'temperature jump' across the friction interface. This is discussed later in the context of predictive simulations involving thermal contact resistance at the interface. The 'five-phase' model of a brake friction pair (resin-bonded composite friction material and cast iron; see Figure 2.2) indicates how a temperature jump can occur. The surface of the friction material comprises a 'char' layer of friction material that has experienced high temperature arising from frictional energy transformation at the friction interface. The surface of the rotor is coated by a transfer film that is deposited during the process of frictional sliding, from the friction material. In between these two surfaces is wear debris. The physical mechanism of frictional energy transformation is that heat is generated by work done as the friction material slides over the rotor surface. As indicated in Chapter 2, the three main mechanisms of friction include adhesion, abrasion and deformation, and since the friction material is the weaker of the friction pair, it can be assumed that the majority of the work is done by adhesion, abrasion and deformation of the friction material at or near the friction surface, and the frictional heat energy is therefore generated in the surface layers of the friction material. The heat flows back into the friction material and across the friction interface into the rotor, and because the surface layers and the interfacial wear debris create a contact resistance to heat transfer, a temperature difference or 'jump' is created. The partition of heat between the friction material and the rotor is therefore determined by the thermal properties of the materials and the boundary conditions, but in practice for FE-based brake thermal analysis, realistic results have been obtained by assuming the heat is generated at the surface nodes of the friction material and transferred to the rotor through an 'interface contact conductance' (i.e. the reciprocal of contact resistance) of typically 2 kW/m² K.

Frictional contact at the microscopic scale is fundamental to the study of friction and wear and relates to the basic tribological characteristics of the selected friction pair. Much research is being undertaken to better understand and model these 'micro' effects, which are thus a research topic and not considered further here.

FE analysis starts with the generation of a 'mesh' of 'elements', which in modern FE analysis software systems is generated automatically. The type of element used must follow the FE system supplier's recommendations and verification examples should be completed before predictions are started. The element size and refinement of the mesh must be appropriate to the geometry of the part(s) being modelled, and their thermophysical properties, and the quality of the mesh must always be checked. For example, because the friction material has a low thermal conductivity, smaller elements (in the predominant direction of heat flow) must be specified to accommodate the steeper temperature gradient. This can be checked using the 'Fourier number' (Fo):

(7.22) $Fo=\alpha \mathrm{\Delta }t/\mathrm{\Delta }{x}^2$

where α is the thermal diffusivity of the friction material, Δt is the FE solution time increment and Δx is the element size thickness in the predominant direction of heat flow, usually normal to the friction interface. It is recommended that in FE mesh designs for heat transfer analysis and transient temperature prediction, Fo ≤ 0.5 for a 1D heat transfer model and Fo ≤ 0.17 for 3D models. As a general rule, for metal rotors the element size in the proximity of the friction surface should be approximately 1 mm thick, which with a time increment of 0.1 s usually gives good results. Setting the mesh design and the associated time increment values is much more difficult for the friction material, e.g. for FE modelling the complete friction pair. There is a large difference in the diffusivity (α) values between the rotor (disc/drum) material and the friction material (pad/lining). Because the solution time increment must be the same for the whole model, large differences in the preferred element size (Δx) can be required, and it is usually necessary to work to a satisfactory compromise.

Correct and meaningful material thermophysical properties are essential for successful and efficient FE modelling. The starting point is to assume constant linear, isotropic elastic properties for all FE thermal modelling, and then moving to include more advanced features such as temperature and time dependence, non-linearity and anisotropy as the detail of the exercise requires.

Boundary conditions in the form of surface heat transfer coefficients must be specified on all free surfaces. The heat energy input is specified (as calculated using the methods explained above). It is usual to specify an initial temperature distribution in the components that is the same as the surrounding or ambient temperature. Two types of temperature prediction analysis can be generated:

•

Steady-state temperature prediction, which calculates the temperature distribution when equilibrium has been reached.

•

Transient temperature prediction, which calculates the temperature distributions at different times as the heat flows through the components.

The most important is transient thermal analysis because single and repeated brake applications do not normally reach equilibrium conditions. Steady-state analysis is useful for drag braking where thermal equilibrium is reached. 2D axisymmetric or 3D approaches to FE modelling can be adopted. 2D axisymmetric FE modelling presents a mesh of a radial section through the disc, which is mathematically 'swept' through 360°; thus, any circumferential detail or variation such as the vanes in a ventilated disc or localised contact or heat flux cannot be included and so this approach is generally only useful for basic thermomechanical predictions. In 3D FE modelling, a sector of a ventilated brake disc bounded by planes of symmetry may be modelled to reduce the computational time requirement. In the past, computer power and availability were limited, and 2D axisymmetric modelling was preferred, but modern systems are sufficiently powerful that 3D FE modelling has become the standard.

The method for the prediction of brake disc temperatures, stresses and deformations using the FE method was demonstrated by Day et al. (1991). For a single brake application on a solid disc rotor of a passenger car front brake (vehicle and brake data are summarised in Table 7.1), the calculated maximum heat flux over the whole rubbing path (both sides of the disc) is approximately 1.5 MW/m² and the SSTR is 171°C. This is a simple design that was quickly modelled using the 2D axisymmetric FE mesh design illustrated in Figure 7.6 (Day et al., 1991). Detailed geometric features such as the 'undercut' where the friction ring joins the hat section can usually be ignored for the purposes of temperature prediction, being only required for stress analysis where they are known to be important, e.g. in areas of high stress concentration. However, it is usually convenient to use the same FE mesh for thermal and stress analyses as the computational time costs are relatively insignificant, so in practice detailed FE models as illustrated in Figure 7.7 are now used for these purposes.

Table 7.1. Vehicle and Brake Data

Vehicle Characteristics
Vehicle mass	1300 kg
Brake force distribution X 1/X 2	66:34
Brake Application Characteristics
Initial vehicle speed	100 km/h (28 m/s)
Final vehicle speed	0
Duration of brake application	4 s
Average deceleration	7 m/s2 (0.7 g)
Initial brake temperature	20°C
Front Brake Discs (Solid)
Disc outer diameter	227 mm
Disc ring inner diameter	127 mm
Disc thickness	11 mm
Total disc mass	3 kg
Disc ring mass	2.2 kg

Figure 7.6. 2D Axisymmetric FE Brake Disc Model (Day et al., 1991).

Figure 7.7. (a) 3D Sector Model of a Brake Disc Showing the Predicted Transient Temperature Distribution for a 0.4 g Brake Application. (b) The Same Transient Temperature Distribution Showing the Coning Distortion of the Friction Ring (Tang, private communication).

The frictional heat input for this brake was modelled assuming 98.4% of the total braking energy enters the disc. The braking power profile with time for a single brake application was modelled as in Figure 7.2(b), from which the local heat flux was calculated assuming uniform interface pressure (and thus frictional heat generation) on the friction surfaces of the disc. Surface heat transfer was neglected for a single brake application. Starting from an initial uniform temperature (the ambient temperature 20°C), the peak temperature predicted on the disc surface (rubbing path) was about 215°C after approximately 3 s of the 4 s brake application, which compares with 127°C calculated using Equation (7.19). The reason for the difference is that FE analysis predicted a significant radial temperature gradient from 215°C near the outer periphery to around 80°C at the inner radius of the rubbing path, while the analytical method assumes a uniform distribution of heat. High thermal gradients from the friction surfaces into the thickness of the disc were predicted by the FE analysis, and entirely disappeared towards the end of the brake application. Negligible temperature rise in the hat section of the disc indicates that the assumption that only the disc ring mass need be considered for the SSTR calculation was justified. At the end of the brake application the temperature predicted by the FEA was approximately 200°C compared with the calculated SSTR of 137°C.

Once the temperature distribution (steady-state or transient) at any point in time has been predicted, the FE analysis can be continued to predict the thermal stresses, deformations and displacements of a brake disc or drum in use. For this it can be assumed that the disc or drum flange (which is bolted to the hub) is a rigid boundary, i.e. there is no displacement of any node on this flange. Usually no mechanical loading is included in thermal deformation predictions for brake discs (clamp/actuation, friction or centrifugal forces) because the thermal stresses are the most significant, but for brake drums it is wise to check the effect of centrifugal forces. In all thermal stress and deformation analyses the thermal analyses are performed first and the predicted temperature data are usually stored for subsequent thermal stress and deformation analyses at selected braking times.

Example results from the FE analysis of a ventilated brake disc (from the front axle of a medium-sized passenger car) are shown in Figure 7.7, based on a 3D FE modelling procedure that utilises a 'sector' model, with symmetrical boundaries on radial planes either side of a complete vane region. For thermal analysis the assumption is that there is no heat flow across the boundary planes, and for stress analysis a plane strain boundary condition applies. Figure 7.7(a) shows the FE model with the predicted transient temperature distribution halfway through a single 0.4 g brake application, and Figure 7.7(b) illustrates the same temperature distribution superimposed on the deflected shape, i.e. the shape of the thermal deformation predicted for this temperature distribution. The detail of the junction between the friction ring joining the hat section is modelled in detail, and 'coning' distortion of the friction ring is clearly visible. These results illustrate the advantages of computational modelling methods like the FE method over classical analytical methods.

Although this analysis uses a 3D FE sector model, the thermal loading and boundary conditions assume circumferential symmetry, so full 3D modelling, which includes circumferential variation of contact and pressure at the friction interface, frictional heat generation, temperature distribution, thermal stresses and displacement, is needed to provide complete realism in FE-based thermal and thermomechanical predictions for the competitive design of high-duty brake rotors. There are many examples of this type of sophisticated modelling in published research; the frictional power can be modelled as being generated at the pad friction surface so that heat flows across the friction interface (via the interface contact resistance/conductance) to the disc, and as explained previously this avoids the need for heat partition. Because the heat flow into the disc varies circumferentially as it rotates, the disc surface temperature also varies circumferentially, as illustrated in Figure 7.8 (Tirovic, 2013).

Figure 7.8. Circumferential Variation in the Surface Temperature of a Rotating Brake Disc (Tirovic, 2013).

Results from thermal stress prediction by FE analysis are often quoted in terms of von Mises stress values. Although the von Mises criterion is not representative for cast-iron materials, the resultant stress values can give a broad indication of where temperature gradients may cause high thermal stresses and thermal cracking on the surface of the brake disc. Radial cracks or crazing on the rubbing path of a brake disc is termed thermal cracking (see Figure 7.9) and is a brake disc failure mode resulting from stresses created by elastoplastic deformation of the cast iron during the heating cycle, which in turn generates residual tensile stresses during cooling. Even if thermal cracking is avoided, experience has indicated that hot and cold judder can result from high temperature and high thermal stress gradients normal to the disc friction surface (see Chapter 10).

Figure 7.9. Disc Thermal Cracking.

Thermal deformation and distortion of brake discs is undesirable for many reasons. Coning can cause the pad/disc contact and pressure distribution to be biased towards the outer radius on the outboard side and the inner radius on the inboard side, which causes increased pedal travel, taper wear and non-uniform heat flux input into the disc. This leads to an increased risk of the formation of hot bands on the disc and hot spots or 'blue spots' associated with localised surface overheating and the localised change in phase of the cast-iron disc (or drum) material from pearlite to martensite. The overall effect is an increased incidence of thermal cracking and hot judder (Chapter 10). Good disc brake design should therefore aim to keep the braking energy input as even as possible, radially and circumferentially, so particular attention should be paid to the following factors:

•

Design the disc for minimal coning.

•

Design the disc and hub assembly for minimal runout and bolt-up distortion.

•

Use a good quality and consistent disc material.

•

Design the calipers and their mounting to be as rigid as possible.

As previously mentioned, high temperature gradients are a major cause of thermal cracking in brake discs. In a ventilated brake disc, in order to allow airflow through the vanes, only one cheek is connected to the hat region, which can therefore be either the inner or the outer cheek. When the disc is designed so the outboard face is connected to the hat section (as in Figure 7.7) it has been observed that cracking occurs more on that (outboard) face. Conversely, if the inner disc cheek is attached to the hat section, then cracks occur on the inboard face. This may be the result of the additional heat transfer from the outer cheek to the hub causing higher temperature gradients through its thickness. To minimise coning distortion the hat section and the undercut should be designed to be as flexible as possible, consistent with structural integrity and the casting process, and since this also has the effect of restricting heat flow to the hub then disc cracking propensity may also be reduced. When the inner disc cheek is attached to the hat section this is known as 'reverse ventilated' disc design, which provides a longer hat section with consequent increased flexibility and reduced coning distortion of the disc. However, cooling of the disc can be noticeably worse because the airflow path is more restricted between the wheel and the disc.

Brake disc materials have recently seen a trend towards higher carbon contents in order to maximise rotor conductivity at the expense of tensile strength, and this demands even more sophisticated methods of design prediction because of the safety-critical nature of the brake disc. The FE analyses shown here assumed linear elastic properties, but more advanced non-linear elastic–plastic thermomechanical FE models can be used to predict the thermal performance more accurately, e.g. thermal cracking under extreme braking conditions. Grey cast iron is much weaker in tension than in compression, which needs to be included in FE analysis to take account of the cyclic loading behaviour of the material. Koetniyom et al. (2000) generated a model based on uniaxial cyclic tests, which was used in a non-linear 3D FE thermomechanical analysis of two variants of a ventilated brake disc, conventional and reverse ventilated. The results indicated which areas of the disc might experience high temperatures and plastic stresses that might lead to thermal cracking under repeated thermal cycling.

For a full FE simulation and prediction of temperatures and thermal stresses in brakes, both the rotor and the stator should be included, together with some consideration of thermoelastic instability (TEI). The maximum temperatures where the frictional contact is concentrated at 'hot spots' on the rotor surface can be expected to be much higher, perhaps by an order of magnitude, depending on many factors. Such advanced modelling is well established but is beyond the scope of this book. TEI modelling will inevitably predict increased stresses at the surface of a brake disc or drum as a result of hot spotting, even causing local material yielding. A brake disc has frictional heat flux input from both sides and away from the friction surface; the effects tend to 'even out' so that TEI has relatively little influence on the stresses in the friction disc ring–hat transition area where the detail design is critically important. Therefore, 2D circumferentially symmetric disc models can still be very useful and quick for the optimisation of this region, before more complex analyses are performed.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123973146000073

THERMAL NONEQUILIBRIUM FORCED CONVECTION IN POROUS MEDIA

A.V. KUZNETSOV , in Transport Phenomena in Porous Media, 1998

62. Optimization of the initial temperature distribution - rate of "forgetting" the initial conditions

The second problem we have considered is that of a porous packed bed which is cooled by the flow of an incompressible fluid. The initial temperature of the packed bed is considered as the control and the amount of heat transferred to the fluid phase is utilized as the optimization criterion. It is necessary to maximize this quantity under the following constraints: (a) a given amount of heat is initially stored in the packed bed, (b) there is a given duration of the process, t_f, and (c) the initial temperature is subject to a variation between a given minimum value, u_min, and a given maximum value, u_max.

Again, a one-dimensional porous slab of the length L′ is considered and the dimensionless length of the slab is then defined by equation (71) and it is assumed that the inlet fluid phase temperature is constant during the process. In this case it is always possible to select the reference temperatures, T₁ and T₂, in equations (63) and (64), so that the dimensionless inlet temperature is zero. This essentially simplifies equations (69) and (70), because in this case the first term on the right-hand side of equation (69) and the second term on the right-hand side of equation (70) are zero.

The mathematical formulation of this problem is as follows. It is necessary to determine the optimal initial temperature distribution ${\hat{\theta }}_0\left(z\right)$ that maximizes the following functional

(82) $\Phi \left({\theta }_0\right)=\underset{\Lambda \text{L}}{\overset{t_f}{\int }}\varphi \left(\text{L},\text{t}\right)\text{dt}\to max$

where the function ϕ(L,t) is determined by equation (70), under the following constraints:

(83) $\underset{0}{\overset{\text{L}}{\int }}{\theta }_0(\xi )\text{d}\xi =\text{E}=\text{constant}$

and

(84) ${\text{u}}_{min}\le {\theta }_0\left(\text{z}\right)\le {\text{u}}_{max}$

The lower integration limit in equation (82) corresponds to the instant of time when the fluid front reaches the outlet boundary, and the upper limit corresponds to the duration of cooling. We consider the case when t_f > ΛL.

The problem given by equations (82)–(84) was rearranged in Kuznetsov [41] in a similar manner to that presented in Section 6.1 and it was transformed to an optimal control problem. It was then solved in Kuznetsov [41] by the minimum principle of Pontryagin which leads to the following requirement

(85) ${\hat{\theta }}_0\left(\text{z}\right)\left[{\lambda }_1-\Psi \left(\text{z}\right)\right]\to min$

where

(86) $\Psi \left(\text{z}\right)=exp\left(\Lambda \text{L}+\text{z}-\text{L}\right)\underset{\Lambda \text{L}}{\overset{{\text{t}}_{\text{f}}}{\int }}exp\left(-\text{t}\right){\text{I}}_0\left[{\left\{4\left(\text{L}-\text{z}\right)\left(\text{t}-\Lambda \text{L}\right)\right\}}^{1/2}\right]\text{dt}$

and λ₁ is the Lagrange multiplier.

Equation (85), when applied to the constraint (84), makes it possible to determine the optimal control, ${\hat{\theta }}_0\left(z\right)$ , as

(87) $\begin{array}{l}{\hat{\theta }}_0\left(\text{z}\right)={\text{u}}_{min}\text{if}{\lambda }_1-\Psi \left(\text{z}\right)>0\hfill \\ {\hat{\theta }}_0\left(\text{z}\right)={\text{u}}_{max}\text{if}{\lambda }_1-\Psi \left(\text{z}\right)<0\hfill \end{array}$

which requires the calculation of the value of the Lagrange multiplier, λ₁. To do this, the transcendental equation (83) needs to be solved taking into account equation (87) and to solve this problem we first took a segment that unequivocally contains the desired value of λ₁ Then an algorithm for finding a root of a transcendental equation on a given segment was applied to equation (83).

It is interesting to compare the value of the performance functional Φ(θ₀) at the optimal functions determined by equations (87) and (83) and at the functions ${\hat{\theta }}_0^{*}\left(\text{z}\right)\equiv \text{E}/\text{L}=\text{constant}\text{.}$ . These functions, ${\hat{\theta }}_0^{*}\left(\text{z}\right)$ , correspond to a constant uniform initial temperature distribution in the slab and clearly satisfy constraint (83). The ratio $\left[\Phi \left({\hat{\theta }}_0\right)-\Phi \left({\theta }_0^{*}\right)\right]/\Phi \left({\theta }_0^{*}\right)$ characterizes the gain in the amount of the heat energy transferred to the fluid phase when the optimal initial temperature distribution is utilized, instead of the constant uniform initial temperature. The dependence of this ratio on the duration of the process for u_min = 0, u_max = 1, E = L/2, Λ = 0.05 is depicted in Figure 3. It can be seen that the optimal initial temperature distribution makes it possible to considerably improve the discharging characteristics of the packed bed (up to approximately 25% for L=1) for short durations of the process. The influence of the initial temperature on the discharging efficiency quickly decreases with an increase of the duration and for t_f > 10 it becomes so small that it is hardly visible in Figure 3. This conclusion agrees with the well known property of heat transfer processes to forget initial conditions, see Prigogine [42]. According to Figure 3, the longer is the duration of the process, the smaller the amount of heat transferred to the fluid phase by t = t_f is influenced by the choice of the initial temperature distribution in the packed bed. For t_f >10 this quantity, which can be interpreted as a characteristic of a final state of the process, is practically independent on whether the optimal or the uniform temperature distribution is utilized as the initial condition. The rate at which the curves depicted in Figure 3 tends to zero characterizes the rate at which the process forgets its initial condition.

Figure 3. The time dependence on the gain in the amount of heat energy transfered to the fluid phase when the optimal initial distribution is utilized instead of the constant initial temperature of the duration of the process

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780080428437500052

Distribution of the temperature field in a pavement structure

Lijun Sun , in Structural Behavior of Asphalt Pavements, 2016

2.10.2 Feature of Pavement Temperature Field Under Continuously High Air Temperature

The longer the period of continuously high air temperature, the time span of the zigzag line representing the pavement temperature rising is greater and the difference between real measured pavement temperature and predicted pavement temperature is greater as well. In this section, the features of the pavement temperature field under a continuously high air temperature will be analyzed under the same and repeated influence of environmental factors. At the same time, the influence of some main factors such as a heat convection coefficient, thermal diffusivity coefficient, air temperature, solar radiation, and initial temperature distribution will be also analyzed.

In the short term, such as 6–7 days, the theoretical value of solar radiation should be similar. When the weather condition is sunny, the radiation energy absorbed by pavement should be similar, too; therefore, during the analysis phase, the daily value of solar radiation is regarded as constant.

Normally, when the pavement temperature is continuously rising, the air temperature should be also increased. The detailed increasing difference will be determined by the flow of air and distribution of moisture, etc., and so it is hard to carry out quantitative analysis. Therefore, the daily variations of air temperatures during the analysis phase are also regarded as the same (this situation has been observed).

The air temperature, solar radiation, and initial temperature distribution used in analysis are coming from the collected data from Wuhan at Aug. 11, 2005, as illustrated in Fig. 2.57.

Fig. 2.57. Climatic data in Wuhan, Aug. 11, 2005.

The initial temperature distribution at a lower position will be predicted by the earth temperature, and the initial temperature at the pavement surface will be determined by regression equations. The initial time for analysis is 0 O'clock, when the correlation of air temperature and pavement temperature is highest and the prediction error is lowest.

The parameters used in analysis are: initial temperature is expressed as 33.54-3x, the thickness for asphalt layer is 20   cm, the thermal conductivity is 0.0023   m²/h, and the thermal diffusivity is 1.6   W/mh. The heat convection coefficient of pavement surface is 10, the thickness of base layer is 40   cm, and the thickness for subgrade is 3   m.

The predicted temperature and measured temperature are plotted in Fig. 2.58. Because Aug. is the hottest month in Wuhan, then the maximal pavement temperature in 115   mm depth could be 48°C. From Fig. 2.58E, it can be seen that the predicted value fits the measured temperature well.

Fig. 2.58. Measured and predicted pavement temperature.

Under the application of repeated the climatic data as above, the variation of pavement temperature are shown in Fig. 2.59.

Fig. 2.59. Variation of pavement temperature under continuous heat flux.

Under the continuous heat flux, the pavement temperatures at difference depths are all increasing. But the value of temperature increasing are decreasing with time. After several cycles, the daily pavement temperature variations are becoming stable. The pavement temperature at the surface is the first to reach stability, and after 17 cycles, the temperature increase is only 1.9°C, of which the increase after the first cycle is around 0.8°C. After five cycles, the variations of pavement temperature at 5   cm depth are becoming stable. After three cycles, the maximal pavement temperature at 10   cm depth will be greater than 50°C; at this time, the pavement temperature at the surface is only slightly greater than 50°C. Thus, the pavement temperature within the 10   cm depth is greater than 50°C in a short period of time.

As illustrated in Fig. 2.60, the variation of temperature differences between the surface and the bottom of the asphalt layer under continuous heat flux are plotted. It can be determined that this temperature difference is decreasing periodically. Compared with the temperature at the pavement surface, the temperature of the internal pavement is relatively increasing. During a high-temperature period, the temperature difference (between 0 and 20   cm pavement depth) is decreasing from 25°C to 20°C. During the low-temperature period, the temperature difference is decreasing from −   5°C to −   10°C.The time when a pavement temperature difference is changing from positive value to negative value or from negative value to positive value is also postponed.

Fig. 2.60. Variation of temperature gradient in surface course and base course.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B978012849908500002X

THERMAL PROCESSES | Pasteurization

F.V.M. Silva , ... R. Simpson , in Encyclopedia of Food Microbiology (Second Edition), 2014

Heat Transfer Model for Pure Conduction

Heat transfer for pure conduction is based on Fourier's equation and can be written as follows:

[5] $\rho C_{\text{p}}\frac{\partial T}{\partial t}=\nabla k\nabla T$

If thermal conductivity (k) is independent of temperature and the food material is assumed to be isotropic, as is the case for most foods within the sterilization temperature range, then eqn [1] becomes

[6] $\frac{\partial T}{\partial t}=\alpha {\nabla }^{2}T$

The solutions for different geometries are not necessarily straightforward. In general, however, the dimensionless temperature ratio for constant retort temperature can be expressed as follows (Carslaw and Jaeger, 1959):

[7] $\frac{\mathrm{PT}-T_{\text{C}\text{.P}\text{.}}}{\mathrm{PT}-\mathrm{IT}}=f\left(\text{initial}\text{temperature}\text{distribition,}\text{geometry,}\text{thermal}\text{properties,}\text{time}\right)$

Therefore, if the initial temperature distribution, geometry, product (thermal properties), and time are maintained constant (only changing PT or IT), then the dimensionless temperature ratio of the solid must be the same at different PT or IT:

[8] $\frac{\mathrm{PT}-T_{\text{C}\text{.P}\text{.}}}{\mathrm{PT}-\mathrm{IT}}=\frac{\mathrm{P}{\mathrm{T}}^{\prime }-T_{\text{C}\text{.P}\text{.}}^{\prime }}{\mathrm{P}{\mathrm{T}}^{\prime }-\mathrm{I}{\mathrm{T}}^{\prime }}=\text{Constant}$

It is important to note that eqn [8] is valid for constant process temperatures (PT). A simplified analytical solution for homogeneous solids confined to a finite cylinder is presented in eqn [9] (Merson et al., 1978). This simplified solution is only valid for long periods of time (after the initial lag period when Fourier number > 6) and assumes a Biot number > 40 (the external heat resistance is negligible when compared to the internal resistance).

[9] $\frac{\mathrm{PT}-T_{\text{C}\text{.P}\text{.}}}{\mathrm{PT}-\mathrm{IT}}=2.0396\exp \left[-\left(\frac{\left({2.4048}^2\right)}{R^2}+\frac{{\mathrm{\pi }}^2}{1^2}\right)·\frac{k}{\rho C_{\text{p}}}·t\right]$

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123847300004043

Variational Methods in Nonconservative Phenomena

In Mathematics in Science and Engineering, 1989

5.9 The Temperature Distribution in a Finite Rod with a Nonzero Initial Temperature Distribution

In the previous sections, which were devoted to linear and nonlinear heat transfer problems, all examples were for semi-infinite, one-dimensional, insulated solids. In this and the next section, we give examples of heat transfer problems that do not fall into that class. The main characteristics of the problems discussed in the text that follows is that the notion of penetration depth cannot be used.

Consider the problem of finding the temperature field in a rod of length 2l. The rod is thermally insulated and its ends are kept at a constant temperature T = 0. The initial temperature distribution is parabolic (see Figure 5.9.1):

Fig. 5.9.1. Heat flow in an insulated, thermodynamically linear rod with the given initial temperature distribution.

(5.9.1) $T\left(0,x\right)=T_0\left[1-{\left[\frac{x}{l}\right]}^2\right].$

The problem can be stated as one of finding an approximate solution of the differential equation

(5.9.2) $\frac{\partial T}{\partial t}=a\frac{{\partial }^{2}T}{\partial x^2}$

subject to the initial condition (5.9.1) and the given boundary conditions

(5.9.3) $T\left(-l,t\right)=T\left(l,t\right)=0$

Let us consider the action integral

(5.9.4) $I={\displaystyle {\int }_{t_0}^{t_1}{\displaystyle {\int }_{-l}^l\left[\frac{\tau }{2}{\left(\frac{\partial T}{\partial t}\right)}^2-\frac{1}{2}a{\left(\frac{\partial T}{\partial x}\right)}^2\right]}}{e}^{t/\tau }dxdt,$

where the time interval (t ₀, t ₁) is fully specified. For the sake of simplicity, we assume a temperature distribution of the form

(5.9.5) $T=T_0\left[1-{\left(\frac{x}{l}\right)}^2\right]f\left(t\right),$

where f(t) plays the role of the generalized coordinate mentioned previously. In order to satisfy the initial condition (5.9.1), we take f(0) = 1.

Substituting (5.9.5) into the action integral (5.9.4) and integrating with respect to x, we obtain

(5.9.6) $I=\lambda {\displaystyle {\int }_{t_0}^{t_1}\left(\frac{2}{5}\tau l{\dot{f}}^2-\frac{a}{l}{f}^2\right)e^{t/\tau }dt\equiv \lambda {\displaystyle {\int }_{t_0}^{t_1}{L}_1\left(t,\tau ,f,\dot{f}\right)dt,}}$

where $\lambda =\frac{4}{3}{T}_0^2$ . The variational equation δI = 0 yields the Euler–Lagrangian equation:

$\frac{d}{dt}\left(\frac{\partial L_1}{\partial \dot{f}}\right)-\frac{\partial L_1}{\partial f}=0,$

which, in our case, after dividing by e ^t/τ , leads to

(5.9.7) $\frac{4}{5}\tau l\ddot{f}+\frac{4}{5}l\dot{f}+2a\frac{f}{l^2}=0.$

Letting τ → 0 and integrating the resulting equation of the first order, we find that

(5.9.8) $f=e^{-\left(5/2\right)\left(a/l^2\right)t}.$

Thus, the approximate solution is

(5.9.9) $T=T_0\left[1-{\left(\frac{x}{l}\right)}^2\right]e^{-5at/2l^2}.$

The exact solution of this problem is given in the form of a series (see [28], p. 98):

(5.9.10) $T=T_0\frac{32}{{\pi }^3}{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(2n+1\right)}^3}\cos \frac{\left(2n+1\right)x}{2l}}\exp \left[-\frac{a{\left(2n+1\right)}^2{\pi }^{2}t}{4l^2}\right].$

The graphical comparison between this exact expression and the approximate solution (5.9.9) is sketched in Figure 5.9.2 for various values of the parameter at/l ². This diagram reveals, in spite of the rather crude trial solution, a good agreement with the exact solution.

Fig. 5.9.2. The temperature distribution in a rod 0 &lt; x &lt; l. The numbers on the curves given in parentheses are the values of at/l ². The solid lines represent the exact solution. Small circles denote the approximate solution.

This simple problem has been treated by numerous authors. For example, in Reference [6], the same result was obtained by means of the local potential method. In Reference [21], it was shown that a modification of the variational principle with vanishing parameter, whose Lagrangian function is given by (5.3.11), generates the same result as obtained here. A two-component version of the same problem (i.e., simultaneous heat and mass transfer) has been discussed, starting from the Lagrangian function (5.3.15), in [22]. The same approximate solution is also reported there.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/S0076539208618056

Mathematical Models and Numerical Solutions for Thermal Storage Processes

Pei-Wen Li , Cho Lik Chan , in Thermal Energy Storage Analyses and Designs, 2017

4.2.4 Numerical Solution for the One-Dimensional Model for Dual-Media Sensible Thermal Storage

A number of analyses and solutions to the heat transfer governing equations of a working fluid flowing through a filler packed bed have been presented in the past (Schumann [17]; Shitzer & Levy [18]; McMahan [19]; Beasley and Clark [20]; Zarty & Juddaimi [21]). As the pioneering work, Schumann [17] presented a set of equations governing the energy conservation of fluid flow through porous media. Schumann's equations have been widely adopted in the analysis of thermocline heat storage utilizing solid filler material inside a tank. His analysis and solutions were for the special case where there is a fixed fluid temperature at the inlet to the storage system. In most solar thermal storage applications, this may not be the actual situation. To overcome this limitation, Shitzer and Levy [18 ] employed Duhamel's theorem on the basis of Schumann's solution to consider a transient inlet fluid temperature to the storage system. The analyses of Schumann and Shitzer and Levy, however, still carry with them some limitations. Their method does not consider a nonuniform initial temperature distribution. For a heat storage system, particularly in a solar thermal power plant, heat charge and discharge are cycled daily. The initial temperature field of a heat charge process is dictated by the most recently completed heat discharge process, and vice versa. Therefore, nonuniform and nonlinear temperature distribution is typical for both charge and discharge processes. To consider a nonuniform initial temperature distribution and varying fluid temperature at the inlet in a heat storage system, numerical methods have been deployed by researchers in the past.

To avoid the long mathematical analysis necessary in analytical solutions, numerical methods used to solve the Schumann equations were discussed in the literature by McMahan [19,22] and Pacheco et al. [23], and demonstrated in the TRNSYS software developed by Kolb and Hassani [24]. Based on the regular finite-difference method, McMahan gave both explicit and implicit discretized equations for the Schumann equations. Whereas the explicit solution method had serious solution stability issues, the implicit solution method encountered an additional computational overhead, thus requiring a dramatic amount of computation time. The solution for the complete power plant with thermocline storage provided by the TRNSYS model in the work of Kolb and Hassani [24] cites the short time step requirement for the differential equations of the thermocline as one major source of computer time consumption. To overcome the problems encountered in the explicit and implicit method, McMahan et al. [22] also proposed an infinite-NTU method. This model, however, is limited to the case in which the heat transfer of the fluid compared to the heat storage in fluid is extremely large.

The present study has approached the governing equations using a different numerical method [25]. The governing equations have been reduced to dimensionless forms, which allow for a universal application of the solution. The dimensionless hyperbolic type equations are solved numerically by the method of characteristics. This numerical method overcomes the numerical difficulties encountered in McMahan's work—explicit, implicit, and the restriction on infinite-NTU method (McMahan [19] and McMahan et al. [22]). The current model gives a direct solution to the discretized equations (with no iterative computation needed) and completely eliminates any computational overhead. A grid-independent solution is obtained at a small number of nodes. The method of characteristics and the present numerical solution has proven to be a fast, efficient, and accurate algorithm for the Schumann equations.

The nondimensional energy balance equations for the HTF and rocks can be solved numerically along the characteristics [26–28]. Eq. (4.8) can be reduced along the characteristic $t^{*}=z^{*}$ so that:

(4.53) $\frac{D{\theta }_f}{D{t}^{*}}=\frac{1}{{\tau }_r}\left({\theta }_s-{\theta }_f\right)$

Separating and integrating along the characteristic, the equation becomes:

(4.54) ${\displaystyle \int }d{\theta }_f={\displaystyle \int }\frac{1}{{\tau }_r}\left({\theta }_s-{\theta }_f\right)d{t}^{*}$

Similarly, Eq. (4.11) for the energy balance of rocks is reposed along the characteristic $z^{*}=$ constant so that:

(4.55) $\frac{d{\theta }_s}{d{t}^{*}}=-\frac{H_{CR}}{{\tau }_r}\left({\theta }_s-{\theta }_f\right)$

The solution for Eq. (4.55) is very similar to that for Eq. (4.53) but with the additional factor of H _CR . The term H _CR is simply a ratio of fluid heat capacitance to rock heat capacitance. Therefore, the equation for the solution of θ _s will react with a dampened speed when compared to θ _f , as the filler material must have the capacity to store the energy being delivered to it, or vice versa. Finally, separating and integrating along the characteristic for Eq. (4.55) results in:

(4.56) ${\displaystyle \int }d{\theta }_s=-{\displaystyle \int }\frac{H_{CR}}{{\tau }_r}\left({\theta }_s-{\theta }_f\right)d{t}^{*}$

There are now two characteristic equations bound to intersections of time and space. A discretized grid of points, laid over the time-space dimensions, will have nodes at these intersecting points. A diagram of these points in a matrix is shown in Fig. 4.5. In space, there are i  =   1, 2, …, M nodes broken up into step sizes of Δz* to span all of z*. Similarly, in time, there are j  =   1, 2, …, N nodes broken up into time-steps of Δt* to span all of t*. Looking at a grid of the ϑ nodes, a clear picture of the solution can arise. To demonstrate a calculation of the solution we can look at a specific point in time, along z*, where there are two points, ϑ _1,1 and ϑ _2,1. These two points are the starting points of their respective characteristic waves described by Eqs. (4.53) and (4.56). After the time Δt* there is a third point ϑ _2,2 which has been reached by both wave equations. Therefore, Eq. (4.54) can be integrated numerically as:

Fig. 4.5. Diagram of the solution matrix arising from the method of characteristics.

(4.57) ${\displaystyle {\int }_{{\vartheta }_{1,1}}^{{\vartheta }_{2,2}}d{\theta }_f}={\displaystyle {\int }_{{\vartheta }_{1,1}}^{{\vartheta }_{2,2}}\frac{1}{{\tau }_r}}\left({\theta }_s-{\theta }_f\right)d{t}^{*}$

The numerical integration of the right-hand side is performed via the trapezoidal rule and the solution is:

(4.58) ${\theta }_{f_{2,2}}-{\theta }_{f_{1,1}}=\frac{1}{{\tau }_r}\left(\frac{{\theta }_{s_{2,2}}+{\theta }_{s_{1,1}}}{2}-\frac{{\theta }_{f_{2,2}}+{\theta }_{f_{1,1}}}{2}\right)\mathrm{\Delta }{t}^{*}$

where ${\theta }_{f_{1,1}}$ is the value of θ _f at ϑ _1,1, and ${\theta }_{f_{2,2}}$ is the value of θ _f at ϑ _2,2, and similarly so for θ _r .

The integration for Eq. (4.36) along $z^{*}=\mathit{constant}$ constant is:

(4.59) ${\displaystyle {\int }_{{\vartheta }_{2,1}}^{{\vartheta }_{2,2}}d{\theta }_s}={\displaystyle {\int }_{{\vartheta }_{2,1}}^{{\vartheta }_{2,2}}-\frac{H_{CR}}{{\tau }_r}}\left({\theta }_s-{\theta }_f\right)d{t}^{*}$

The numerical integration of the right-hand side is also performed via the trapezoidal rule and the solution is:

(4.60) ${\theta }_{s_{2,2}}-{\theta }_{s_{2,1}}=-\frac{H_{CR}}{{\tau }_r}\left(\frac{{\theta }_{s_{2,2}}+{\theta }_{s_{2,1}}}{2}-\frac{{\theta }_{f_{2,2}}+{\theta }_{f_{2,1}}}{2}\right)\mathrm{\Delta }{t}^{*}$

Eqs. (4.58) and (4.60) can be reposed as a system of algebraic equations for two unknowns, ${\theta }_{f_{2,2}}$ and ${\theta }_{r_{2,2}}$ , while θ _f and θ _r at grid points ϑ _1,1 and ϑ _2,1 are known.

(4.61) $\left[\begin{array}{cc}\hfill 1+\frac{\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\hfill & \hfill -\frac{\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\hfill \\ \hfill -\frac{H_{CR}\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\hfill & \hfill 1+\frac{H_{CR}\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\theta }_{f_{2,2}}\hfill \\ \hfill {\theta }_{s_{2,2}}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {\theta }_{f_{1,1}}\left(1-\frac{\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\right)+{\theta }_{r_{1,1}}\frac{\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\hfill \\ \hfill {\theta }_{f_{2,1}}\frac{H_{CR}\mathrm{\Delta }{t}^{*}}{2{\tau }_r}+{\theta }_{r_{2,1}}\left(1-\frac{H_{CR}\mathrm{\Delta }{t}^{*}}{2{\tau }_r}\right)\hfill \end{array}\right]$

Cramer's rule [28] can be applied to obtain the solution efficiently. It is important to note that all coefficients/terms in Eq. (4.61) are independent of z*, t*, θ _f , and θ _s , thus they can be evaluated once for all. Therefore, the numerical computation takes a minimum of computing time, and is much more efficient than the method applied in references [19,22].

From the grid matrix in Fig. 4.5 it is seen that the temperatures of the rock and fluid at grids ϑ _i,1 are the initial conditions. The temperatures of the fluid and rock at grid ϑ _{1, 1} are the inlet conditions, which vary with time. The inlet temperature for the fluid versus time is given. The rock temperature (as a function of time) at the inlet can be easily obtained using Eq. (4.11), for which the inlet fluid temperature is known. Now, as the conditions at ϑ _1,1, ϑ _{1, 2}, and ϑ _{2, 1} are known, the temperatures of the rocks and fluid at ϑ _{2, 2} will be easily calculated from Eq. (4.61).

Extending the preceding sample calculation to all points in the ϑ grid of time and space will give the entire matrix of solutions in time and space for both the rocks and fluid. While the march of Δz* steps is limited to $z^{*}=1$ the march of time Δt* has no limitation.

The previous numerical integrations used the trapezoidal rule; the error of such an implementation is not straightforwardly analyzed but the formal accuracy is on the order of O(Δt*   2) for functions [28] such as those solved in this study.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128053447000043

The BEM for Time Dependent Problems

John T. Katsikadelis , in The Boundary Element Method for Engineers and Scientists (Second Edition), 2016

Examples

The FORTRAN code employed for the hyperbolic equation has been adjusted to solve the initial boundary value problem (9.30), (9.31). Numerical results are presented to illustrate the efficiency of the AEM.

Example 9.3

In this example the heat conduction equation which describes the unsteady temperature distribution in a two-dimensional homogeneous orthotropic body occupying the domain $\mathrm{\Omega }$ is studied. The temperature distribution is governed by the following initial value problem (cf. Eq. 8.117)

(9.32) $\rho (\mathbf{x})c\dot{u}=\nabla \cdot (\mathbf{D}\nabla u)+f(\mathbf{x},t)\mathrm{in}\mathrm{\Omega }$

(9.33a) $u=\overline{u}\mathrm{on}{\mathrm{\Gamma }}_1$

(9.33b) $\nabla u\cdot \mathbf{m}=-{\overline{q}}_m\mathrm{on}{\mathrm{\Gamma }}_2$

(9.33c) $\nabla u\cdot \mathbf{m}=a(u-u_a)\mathrm{on}{\mathrm{\Gamma }}_3$

and

(9.34) $u(\mathbf{x},0)=g_1(\mathbf{x})\mathrm{in}\mathrm{\Omega }$

where $\mathrm{\Gamma }={\mathrm{\Gamma }}_1\cup {\mathrm{\Gamma }}_2\cup {\mathrm{\Gamma }}_3$ ; $u(\mathbf{x},t)$ is the temperature in the solid; $\mathbf{D}$ is the conductivity matrix; $c$ the specific heat; $\rho$ the mass density; $f(\mathbf{x},t)$ the external heat flux; ${\overline{q}}_m$ the intensity of the heat input on ${\mathrm{\Gamma }}_2$ ; $\overline{u}$ the specified temperature on ${\mathrm{\Gamma }}_1$ ; $a$ the heat transfer coefficient and $u_a$ the ambient temperature. Finally, $g_1(\mathbf{x})$ is the initial temperature distribution; $\mathbf{m}=k_x\mathbf{i}+k_y\mathbf{j}$ the vector in the direction of the conormal with $k_x,k_y$ being the specified thermal conductivities.

For the problem at hand the conductivity matrix is

(9.35) $\mathbf{D}=\left[\begin{array}{cc}{k}_x& 0\\ 0& k_y\end{array}\right]$

and Eq. (9.32) reduces to

(9.36) $\rho c\dot{u}=k_{x}u{,}_{xx}+k_{y}u{,}_{yy}+f(\mathbf{x},t)\mathrm{in}\mathrm{\Omega }$

The problem has been solved for the rectangular body $\mathrm{\Omega }:\{a\times b\}$ with $f(\mathbf{x},t)=0$ ,

(9.37a) $u(0,y,t)=u(x,0,t)=u(a,y,t)=u(x,b,t)=0\mathrm{on}{\mathrm{\Gamma }}_1=\mathrm{\Gamma }$

and

(9.37b) $u(x,y,0)=u_0,\mathrm{in}\mathrm{\Omega }$

Equation (9.36) admits an analytic solution [9]

(9.38a) $u_{ex}(x,y,t)=\sum _{n=1}^{\infty }\sum _{j=1}^{\infty }{A}_n\sin \frac{n\pi x}{a}\sin \frac{j\pi y}{b}\exp \left[-\left(\frac{k_x{n}^2{\pi }^2}{a^2}+\frac{k_y{j}^2{\pi }^2}{b^2}\right)t\right]$

where

(9.38b) $A_n=\frac{4u_0}{nj{\pi }^2}[{(-1)}^n-1][{(-1)}^j-1]$

The adopted numerical values are: $a=b=3$ , $k_x=k_y=1.25$ $c=1$ , $\rho =1$ , and $u_0=30$ .

The solution was computed using $\mathrm{N}=200$ constant boundary elements and $\mathrm{L}=121$ domain nodal points distributed uniformly using MQs with shape parameter $c=0.2$ . The time history of the temperature at point (1.5,1.5) is shown in Fig. 9.9, while Fig. 9.10 shows the relative error $[u(1.5,1.5,t)-u_{ex}(1.5,1.5,t)]/u_{ex}(1.5,1.5,t)$ with $\mathrm{\Delta }t=0.005$ . Moreover, Fig. 9.11 shows the contours of the temperature distribution at $t=0.5$ . The semidiscretized parabolic equations were solved using the algorithm given in Table 9.3.

Figure 9.9. Time history of the temperature at point ( $1.5,1.5$ ) in Example 9.3.

Figure 9.10. Relative error at point $(1.5,1.5)$ in Example 9.3.

Figure 9.11. Contours of the temperature distribution at $\mathit{t}=0.5$ in Example 9.3.

Example 9.4

In this example the heat conduction in the inhomogeneous anisotropic plane body shown in Fig. 9.12 is studied. The conductivity matrix is given as

Figure 9.12. Geometry of the domain and boundary conditions in Example 9.4.

(9.39) $\mathbf{D}=\left[\begin{array}{cc}{(2x+y+2)}^2& {(x-y)}^2\\ {(x-y)}^2& {(x+2y+2)}^2\end{array}\right]$

The transient response is described by Eq. (9.32), which by virtue of Eq. (9.39) becomes

(9.40) $\begin{array}{c}\rho (\mathbf{x})c\dot{u}={(2x+y+2)}^2{u,}_{xx}+2{(x-y)}^2{u,}_{xy}+{(x+2y+2)}^2{u,}_{yy}\\ +(6x+6y+8){u,}_x+(6x+6y+8){u,}_y+f(x,y,t)\end{array}$

The considered boundary and initial conditions are (see Fig. 9.12)

(9.41a) $u=\overline{u}\mathrm{on}{\mathrm{\Gamma }}_u$

(9.41b) $\nabla u\cdot \mathbf{m}=-{\overline{q}}_m\mathrm{on}{\mathrm{\Gamma }}_q$

(9.42) $u(\mathbf{x},0)=g_1(\mathbf{x})\mathrm{in}\mathrm{\Omega }$

Equation (9.40) for

$c=1,\rho =1,f(x,y,t)=[17(x^2+y^2)+5xy-16]\exp (-t)$

$\overline{u}=(x^2+y^2-5xy)\exp (-t),$

${\overline{q}}_m=[(2x-5y){(2x+y+2)}^2+(2x-5y)(x-y)]\exp (-t),\mathrm{and}$

$g_1(x,y)=(x^2+y^2-5xy),$

admits the analytical solution

(9.43) $u_{exact}=(x^2+y^2-5xy)\exp (-t)$

The AEM solution was computed using $\mathrm{N}=200$ constant boundary elements and $\mathrm{L}=103$ domain nodal points distributed uniformly (Fig. 9.13) using MQs with shape parameter $c=0.1$ . The time history of the temperature at point (0.5,0.5) is shown in Fig. 9.14, while Fig. 9.15 shows the error $u(0.5,0.5,t)-u_{exact}(0.5,0.5,t)$ for $\mathrm{\Delta }t=0.005$ . Moreover, Fig. 9.16 shows the contours of the temperature distribution at $t=1$ .

Figure 9.13. Boundary and domain nodal points.

Figure 9.14. Time history of the temperature at point (0.5, 0.5) in Example 9.4.

Figure 9.15. Error $\mathit{u}-{\mathit{u}}_{\mathit{e}\mathit{x}\mathit{a}\mathit{c}\mathit{t}}$ at point (0.5, 0.5) in Example 9.4.

Figure 9.16. Contours of the temperature distribution at $\mathit{t}=1$ .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128044933000096

mccaffertysaink1971.blogspot.com

Source: https://www.sciencedirect.com/topics/engineering/initial-temperature-distribution

Share this post