Model Definition with Hysteresis - Nemlineáris hőterjedés és hőképanalízis

In this study a 2D transient heat transfer problem is assumed with temperature-dependent thermal diffusivity κ (it is assumed that the temperature dependence of the volumetric heat capacity c_v can be neglected and the expression of thermal diffusivity

κ = is valid, where λ is the thermal conductivity). The heat equation is written with a dimensionless temperature θ in the form

( ) ( )

( )

where the thermal diffusivity κ

( )

θ is a function of the temperature, the spatial coordinates x,y∈Ω, Ω=

[ ] [ ]

0,L × 0,L are held. The source term s

( )

θ t, is not detailed in this study. The initial condition of θ

(

x,y,0

)

=θ_init, x,y∈Ω and the first order boundary conditions θ

(

x,y t,

)

=θ_w, x,y∈∂Ω are assumed. The nonlinearity of the system originates from the composite material properties, since the diffusivity shows hysteresis like temperature-dependence. The thermal diffusivity as a function of the temperature is shown in Fig. 1 and corresponds to the measured data presented in [113].

Santos et al. [113] observed that during a heating cycle, in the range of lower temperatures, there is a sudden change in the thermal diffusivity, according to the microstructure modification, in the referred case, the phase transition of α →β cristobalite. Furthermore the thermal diffusivity upon heating remains relatively unaltered until a higher temperature range, where it presents another steep increase.

During cooling a reverse process takes place, but higher diffusivity values can be

expected. Furthermore, experiments showed that thermal cycling can reduce the hysteresis of diffusivity in composites.

In the following, I assumed that the temperature dependent diffusivity hysteresis could be modelled with a temperature dependent hysteresis operator [61], [69]. This assumption, i.e. that hysteresis of diffusivity has been invariably present, has been established by a series of steady state measurements in [113]. In this study an appropriate parameterized Preisach-type hysteresis model is introduced. Preisach-type of hysteresis models have been widely used to describe the ferromagnetic hysteresis [57].

0 0.2 0.4 0.6 0.8 1

25 30 35 40 45 50 55 60 65 70 75

θ κ ( 10-8 m2 /s)

cycle 2 cycle 1

cycle 3

Fig. 1. Hysteresis of thermal diffusivity with experimentally observed damping effects of thermal cycling [113]

The main assumption of the Preisach model is that the material can be modelled by a parallel summation of weighted elementary (bistable) hysteresis relays. Each relay is represented mathematically by a hysteresis operator γ_αβ , where α , β represent the state of relay ‘0’ and ‘1’. The Preisach hysteresis model can be applied as a general hysteresis function that has single-input single-output. When modelling the diffusivity hysteresis, the input is the temperature θ and the thermal diffusivity κ

( )

θ is the output.

In this macroscopic analysis, the physical interpretation of the elementary relay is unnecessary. A detailed microscopic explanation of the hysteresis phenomena in composites can be referred to in [19]. It can generally be considered that a relay corresponds to a homogeneous cluster of material, which has commonly changing physical properties, two states of an elementary bistability (cluster) corresponding to two different solid phases, at a low temperature stable phase and at a high temperature stable phase, respectively, and the hysteresis function represents the statistics on the probability of the state of those clusters. It has to be mentioned that in the referred literature [113] it is explained that there are some other microcrystal changes with changing temperature, for instance opening and closing of the internal flaws that could cause diffusivity variations as well.

The Preisach plane with coordinates of θ_α, θ_β are shown in Fig. 2a. Considering that each elementary hysteresis relay has a weight factor μ

(

θα,θβ

)

, by knowing the input history θ

( )

t , the output κ

( )

θ can be obtained by counting how many relays are in state

‘1’. Mathematically this counting process can be written as

( ) ( ) ( )

Fig. 2. a) The Preisach triangle with symbolic relays and b) the corresponding weight function interpreted over the discretized triangle

Repeated heating-cooling cycles can decrease the magnitude of the thermal diffusivity and decrease the area of the hysteresis loops as well. This can be introduced into the model through a cycle factor, ν . So, the continuous form of the modified Preisach operator is where ν is the heating-cooling cycle number, the bias thermal diffusivity κ^ν₀ depends on ν , b^ν is an appropriate multiplier that can vary with ν as well.

The numerical implementation of the Preisach scalar model is based on the Preisach triangle of elementary hysteresis relays. The normalized distribution function

