Avalanche dynamics in higher dimensional fiber bundle models
Zsuzsa Danku Geza Odor Ferenc Kun





 
 
  
 Corresponding author: ferenc.kun@science.unideb.hu
  ^1Department of Theoretical Physics, University of Debrecen,
 P.O. Box 5, H-4010 Debrecen, Hungary
 ^2MTA-MFA-EK,  Center  for  Energy  Research  of  the  Hungarian  Academy  of  Sciences,
H-1121  Budapest,  P.O. Box  49,  Hungary


We investigate how the dimensionality of the embedding space affects 
the microscopic crackling dynamics and the macroscopic response
of heterogeneous materials. Using a fiber bundle model with localized load sharing 
computer simulations
are performed from 1 to 8 dimensions slowly increasing the external load up to failure.
Analyzing the constitutive curve, fracture strength and avalanche statistics of bundles 
we demonstrate that a gradual crossover emerges from the universality class 
of localized behavior to the mean field class of fracture as the embedding dimension 
increases. The evolution between the two universality classes 
is described by an exponential functional form.
Simulations revealed that the average temporal profile of crackling avalanches 
evolves with the dimensionality of the system from a strongly asymmetric shape to 
a symmetric parabola characteristic for localized stresses and homogeneous stress 
fields, respectively.
 
Avalanche dynamics in higher dimensional fiber bundle models
    Ferenc Kun^1
    
============================================================





§ INTRODUCTION

A large variety of heterogeneous materials respond to a slow external driving 
in a jerky way where sudden outbreaks of activity are separated by silent periods 
<cit.>.
From the propagation of imbibition fronts in heterogeneous materials <cit.>, 
through dislocation bursts of plastically deforming crystals <cit.>,
and Barkhausen noise in ferromagnets <cit.>, to fracture phenomena 
<cit.>
and earthquakes <cit.>, crackling noise has been observed over 
a broad range of length scales.
It was found that crackling noise
is characterized by scaling laws, i.e. the statistics of the quantities of single 
bursts is described by power law distributions which may be the fingerprint of an underlying phase 
transition <cit.>. 

Recently, it has been pointed out that the average temporal profile is a fundamental 
feature of crackling avalanches. Experimental and theoretical studies have revealed that 
the precise shape of the average profile 
of bursts is sensitive to the details of the physics of the system and it encodes 
valuable information about the underlying intermittent dynamics of pulse generation 
<cit.>.
For the fracture of heterogeneous materials careful experiments have been performed
where the temporal evolution of individual bursts formed at a propagating 
crack front was determined by direct optical
observation using high speed imaging <cit.>. 
These investigations provided symmetric 
parabolic profiles mainly attributed to long range elastic forces acting along 
the crack front <cit.>. 
Measurements of magnetic noise induced by the dynamic fracture of steal revealed similar
pulse profiles, however, with a right handed asymmetry <cit.>.
The front propagation was modeled 
as the driven motion of an elastic line in a disordered environment of pinning centers.
Varying the range of interaction it was found that the degree of asymmetry depends on the 
range of stress redistribution, i.e. profiles evolve with the universality class 
of fracture from a strongly asymmetric shape (localized interaction) to a symmetric parabola 
(long range interaction) <cit.>. Simulation studies of the dynamics 
of breaking bursts in the 
fiber bundle model underlined the general validity of this behavior 
<cit.>.

In the present paper we take the opposite strategy and address the question 
how the dimensionality of the sample affects the fracture process when the range of 
interaction is kept constant. 
We performed computer simulations in the framework of a fiber bundle model 
with nearest neighbor load sharing after local failure events varying the 
dimensionality of the system from 1 to 8. Both on the macro- and micro-scales
the system exhibits a crossover between the universality classes of localized 
behavior and the mean field class of fracture phenomena. We show that 
this evolution is described by a genuine exponential form.
The temporal profile of breaking avalanches can be well described
by the scaling form suggested in Ref. <cit.> where 
the parameters clearly confirm the crossover between the two universality classes.
Our study shows that the upper critical dimension of the fracture of heterogeneous materials 
is infinite in agreement with a recent theoretical prediction  <cit.>.
We give numerical evidence that the critical exponents change as an exponential function of the 
dimension.



§ LOCAL LOAD SHARING FIBER BUNDLE MODEL IN 1 TO 8 DIMENSIONS