(

^θ ^θ

)

¹ ₂

(

₁ ₂ ^σ₁² ^σ₂²

)

μ _α, _β = / N _D m ,m , , is discretized on a 500×500 grid, N₂_D represents the

variances σ₁², σ₂². The graphical representation of the applied distribution function can be seen in Fig. 2b. The diffusivity at grid points ^κ

( )

^θ^k⁺¹ is calculated by a recursive formula provided by discretization of the Preisach operator in (3)

( ) ( ) ( ) ( )

^, ⁱ ^j ^m^, ⁱ^,^j ^P

where T_P is the area of the domain over which in the last input step the hysteresis relays have changed their state and m is the number of relays in the domain P. The sign of the second part of the right-hand side depends on the direction of the temperature change.

The hysteresis cycles in Fig. 1 can be provided by the following parameters:

525 1.2.2 Discretization of the Nonlinear Model

The solution domain Ω has been discretized by a Cartesian mesh with a uniform mesh size h=L n. Time differentiation is discretized by the Crank-Nicholson method that is second-order accurate [121]. The right hand side of (1) is discretized in two steps. First the thermal flux differences are calculated, considering cell interfaces at the distance of

Depending on the time averaging of the diffusivity, various approaching methods can be obtained to (1). The Crank-Nicholson method with space averaged, time-instantaneous diffusivity has the form of

( ) ( ) ( )

where the first superscript k refers to the time step and the second superscript l refers to the iteration cycle. The weighting factor is σ =0.5 for the Crank-Nicholson method, J is the heat flux and ΔJ denotes the flux differences. Each flux difference has been computed in the same manner, for example

The thermal fluxes in the iteration step l+1 and at time step k+1 are determined by the diffusivity hysteresis model k⁺¹l,⁺¹⁼

(

ij^k⁺¹^l,⁺¹

)

ij Hθ

κ . At time step k, the diffusivity hysteresis ^κ_ij ⁼^H

( )

^θ_ij^k is not involved in the inner iteration marked by l.

Another approaching scheme can be evaluated using time-averaged temperatures to determine the diffusivity hysteresis. In this case the discretized form of (1) is

( ) ( ) ( )

calculation is shown in Fig. 3.

( )

⎟

Fig. 3. Semi-implicit finite difference method, diffusivity calculation with hysteresis operator at grid points of half time steps. Intergrid averaged diffusivities are used for the calculation of fluxes

Diffusivity at cell interfaces is interpolated linearly by

(

1¹² ¹ ¹² ¹

)

The thermal flux across the cell interfaces at time level k+1 and in iteration cycle l+1 is calculated according to

( )

^h fluxes could be derived in the same manner.

Substituting thermal fluxes (9) into (7) produces the so-called time-centred discretization scheme

( ) ( )

Assembling (10) for each interior node, it yields a system of equations

2 multiplied by Δt. Subscript h refers to the level of the spatial resolution.

The semi-implicit finite difference method (10) resembles the 2D variation of the DuFort-Frankel scheme, which is consistent with the original PDE only if

2 <<1

Δtκ_max h , where κmax=max

( )

κij , 1i,j=1,L,n− [121]. To avoid solutions that do not satisfy the physical nature of the modelled system, the upper limit to the time step can be evaluated from the following inequalityΔt≤h²

(

1−σ

)

κ_max

)

1.2.3 The Coarse Level Iterating Multigrid Solver

The nonlinearity of κ in method (10) needs iteration. The Newton iteration is very expensive in CPU time because of the large set of unknown variables [121]. Hysteresis-like temperature dependence of the diffusivity dramatically increases the computational cost by handling, at each grid point, the local memory in the discretized model (12).

Treatment of this memory is very resource and time consuming. The finer the resolution, the more hysteresis memory has to be treated. Using coarser resolution decreases the computational costs, but produces a coarser solution as well. In this study I suggest using a coarse level iterating multigrid algorithm, in which the fine resolution of the temperature field is preserved, but the hysteresis of diffusivity is reduced to a coarser grid, by utilising the multigrid algorithm that works on several resolutions [61], [65].

In a general case, the multigrid solver assumes an initial guess on the fine resolution for the diffusivity k ,

( )

^kh

h H θ

κ ⁺¹²⁰ = and for the temperature field θ^k_h⁺¹^,⁰ =θ^k_h. In every iteration cycle the diffusivity coefficients have to be re-calculated for all fine grid