The fiber bundle model provides an efficient modeling framework
for the fracture of heterogeneous materials <cit.>. 
In spite of its simplicity it captures 
the essential ingredients of fracture phenomena
allowing also for analytical solutions for the most important quantities 
<cit.>.
The classical fiber bundle model consists 
of N parallel fibers which are organized on a regular lattice. 
In D=1 fibers are placed
equidistantly next to each other along a line, while in D=2 the fibers are 
assigned to the sites of a square lattice of side length L. 
The fibers are assumed to have a perfectly brittle behavior, i.e. 
they exhibit a linearly elastic response with a Young modulus E up to breaking 
at a threshold load σ_th. The Young modulus is assumed to be 
constant E=1 such that the disorder of the material is solely represented 
by the randomness of the breaking threshold σ_th:
to each fiber a threshold value is assigned σ_th^i, i=1,… , N 
sampled from the probability density p(σ_th). In the present calculations 
we used exponentially distributed breaking thresholds 

p(σ_th) = 1/λ e^-σ_th/λ


over the range 0≤σ_th<+∞. 
The scale parameter was fixed to λ = 2 in all the calculations.

In D≤ 2 the bundle is loaded in the direction parallel to the fibers,
which represents the uniaxial loading of a bar shaped specimen. 
Under stress controlled 
loading, when the local load on a fiber reaches its failure strength the fiber breaks
and its load has to be redistributed over the remaining intact fibers. We assume 
localized load sharing (LLS) so that the load of a broken fiber is redistributed 
equally over its intact nearest neighbors in the square lattice 
<cit.>. When the breaking 
fiber is entirely surrounded by intact ones in the square lattice, four fibers share the load, 
however, when the breaking fiber is at the perimeter of a growing broken cluster 
(crack) typically 2-3 fibers receive the excess the load. 
As a consequence, high stress concentration builds up along 
the perimeter of cracks and local stress fluctuations develop.

In higher dimensions D>2 the generalization of the model is straightforward 
although it does not have a direct relevance for practical applications: the 
fibers are assigned to sites of cubic lattices and the load on them is represented by 
a scalar variable. 
After failure events the load is redistributed
over the intact nearest neighboring sites along the edges of the lattice.
The emerging stress concentration is controlled by the coordination number z 
of the underlying lattice which depends on the embedding dimension as z=2D in our setup.
In all dimensions periodic boundary conditions are implemented in all lattice directions.

Computer simulations were performed by quasi-statically increasing the external 
load σ which is realized by increasing σ to provoke the breaking 
of a single fiber. After the
fiber has been removed its load gets redistributed according to the rules discussed above.
The enhanced load on the neighboring fibers may induce further breaking which is followed 
again by a local stress redistribution. As the result of this repeated breaking and stress
redistribution the breaking of a single fiber can give rise to an avalanche of breakings.
The loading process stops when a catastrophic avalanche is triggered in which all remaining
intact fibers break.

In order to keep the problem numerically tractable for the dimensions 
D=1, 2, 3, 4, 5, 6, 7, 8 the lattice size was set to L=4084101, 2021, 159, 45, 21, 13, 9, 7,
which ensures nearly the same number of fibers in all dimensions. 
To obtain reliable results statistical averaging was done over K=5000 simulations.




§ MACROSCOPIC RESPONSE

The mean field limit of FBMs is realized by the equal load sharing (ELS) 
of the load of broken fibers over the intact ones. Under ELS conditions all intact 
fibers receive the same amount of load irrespective of their distance from the 
broken one. It follows that no stress fluctuations can emerge, all fibers keep
the same load during the entire loading process.

Hence, in the mean field limit the random strength of fibers is the only source 
of disorder in the system.
For ELS the macroscopic stress-strain relation σ(ε) 
of the bundle can simply be obtained from 
the general expression σ=Eε[1-P(Eε)]. Here P(x) denotes 
the cumulative distribution of the failure thresholds so that the term 1-P(Eε)
yields the fraction of intact fibers which all keep the same load Eε. 
Substituting the exponential distribution Eq. (<ref>) we obtain 

σ= Eεe^-Eε/λ.


In Fig. <ref> the curve of Eq. (<ref>) is presented 
up to the maximum where catastrophic failure occurs under stress controlled loading.
The constitutive curve has a quadratic maximum the value σ_c and position
ε_c of which determine the fracture stress and strain of the bundle, respectively

σ_c^MF       =    λ/e, 

ε_c^MF    =    λ/E.

Note that the fracture strength σ_c and ε_c of FBMs depend on the 
system size N even in the mean field limit 
<cit.>. 
Although the convergence is rapid, strictly
speaking the above expressions give the bundle strength in the limit of infinite system size.

For finite dimensions D the σ(ε) curves were determined by computer 
simulations averaging 
over a large number of loading processes with different realizations of the threshold 
disorder. It can be seen in Fig. <ref> that for all dimensions D the mechanical 
response σ(ε) of LLS bundles follows the mean field solution Eq. (<ref>). 
However, for low dimensional bundles the curves stop significantly earlier implying a 
lower fracture strength and a higher degree of brittleness. As D increases 
the LLS constitutive response completely recovers the mean field behavior.
This tendency becomes more transparent when analyzing the average fracture stress <σ_c>
and strain <ε_c> as a function of D. In
Figure <ref> both quantities gradually converge to their mean field 
counterpart, however, the convergence is somewhat faster for <σ_c>.
For the highest dimension D=8 considered only a few percent difference is observed
from the value of σ_c^MF of Eq. (<ref>).

The most remarkable result is that the convergence to the mean field limit 
can be described by an exponential form. The inset of Fig. <ref> 
demonstrates that subtracting a constant σ_c^* from <σ_c> a straight 
line is obtained on a semi-logarithmic representation which implies the functional form 

<σ_c>(D) = σ_c^* +Aexp(-D/D^*).


Here D^* denotes a characteristic value of the dimension.

Formally, σ_c^* is a free parameter in Eq. (<ref>) which 
was tuned to σ_c^*=0.725± 0.015 to obtain the best straight line 
in Fig. <ref>. Note that this value falls very close to σ_c^MF=0.735.
The characteristic dimension was obtained by fitting D^*=1.65± 0.06.

The results show that in higher dimensions the role of stress fluctuations is
diminishing in the fracture process and the behavior of the system gradually 
approaches the one of the completely homogeneous stress field 
of the infinite dimensional ELS solution.
It follows that, in spite of the highly localized stress redistribution,
in high D the stochastic breaking process is completely controlled by the 
quenched disorder of the failure strength of fibers.



§ AVALANCHE DYNAMICS

The bundle is loaded in a quasi-static way such that the external load is increased to provoke
solely the breaking of a single fiber. In the simulations the cascade of breaking fibers, emerging 
due to the repeated steps of load redistribution and breaking, is followed 
until it stops while keeping the external load fixed. The size Δ of
avalanches is characterized by the number of fibers breaking in the avalanche.
These breaking avalanches are analogous to crackling bursts measured in experiments 
with acoustic <cit.> or electromagnetic 
<cit.> 
techniques. Under simple 
geometrical conditions such as during the slow propagation of a planar crack direct 
optical observation has also proven very successful 
<cit.>.

The probability distribution p(Δ) of avalanche sizes is presented in Fig. 
<ref>.

In all dimensions D the probability density p(Δ) is described by the same 
functional form, i.e. power law distributions are obtained followed by an exponential 
cutoff

p(Δ) ∼Δ^-τexp(-Δ/Δ^*),


where both the exponent τ and the characteristic burst size Δ^*
depend on the dimensionality D of the bundle. In Fig. 
<ref> Eq. (<ref>) provides excellent 
fits of the numerical data where the exponent τ decreases while 
Δ^* increases with D. The result implies that in higher dimensions 
the system can tolerate larger avalanches without catastrophic collapse.
It has been shown by analytical calculations that in mean field FBMs the 
size distribution exponent takes the value τ_MF=5/2 which has proven to be 
universal for a broad class of threshold distributions 
<cit.>. 

For our LLS system it can be observed in Fig. <ref> that τ has high values 
in low dimensions but with increasing D it approaches the mean field exponent.
The inset of Fig. <ref> demonstrates that the convergence 
is again described by an exponential form similar to Eq. (<ref>) 

τ(D) = τ^*+Bexp(-D/D^*),

where best fit was obtained with the same value of the characteristic dimension 
D^*=1.65± 0.06 as for the fracture strength Eq. (<ref>). 
The value of τ^* providing the best straight line is τ^*=2.52± 0.04 
which falls very close to the mean field exponent τ_MF.



§ TEMPORAL PROFILE OF AVALANCHES