points. To reduce the solution time, before the iteration, problem (12) should be reduced to the next coarser grid, with mesh size H =2h, in the form of

2 1 1

1 1 1 1

1 + + + + + +

+ ^l, ^k_H ^l, = ^k_H ^l, ^k_H+ ^k_H

H θ B θ S

A , (13)

where the initial temperature θ^k_H⁺^1,⁰ of time step k+1 has been determined by the restriction of the fine grid solution at time step k, θ^k_H⁺¹^,⁰ =θ^k_H =Rθ^k_h. The initial diffusion coefficient field has been defined with the hysteresis operator as κH⁰ =H

( )

θ^kH . The source term S^k_h⁺¹² has been restricted to S^k_H⁺¹², or S^k_H⁺¹² could be calculated directly as well. The temperature field approximation for time step k+1 and mesh size

H could be evaluated by iteration (13) with the new values of diffusivity

( )

(

¹ ¹⁰

)

1 2

1 ^l, 05 ^k_H ^l, ^k_H ^,

H⁺ ⁺ =H . θ ⁺ +θ ⁺

κ , until the absolute difference decreases below a

preliminarily defined ε limit, θ^k_H⁺¹^l,⁺¹−θ^k_H⁺¹^l, <ε. After this so-called inner-iteration, the diffusivity and the temperature are interpolated to the fine grid. By denoting P as the prolongation operator and I^h_H as the interpolation operator, the fine grid diffusivity is κ^k_h⁺¹² =S_GSPκ^k_H⁺¹²^l,⁺¹, the initial temperature is θ^k_h⁺¹^,⁰ =I^h_Hθ^k_H⁺¹^l,⁺¹. The prolonged diffusivity needs smoothing, S_GS denotes the Gauss-Seidel smoothing operator corresponding to ( A-14) in Appendix A. The interpolated temperature could be seen as a good initial solution to the grid level h. To evaluate the fine grid solution θ^k_h⁺¹, some multigrid V-cycles are applied, until the convergence criterion is reached. During the fine grid V-cycles the diffusivity field remains unaltered. The graphical scheme of the iteration method can be seen in Fig. 4.

Fig. 4. Schematic representation of the coarse level iterating multigrid method

This proposed coarse level iterating hierarchical algorithm consists of two types of multigrid methods. The inner-iteration is on the coarser grid, the solver is the so-called full multigrid algorithm. The advantages of this FMG method are the exact initial

the coarser grid. The stencils of A^k_H⁺¹^l,⁺¹ and B^k_H⁺¹^l,⁺¹ are changed point-by-point,

respectively, where β =Δt 2H². In the proposed algorithm the diffusivity fields κ, the temperature fields θ and the source terms S are stored at every grid level. The FMG determines the solutions at each grid level, so the approximated truncation error can be the stop criterion of the FMG cycles. The hysteresis memories are treated and stored on the coarser grid level.

To refine the temperature field after the inner iteration, some simple multigrid V-cycles run on the finest grid. The temperature field is iterated on the coarser grid because of nonlinearity, and after interpolation it is iterated on the fine grid to increase the accuracy. The diffusivity is iterated on the coarser grid as well, but after prolongation remains unchanged on the fine grid. The high level V-cycles (MG-V) coarsening technique and stencil determination are similar to the techniques used in FMG, but the temperature field is calculated only on the finest grid. On the coarse levels MG-V works with the corrections of the temperatures. The V-cycles iteration could be stopped when the residual norm has been reduced below the η times of the initial norm, i.e.

solved in both multigrid algorithms by the successive over-relaxations method (SOR).

The schematic flowchart of the proposed coarse level iterating multigrid algorithm is shown in Fig. 5.

− θ^k_H⁺¹^,^l⁺¹ ^k_H⁺¹^,^l

1 k h

θ+ k

θh

reduction norm

Fig. 5. Flow Chart of the Proposed Coarse Level Iterating Multigrid Algorithm

1.2.4 Numerical Results

In this section, first the proposed method and algorithm are analysed, then through two test problems, some characteristics of the diffusion hysteresis are presented. In all the test problems the source term is considered to be zero. The test parameters are

m 5 0.