Recently, we have shown that breaking avalanches in fiber bundles have a complex 
time evolution <cit.>: an avalanche typically starts
with the breaking of a single fiber which in turn triggers the breaking of 1-2 additional
fibers after load redistribution. The subsequent load redistribution steps involve a 
larger and larger number of fibers giving rise to a spatial spreading of the avalanche.
The avalanche stops when all the fibers involved in the last redistribution step 
are able to sustain the elevated load. This dynamics implies that avalanches 
are composed of discrete growth steps of size Δ_s, which is the number of fibers 
breaking in a single load redistribution step. The total number of subsequent redistribution -
breaking steps defines the duration W of the avalanche.
The time evolution of a single burst of size Δ = 8705 and duration W=264 
is illustrated in Fig. <ref> for a three-dimensional bundle, where 
cubes represent fibers. The color code corresponding to the growth steps of the avalanche 
facilitates to follow the breaking sequence.

The temporal evolution of an avalanche is characterized by the Δ_s(u) function, 
where u is a time variable taking integer values in the interval u=1,… , W.
Similarly to the size of bursts Δ, their duration W is also a stochastic 
quantity which varies over a broad range. 

It can be observed in Fig. <ref> 
that the probability distribution p(W) of the burst duration W has the same 
functional form as p(Δ), i.e. power law behaviour followed by an exponential 
cutoff is evidenced. In higher dimensions
the system can tolerate larger bursts of longer duration, hence, the power 
law exponent τ_W of p(W) decreases from τ_W=5.7 of D=1 to 
the vicinity of τ_W≈ 4 for D=8, while the cutoff duration gradually increases
with the dimensionality.

For single bursts Δ_s(u) is a stochastic curve, hence,
quantitative characterization of the internal dynamics of avalanches is provided by the
average temporal profile <Δ_s(u,W)>, where the step size
Δ_s is averaged at fixed values of u for avalanches of the same duration W 
<cit.>.
Average profiles are illustrated in Fig. <ref> 
for different durations W for all dimensions considered. It can be observed
that, except for D=1, the shape of avalanches obtained is similar to the experimental findings
<cit.>. 
For D=1 the stress concentration is so high at the tip of growing broken clusters that 
all steps of the breaking sequence have a size Δ_s=1, since a larger number of breaking fibers 
would trigger 
a catastrophic avalanche. Consequently, the emerging pulse profile is completely flat.
At low dimensions D=2,3 the profiles have a strong right handed asymmetry at all durations W, 
which means that these bursts start slowly, they gradually accelerate while
their stopping is more sudden. 
As the dimensionality of the system increases the degree of asymmetry decreases
and eventually a symmetric parabolic shape is obtained which is characteristic for 
mean field avalanches 
<cit.>. 
Recently, similar asymmetric avalanche shapes have been obtained in a fiber bundle model 
of creep rupture with localized load sharing <cit.>. Since in those 
calculations the external load was fixed, bursts were triggered by ageing induced 
slow breaking of fibers. However, the cascading breaking sequence of avalanches had essentially
the same dynamics as in the present study.

It has been shown in Ref. <cit.> for D=2 that the 
short range load sharing and the heterogeneous stress field built up along the
perimeter of cracks are responsible for the right handed asymmetry of pulse profiles.

Comparing the curves in Figs. <ref>(a,b,c,d) for a fixed dimension 
it can be inferred that bursts of a longer duration W have a larger average height <Δ_s^max> 
and average size <Δ>. Figure <ref> demonstrates that 
rescaling the pulse profiles with an appropriate power of W pulses of different duration 
can be collapsed on the top of each other. The good quality data collapse implies
the scaling form

<Δ_s(u,W)>=W^αf(u/W),


where the value of the exponent α and the scaling function f(x) both
depend on the dimensionality of the system D. In Figure <ref>
the value of α was tuned to achieve the best collapse.
Deviations from the scaling function f(x) occur for the shortest durations 
which confirms that Eq. (<ref>) is asymptotically valid.

It follows from the scaling structure Eq. (<ref>) 
that the average avalanche size has a power
law dependence on the duration <cit.>

<Δ> ∼W^1+α, 


which provides an alternative way to determine the exponent α, as well.
Figure <ref> shows that the asymptotic behaviour of the average 
size of bursts <Δ> of the same duration W can be well described 
by a power law. 
The value of the exponent 1+α obtained by fitting 
the <Δ>(W) curves with Eq. (<ref>)
is presented in Fig. <ref> as a function of the embedding dimension D. 
Since in D=1 the step size Δ_s 
does not exceed 1, the total size Δ of an avalanche is proportional 
to its duration so that 1+α=1 and α=0 follows. As D increases, 1+α 
gradually approaches 2, and hence α tends to its mean field value α=1.
Based on Eq. (<ref>) a scaling relation 
can be established between α and the exponents τ and τ_W of the probability 
distribution of the size and duration of bursts

1+α=(τ_W-1)/(τ-1) .

In Figure <ref> 
we compare  the prediction of the scaling law by substituting
the numerical values of τ and τ_W on the right hand side of Eq. (<ref>), 
to the value of 1+α obtained by direct fitting of the average burst size 
in Fig. <ref>.
A reasonable agreement can be observed between the two curves which confirms the consistency
of the results.

Avalanche profiles <Δ_s(u,W)> have an asymmetric functional 
form in all dimensions, however, the degree of asymmetry depends 
on the embedding dimension D.
In Ref. <cit.> the following expression has recently been 
suggested to quantify avalanche shapes

f(x) ∼[x(1-x)]^α[1-a(x-1/2)].


Note that the scaling laws Eqs. (<ref>,<ref>)
are consistent with the generic form of Eq. (<ref>) with the same
value of the exponent α.
The pulse asymmetry is represented by the parameter a
such that zero value of a
implies symmetry, while negative and positive values of a capture right and left 
handed asymmetry, respectively. It can be observed in Fig. <ref>
that Eq. (<ref>) provides excellent fits of pulse profiles for all dimensions.
Figures <ref> and <ref> present the value of the 
parameters α and a 
obtained by fitting for several avalanche durations W. The careful analysis revealed that
the parameter values practically do not depend on the avalanche duration W, except for some
statistical fluctuations the numerical values of α and a agree with each other

at different Ws. In agreement with the results presented in Fig. <ref>, 
for low dimensions the exponent α starts from the vicinity of 0.6-0.7
and it increases to 0.95-1 at high dimensions (see Fig. <ref>). 
At the same time the observed right handed asymmetry of profiles
in Fig. <ref> gives rise to negative values of a in
Fig. <ref>. As the event duration W increases from 10 to 80, 
in Figs. <ref> and <ref> the numerical error of the determination 
of the parameters α and a increases 
from 0.05 to 0.15, and from 0.07 to 0.16, respectively. The reason is that 
due to the rapidly decreasing duration distribution p(W) of bursts, the statistics of
events of higher duration decreases.
In higher dimensions α approaches one
and a increases to the vicinity of zero which imply that the profile shape evolves
towards a simple symmetric parabola of mean field crackling systems 
<cit.>.
The inset of Fig. <ref> demonstrates that subtracting α from 
a proper limit value α^* straight lines are obtained on a semi-logarithmic plot.
For each W the value of α^* was varied in the range 0.96-1.01 
independently until best straight lines were achieved. Similarly, the inset of Fig. <ref>
presents the absolute value of a where again an exponential dependence is evidenced.
The results imply that the behaviour of the pulse parameters α and a is consistent with 
that of the macroscopic strength and avalanche exponent, i.e. they approach their mean field
values with an exponential dependence on the dimensionality. 




§ DISCUSSION

The fracture of heterogeneous materials is strongly affected by the degree of disorder,
which controls the spatial variation of microscopic materials' strength, 
and the stress fluctuations emerging due to localized stress redistribution after 
micro-fracturing events. The competition of these two sources of disorders, 
i.e. strength and stress disorders,
gives rise to a highly complex fracture process which manifests
itself in the variation of macroscopic strength, in the
statistics of crackling bursts, and in the temporal evolution of single crackling events.
At a fixed amount of quenched strength disorder stress fluctuations are mediated
by the range of load redistribution. Former studies have revealed two universality classes
of fracture, i.e. the localized load sharing (LLS) class and the mean field class (ELS),
characterized
by a high degree of brittleness and a quasi-brittle response with a large amount of 
avalanche precursors of failure, respectively.