L= , κ₀ =2.5⋅10⁻⁷ m² s, Δt=250s, the finest level consists of 33×33 mesh points. In all tested FMG and also in MG-V-cycles, full weighting restriction and bilinear interpolation have been applied, except for the diffusivity prolongation at the end of the coarse level iteration. SOR was the coarsest grid solver. The number of pre- and post-relaxation is equal to 3, the coarsest grid resolution has 5×5 points, for nonlinear iteration ε =1⋅10⁻⁴ and the initial norm reduction in MG-V is η=0.1.

Verifying the Coarse Level Iterating Algorithm

The finite difference methods (5) and (7) are compared without reduction of the diffusivity hysteresis, i.e. the diffusivity has been iterated in both cases at the finest resolution. Initial conditions are θ⁰ =1 with constant boundary conditions θ

(

x,y,0

)

=0,

∂

∈ y ,

x . Fig. 6a shows the norm of differences between the solutions obtained by two distinct difference schemes versus time. The highest norm value is less than 2.5⋅10⁻³.

2 4 6 8 10 12 14 16 18

0 0.5 1 1.5 2 2.5x 10^-3

Time step

norm

b c

Fig. 6. Norm tendencies a) differences between schemes (5) and (7), b) between CLI-MG by bilinear interpolation and fine level iteration, c) CLI-MG by prolongation and fine level iteration,

d) CLI-MG with Δt=500s and Δt=50s

Comparing the numbers of iteration cycles, it can be stated that the proposed scheme (7) generally needs two iterations, however formulation (5) at larger temperature differences needs considerably more iteration cycles [61], [65]. When approaching the steady-state, the cycle number could also decrease. In the proposed algorithm, after the coarse level iteration, the diffusivity remains unchanged. The crucial question is how the diffusivity should be prolonged to the finest grid? When dealing with two-dimensional cases the bilinear interpolation is the most popular technique. Another

possibility is the simple prolongation. According to the required rapid iteration, higher order interpolation is not suitable. Fig. 6b and Fig. 6c show the norms of differences between the fine grid iteration and coarse level iteration with two distinct diffusivity projections. With the deviation in the norm with the prolongation algorithm (curve c) all the time remains beneath the norm of bilinear interpolation (curve b). Comparing the results, it can be concluded, that the algorithm with prolongation approximates the results better using the fine grid iteration algorithm, therefore I propose that the temperature field can be interpolated but the diffusivity field should be prolonged to the finest grid.

The correctness of the time-centred diffusivity has been tested by employing different time steps using the same resolutions in space. The reference time step was a tenth of the original time step. Fig. 6d shows the norm of differences between two steps in time. Apart from the first few steps, there is no difference between the two calculations, so time averaging is suitable, until regularity and consistency are ensured.

Comparing the solution time, I measured twenty time steps. Assuming that the solution time for the coarse level iteration is one unit, the solution time with fine grid iteration is four units. Analysing the memory usage, the memory requirement of the multigrid part of the two methods are roughly the same. The discrete Preisach algorithm in the proposed form has one large distribution matrix. The memory character of the Preisach model requires to store the previous input values of all grid points in a matrix of size

(

n+1

) (

× n+1

)

and the return points of the characteristic of θ_ij

( )

kΔt at each grid point, which means that about 30 values have to be stored and treated per cycle and per grid point. If the Preisach algorithm works on the next coarser grid, the memory demand can be calculated in the same way, replacing n by N =n 2. The saved memory size if

>>1

n is about ≈23n² data. With large problems, significant differences in memory usage can be accompanied by reduced CPU-time as well [65].

Hysteresis Induced Heating-Cooling Asymmetry

To investigate the results of simulations using the proposed algorithm, the heating-cooling thermal cycles with constant and temperature dependent hysteretic thermal diffusivity are compared. Numerical simulations well represent the thermal cycling asymmetry induced by hysteresis of diffusivity. In the first test the thermal diffusivity changes on the main hysteresis curves. It can be reached by monotone heating or cooling of the boundaries. In the heating process the initial and boundary conditions are

( )

⁰

θ x,y,0 = at t=0, and ^θ

(

^x^,^y ^t,

)

⁼¹^, ^x^,^y^∈^∂^Ω ^att^>⁰. The initial and boundary conditions in the cooling process are θ

(

x,y,0

)

=1 at t=0; θ

(

x,y t,

)

=0, x,y∈∂Ω and

t . In the compared test cases, the constant diffusivity is considered to be the base diffusivity of the hysteresis diffusion function κ₀.