In order to understand the competing role of different disorder sources between the 
limiting cases of the LLS and ELS classes, in our study we considered a fiber bundle 
model with a fixed amount of strength disorder and varied the dimensionality of the system
from 1 to 8 at a fixed range of load sharing. The strength disorder is represented by an exponential
distribution of the failure threshold of fibers, which falls in the well understood universality
classes both in the ELS and LLS limits of load sharing. In all dimensions nearest neighbor
load redistribution was implemented on cubic lattices with periodic boundary conditions 
in all lattice directions. 


Our study revealed a very interesting dimensional crossover between the two universality 
classes of fracture: Both on the macro and micro-levels fracture characteristics 
evolve with the dimensionality from the highly brittle response of low dimensional systems
controlled by stress fluctuations, to the quasi-brittle behavior in high dimensions where the
strength disorder dominates. For the macroscopic strength of the bundle and for the 
power law exponent of the size distribution of crackling bursts the convergence to the mean 
field limit is described by an exponential functional form.


Avalanche profiles have been found before to be very sensitive to models' details, especially
to the degree of correlations of microscopic events leading to collective avalanches 
of local failures. We showed that the symmetry of avalanche profiles depends on the dimensionality
of the system, gradually shifting from a strongly asymmetric shape at low dimensions to a symmetric 
parabolic form in the mean field limit. The parameters of pulse profiles were found to evolve towards
their mean field values with an exponential dependence on the dimensionality similarly to
the macroscopic strength and avalanche exponent of the system. 

Increasing the dimensionality of the bundle implies a decreasing stress concentration
along broken clusters by increasing the connectivity of the system. 
As a consequence, stress fluctuations
have a diminishing role with increasing dimension giving more room 
for the quenched disorder of fibers' strength. 
The evolution of avalanche temporal profiles 
we observe with the embedding dimension is similar to what has been obtained recently for 
avalanches of a propagating crack front when increasing the range of interaction
in a fracture model of fixed D=2 dimensions <cit.>. 


Recently, the effect of the long range connection 
of fibers on the fracture process has been studied. Instead of a regular lattice,
fibers were assigned to the nodes of a complex network with small world 
properties <cit.>. 
Redistributing the load along the edges of the network simulations revealed 
that a small amount of long range connection is sufficient to converge the system 
from the LLS to the mean field universality class. Eventually, the dominance of 
quenched structural disorder is responsible  for the ELS behavior similarly to our case.


Fiber bundles in higher dimensions have recently been investigated in Ref. 
<cit.>. Focusing on the amount of damage accumulated 
up to failure and on the distribution of avalanches a crossover from the LLS 
to the mean field universality class was pointed out. However, the crossover 
is described by a power law functional form contrary to our exponential 
behavior. This difference may arise due to the different protocols of 
load redistribution used in the simulations. 
In agreement with Ref. <cit.>, our   
results imply that the upper critical dimension of the fracture of heterogeneous materials 
is infinite. Further support of this remarkable feature of fracture phenomena 
could be obtained by a finite size scaling analysis in higher dimensions, however,
it turned to be infeasible due to the overwhelming numerical costs.

In the broader context of critical phenomena, the absence of a finite upper critical
dimension has also been found in Kardar-Parisi-Zhang type surface growth models 
<cit.>,
where the critical exponents present an approximate
exponential dependence on the dimensionality of the system <cit.>.
We conjecture that the quenched disorder of the system and the locally conserving nature 
of the dynamics (i.e. the entire load dropped by broken fibers is redistributed over the intact ones 
without loss) are responsible for the absence of a finite upper critical dimension beyond 
which mean field behavior is attained.
The exponential crossover from the local to the mean field universality class of fracture involves 
a characteristic dimension D^* which falls between 1 and 2. Since avalanche shapes of the one 
dimensional system do not conform with the higher dimensional ones, we propose the interpretation
of D^* as the lower critical dimension of fracture phenomena. Our simulations confirm
that for D>D^* all characteristic quantities of the system evolve through gradual quantitative
changes but the qualitative behaviour remains robust as the embedding dimension increases.



The work is supported by the EFOP-3.6.1-16-2016-00022 project. 
The project is co-financed by the European Union and the European Social Fund.
This research was supported by the National Research, Development and
Innovation Fund of Hungary, financed under the K-16 funding scheme Project no. K 119967
and K 109577. This work was supported through the New National
Excellence Program of the Ministry of Human Capacities
of Hungary.
The research was financed by the Higher Education Institutional
Excellence Programme of the Ministry of Human Capacities in Hungary, 
within the framework of the Energetics thematic
programme of the University of Debrecen.