Temperature fields and diffusivity fields at four different time steps are shown in Fig. 7.

Diffusivity Temperature

with hysteretic diffusivity Temperature

with constant

diffusivity ^s

m²

10⁻7 Temperature Diffusivity

with hysteretic diffusivity Temperature

with constant

diffusivity ^s

m² 10⁻7

Fig. 7. Heating and cooling a rectangular 2D sample with constant and with hysteretic diffusivity

0 50 100 150 200

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Time step

heating

cooling κ=hyst.

κ^=hyst.

κ=const.

Fig. 8. Temperature transients in the middle of the domain with hysteretic and with constant diffusivity

Heating-cooling transients with constant thermal diffusivity are symmetrical. In hysteretic diffusion, the heating-cooling transients differ from each to other. The thermal process speeds up with increased diffusivity. The cooling process begins at higher diffusivity [69]. Differences between constant and hysteretic diffusion processes are expressed by the Euclidean norm. Fig. 9 shows the norm of differences for heating and cooling processes. The norm first increases in time and after reaching a maximum slowly decreases. In the transient time domain, the difference between cooling processes is higher than the difference between heating processes. The cooling process with hysteretic diffusivity approaches the steady-state more rapidly.

0 50 100 150 200 0

1 2 3 4 5 6 7 8 9x 10^-3

Time step

norm

Fig. 9. Norm of differences between hysteretic and normal heat diffusion a) heating and b) cooling along the main hysteresis curves

Periodically Changing Boundary Temperatures

The hysteresis minor loops have been avoided using the presented monotone heating or cooling models. However with periodically heating-cooling of the boundaries, it has to be assumed, that diffusivity could vary along minor loops. In this test case, a rectangular domain with sinusoidal varying temperatures at two neighbouring boundaries, zero temperature boundary conditions on the other two sides of the domain are considered.

The transient thermal characteristic of a material with diffusivity hysteresis is compared to the characteristic of a material with constant diffusivity. Hundred time steps have been evaluated and are now reported. The temperature fields and diffusivities at time steps 100k=10,20,50, can be seen in Fig. 10. In Fig. 10a-d the linear, in Fig. 10e-h the hysteresis heat diffusion time series and in Fig. 10i-l the appropriate diffusivity fields can be seen. Diffusivity minor loops are shown in Fig. 11. The periodical changing of boundary temperatures can be seen in Fig. 12. When determining the norm of differences of the nonlinear and the linear temperature fields in time, it can be concluded, that the norm is increasing rapidly when the boundary temperatures are increasing and is decreasing slowly when the boundary temperatures are decreasing.

The characteristic of the curve approaches a quasi-static state. The asymmetry in the norm cycles well represents the nonlinear effects of the hysteresis.

0.8

0.0 0.4

6.0E-7

2.5E-7 4.2E-7

θ, κ=constant κ(m²/s)

←k=100→

←k=50→

←k=20→

←k=10→

θ, κ=hysteresis

a) e) i)

j) b) f)

c) g) k)

h) l) d)

Fig. 10. Temperature fields with a)-d) linear, e)-h) nonlinear heat diffusion and i)-l) diffusivity at four different time steps

0 0.2 0.4 0.6 0.8 1

2.5 3 3.5 4 4.5 5 5.5 6 6.5

7x 10^-7

θ κ (m2/s)

main loop minor loops 1 minor loops 2

Fig. 11. Minor loops of diffusivity near the boundaries

0 20 40 60 80 100 0

0.5 1 1.5 2 2.5x 10^-3

norm

0 20 40 60 80 100

0 0.5 1

Time step

Fig. 12. Top: Norm of differences between linear and nonlinear heat diffusions. Bottom: Periodically changing boundary temperatures

Nonlinear Diffusion with Inhomogeneous Initial Conditions

Special problems arise from inhomogeneous initial conditions. In these cases the initial diffusivity field depends on how the initial temperature distribution has been reached. A hypothetical initial temperature distribution can be seen in Fig. 13a. To have a correct solution, the initial diffusivity field has to be determined as well. If the initial diffusivity field can be calculated on the main hysteresis curves, different initial diffusivity fields can be found on the increasing and decreasing main hysteresis curves, see Fig. 13b and Fig. 13c, respectively.

In document Nemlineáris hőterjedés és hőképanalízis (Pldal 13-0)