Mathematical Modeling of the Contact Interaction of the Fuel

Element Section, Including up to 100 Pellets,

Taking into Account Creep

PAVEL ARONOV1,2, MIKHAIL GALANIN1,2, ALEXANDR RODIN1,2

1Keldysh Institute of Applied Mathematics (Russian Academy of Sciences)

Miusskaya sq., 4 Moscow, 125047

RUSSIA

2Bauman Moscow State Technical University

Baumanskaya 2-ya st., 5/1, Moscow, 105005

RUSSIA

Abstract: An algorithm for solving axisymmetric contact problems of thermoelasticity taking into account creep

processes is considered. To account for the contact interaction of bodies, the mortar method was used, and the

modiﬁed Jacobi method was used to solve the resulting system of linear equations. For the numerical solution

of the problem simulating the creep process, time discretization based on the implicit Euler method was applied,

and the Newton method was used to linearize the resulting system of equations. An algorithm with automatic step

selection based on obtaining an estimate of the local error of the method is proposed. The results of applying

the proposed algorithm to a demonstration problem simulating thermomechanical processes in a section of a fuel

element comprising from 1 to 100 fuel pellets are presented.

Key–Words: contact problem of the thermoelasticity theory, ﬁnite element method, mortar method, creep, fuel

element

Received: March 29, 2022. Revised: April 15, 2022. Accepted: May 14, 2022. Published: July 4, 2022.

1 Introduction

Taking into account the contact interaction of various

functional components of the equipment is an impor-

tant component in assessing the stress-strain state of

most structures. It is impossible to obtain an analytical

solution for practically important contact problems,

and numerical methods are used to determine the dis-

placement and stress ﬁelds. Therefore, the problem

of creating new effective algorithms and modern ap-

plication software based on them for solving contact

problems of computational thermomechanics is ur-

gent. The numerical solution of such problems can be

carried out using the domain decomposition method

[1], the penalty method [2], various variants of the La-

grange multiplier method [3], in particular, the mortar

method [4].

The paper presents the formulation of a quasi-

stationary multicontact problem and presents an algo-

rithm for the numerical solution of such problems tak-

ing into account changes in time of deformations and

stresses due to the creep effect. The application of the

proposed algorithm to a problem simulating processes

in a section of a fuel element operating in the mode of

constant heat dissipation power is considered. To ac-

count for the creep effect, an algorithm based on the

implicit Euler method in combination with the New-

ton method for linearization is considered. Calcula-

tions with constant and variable time steps were car-

ried out and the results obtained were compared.

2 Mathematical formulation of the

problem

Let us limit ourselves to solving the following axisym-

metric problem simulating thermomechanical pro-

cesses occurring in a fuel element: inside the cylin-

drical cladding GNthere is a column of N−1identi-

cal cylindrical pellets G1, . . . , GN−1stacked on top of

each other, having an inner hole and chamfers at both

ends. Fig. 1 shows the modeling area corresponding

to half of the longitudinal section of the fuel element

section. We introduce the following notation for the

surfaces of bodies: S1is the surface area on which the

1st kind condition for the normal displacement com-

ponent is set, S2is the inner surface of the pellets, S3

is the upper end of the upper pellet, S4is the outer sur-

face of pellets, S5is the inner surface of the cladding,

S6is the outer surface of the cladding.

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

Figure 1: Scheme of the area model

Suppose that the coupling effect can be neglected,

so we will solve the problem of thermal conductivity

separately, and use the resulting temperature ﬁeld to

solve the contact problem of thermomechanics. When

solving the thermal conductivity equation, the inﬂu-

ence of the resulting mechanical deformations is not

taken into account.

Consider the following initial boundary value

problem for a nonlinear heat equation in the domain

α=1

Gα,t > 0[5]:

c(T)ρ∂T

∂t = (kij (T)T,j ),i +q(x, t),x∈G;(1)

T(x,0) = T0(x),x∈G;(2)

T(x, t) = T1(x),x∈S6;(3)

−nikij (T)T,j (x, t)=0,x∈∂Gp\S4;(4)

−nikij (T)T,j (x, t) = α[T(x, t)−

−Tf(¯

x, t)],x∈S4,¯

x∈S5,(5)

where c(T)is the speciﬁc mass heat capacity of the

medium, ρis density of the medium, xis a vector of

spatial coordinates, tis time, kij are the components

of the thermal conductivity tensor, T,j =∂T

∂xj

,q(x, t)

is the power of internal heat sources (different from

zero in pellets), T0(x)is initial temperature, T1(x, t)

is surface temperature S6,T(x, t)is temperature at

time t,niare the components of the unit vector of

the external normal to the boundary ∂G,αis the heat

transfer coefﬁcient on the surfaces S4and S5,Tf(x)

is temperature at a similar point lying on the inner sur-

face of the cladding.

The mathematical formulation of the quasi-

stationary problem of deformable solid mechanics for

the case when there are no volumetric forces includes

[5] the following relations for each body with the

number α(i, j = 1,3), t > 0:

• equilibrium equations

σji,j (u, t) = 0,x∈Gα;(6)

• kinematic boundary conditions

u(x, t) = u0(x, t),x∈SD;(7)

• force boundary conditions

σji(u, t)nj=pi(x, t), x ∈SN;(8)

• Cauchy relations for a linear full strain tensor

εij (x, t) = 1

2(ui,j (x, t) + uj,i(x), t),x∈Gα;

(9)

• deﬁning equations (Hooke’s law)

σij (x, t) = Cijkl(εkl(x, t)−ε0

kl(x, t)),x∈Gα,

(10)

where xiare the coordinates of the vector x∈Gα;

σij are the components of the stress tensor; εkl are

the components of the full strain tensor; ε0

kl are the

components of the inelastic strain tensor; uiare the

components of the displacement vector; Cijkl are the

components of the tensor of elastic constants; piare

the components of the vector of surface forces; njare

the components of the vector of the external normal to

the corresponding surface Sj;SDand SNis the union

of surfaces on which kinematic and force boundary

conditions are set, respectively.

In the calculations carried out, it was assumed that

the lower ends of the cladding and the lower pellet are

ﬁxed vertically, pi(x)differs from zero only on the

upper end of the upper pellet SN−1

3(constant pres-

sure 50 MPa) and on the outer surface of the cladding

4(constant pressure 10 MPa). The model took into

account that each pellet (except G1and GN−1) comes

into contact with two adjacent (top and bottom) pellets

and the cladding (it is assumed that there is no initial

gap between them). Thus, there are N−2pellet/pellet

contact pairs (Sα

1, Sα+1

3)and N−1pellet/cladding

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

contact pairs (Sα

4, SN

2). Sliding conditions without

friction are set on the contact surfaces.

Consider a pair of potentially contact surfaces re-

lated to bodies with numbers α1and α2. To sim-

plify the record, we will use the index ”1” instead

of ”α1” and ”2” instead of ”α2”. For pellet/pellet

contact pairs, the index ”1” denotes the upper end of

the lower pellet, and the index ”2” denotes the lower

end of the upper pellet, and for pellet/cladding con-

tact pairs, the index ”1” denotes the outer surface of

the pellet, and the index ”2” denotes the inner surface

of the cladding. Then the additional conditions on the

surface S1

kfor the case of frictionless contact are as

follows (for the surface S2

kconditions are written in

the same way):

σ1

τ(x1) = 0; (11)

σ1

n(x1) = σ2

n(¯x2)60; (12)

n(x1) + u2

n(¯x2)6δ0n(x1); (13)

σ1

n(x1)u1

n(x1)+u2

n(¯x2)−δ0n(x1)= 0.(14)

Here x1is a some point lying on the surface of S1

k, and

¯x2is a similar point, i.e. located opposite (normal to

k) the point x1on the surface of S2

k,δ0n(x1)>0

is the function that sets the initial gap (the surface

areas could not touch each other at the initial mo-

ment), uα

n=uα·nα,σα

τ= (σ(uα)·nα)·τα,

σα

n= (σ(uα)·nα)·nα.

The conditions (11)–(14) guarantee that if S12

surfaces S1

kand S2

kcoincide (in advance, the conﬁg-

uration and position of this section unknown), then

compressive contact forces will act on the contacting

areas. The conditions (11)–(12) are forceful, the con-

dition (13) is kinematic.

These conditions describe two possible situations:

1) in the ﬁnal conﬁguration, sections of potentially

contacting surfaces are in contact, kinematic condi-

tions (the corresponding coordinates coincide) and

force conditions (contact pressure modules coincide)

are met or 2) areas of potentially contacting surfaces

are not in contact: there is a gap between them, con-

tact pressures are equal to 0.

A characteristic feature of the considered problem

is that it considers many contact pairs, while there is a

strong change in the conﬁguration of the contact sur-

faces during the operation of the fuel element. For

each pellet/pellet contact pair, as a result of heating, a

signiﬁcant number of grid nodes located closer to the

chamfer come out of contact. Due to the large length

of the structure along the zaxis, the fuel pellets, espe-

cially in the upper part of the column, are displaced by

a considerable distance relative to their initial position

(and relative to the cladding).

To simulate the creep process, we will use the

ﬂow theory, the main provisions of which include [6]:

• additive decomposition of the derivative of the

full strain tensor in time

˙εij = ˙εe

ij + ˙εT

ij + ˙εc

ij ;(15)

• incompressibility of creep deformation

˙εc

ii = 0; (16)

• the relation for the derivative of the creep strain

tensor in time

˙εc

ij =µσ0

ij ,(17)

where µ=1

˙εc

σi,σi=q3

2σ0

klσ0

kl,˙εc

i=f(T, σi),

σ0

kl =σkl −1

3σjj δkl.

Solving the problem (6) – (10) at time tmis equiv-

alent to [7] minimizing the functional

Π = 1

Gα

σT(ε−ε0)dG −Z

uTpdS+

λn(u2n(x) + u1n(x)−δ0n)dS (18)

when the kinematic boundary conditions (7) are met,

where λn—Lagrange multipliers that are projections

of stress vectors on the directions of external normals,

un=urnr+uznz.

For spatial discretization of the functional (18),

the ﬁnite element method was used, second-order el-

ements on a quadrangular grid were used in calcula-

tions. The discretization of the integral over the con-

tact surfaces is performed using the mortar method

[8], which is the variant of the Lagrange multiplier

method for the case of inconsistent meshes. The al-

gorithm and some features of the application of the

mortar method are described in [9].

The choice of the most effective method for solv-

ing systems of linear equations arising during the dis-

cretization of the problem was carried out. For this

purpose, the contact problem was solved in a ther-

moelastic formulation for a large number of pellets

(up to 100). The use of modiﬁed iterative methods

(in particular, the modiﬁed Jacobi method) turned out

to be more effective than the use of standard meth-

ods for solving systems of linear equations (for ex-

ample, the Gauss method for sparse matrices or GM-

RES). The selected preconditioner [9] enables us to

automatically take into account the exit from the con-

tact without additional intervention in the algorithm.

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

Since stresses depend on creep deformations, the

equation (17) can be represented as the following sys-

tem of ordinary differential equations in vector form:

˙εc=f(εc, t).(19)

If the implicit Euler method is used to solve (19),

then at time tm+1 the formulas for the local vector are

valid (for the element with the number (e)) increments

of creep strain and local stress vector

∆εc(e)=εc(e)(tm+1)−εc(e)(tm) =

=τmµ(tm+1)σ0(e)(tm+1); (20)

σ(e)(tm+1) = D(e)(ε(e)(tm+1)−εT(e)(tm+1)−

−εc(e)(tm)−∆εc(e)),(21)

where D(e)is the elasticity matrix. The equation (18)

is nonlinear, the Newton method is used for lineariza-

tion [10].

The numerical algorithm can be described as fol-

lows:

1. For the moment of time tm+1, at the iteration

with the number i+ 1, after summing local vectors

and matrices into global ones, the following system of

linear algebraic equations is solved using a modiﬁed

Jacobi method [9]

Ki(tm+1)Ml

MlT0∆ui+1

∆λi+1=R(tm+1)

0

(22)

where ∆ui+1 =ui+1(tm+1)−ui(tm+1)is the vector

of increment of displacements per iteration, ∆λi+1 =

λi+1(tm+1)−λi(tm+1)is the vector of increments of

Lagrange multipliers per iteration, the matrix Mlre-

ﬂects the contribution of the contact surface integral to

(18), the matrix Kitakes into account the contribution

of temperature deformation and creep deformation, it

is calculated anew at each iteration and has the follow-

ing structure [10]:

Ki=





c1c2c20

c2c1c20

c2c2c10

000c3







,(23)

where c1=1

3(2C0+Cm),c2=1

3(Cm−C0),c3=

2C0,C0=1

aE+τkεc,Cm=E

1−2ν,aE=1 + ν

τm=tm+1 −tm.

2. For each ﬁnite element with index (e)by dis-

placement values in nodes u(i+1)(e)(tm+1)the defor-

mation vector ε(i+1)(e)(tm+1)and the stress vector are

calculated as follows:

σ(i+1)(e)(tm+1) = D(e)(ε(i+1)(e)(tm+1)−

−εT(e)(tm+1)−εc(i)(e)(tm+1)).(24)

3. The formula (20) calculates the values of

∆εc(i+1)(e)and new creep strain values are found:

εc(i+1)(e)(tm+1) = εc(e)(tm)+∆εc(i+1)(e).(25)

4. The convergence of the iterative process within

the time step under consideration is estimated: if the

following inequality holds

max

1en ||εc(i+1)(e)(tm+1)−εc(i)(e)(tm+1)|| < δ, (26)

where δis the speciﬁed accuracy, then there is a transi-

tion to the next time step (if the speciﬁed time interval

has not yet been passed (tm+1 < tend)), otherwise —

to the next iteration.

To automatically select the length of the time step

at time tm+1, we will use the following algorithm

[11], based on Richardson extrapolation:

τnew =τ·min 

f1,max 

f2, f3·tol

err 1

p+1 



,

(27)

where f3is the guarantee coefﬁcient, f1,f2are the co-

efﬁcients that ensure not too fast increase or decrease

of the step, in this case p= 1 is the order of accuracy

of the explicit and implicit Euler method, tol is the re-

quired level of accuracy. The implicit Euler method is

A-stable, so when choosing the step length, one only

needs to follow the approximation of the data [11].

The estimation of the local relative error for the

moment of time tm+1 is carried out as follows:

err =1

2p−1max

j



ˆεc

j−εc

εc

j



2

(28)

where ˆεc

j(tm+1)is the component of the global creep

strain vector obtained after from the time point tm−1

the calculation was performed for one time step of

length 2τ,εj(tm+1)is the component of the global

creep strain vector obtained after from the point tm−1

the calculation is performed for two time steps of

length τ. If err 6tol, then the resulting values

are εj(tm+1)are accepted and a new calculation cy-

cle is performed from the point tm+1: 2 time steps

of length τnew and one time step of length 2τnew. If

er > tool, then the resulting values are εj(tm+1)are

considered rejected, there is a return to the point tm−1

and a new calculation cycle is performed, taking into

account that τnew < τ.

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

3 Results of the numerical solution

We will carry out a series of calculations for the fuel

element section, which includes from 1 to 100 fuel

pellets. The inner radius of the pellets is 0.8 mm, the

outer radius is — 3.88 mm, the height is — 10 mm,

the outer radius of the cladding is 4.55 m [12]. We

will consider inconsistent grids: pellets are divided

into 10 elements in the directions of rand z, and the

cladding — into 5 elements in the direction of r, in

the direction of z, each section of the cladding corre-

sponding to the height of one pellet is divided into 10

elements. The average fuel element usage time in the

reactor is several years, so let’s choose the calculation

time tend = 3.2·108c (≈10 years).

Let us consider the creep models used for the ma-

terials of fuel pellets and claddings, which are ele-

ments of the fuel element. The pellets are made of ura-

nium dioxide, the cladding is made of zirconium alloy.

Let’s give these models in more detail, described in

[13].

The rate of stationary creep deformation for fuel

pellets is determined by the expression:

˙εc

i=(A1+A2F)σe−Q1

(A3+D)G2+

+(A4+A7F)σ4,5e−Q2

A5+D+A6σF e−Q3

RT ,(29)

where ˙εc

iis the rate of stationary creep deformation

(1/s), Fis the density of divisions ((divisions/m3)/c),

σis the stress (Pa), R— universal gas constant

(Cal/(mol K)), Tis temperature (K), Dis density (per-

centage of theoretical density), Q1,Q2,Q3is acti-

vation energy (Cal/mol), A1, . . . , A7are the constant

coefﬁcients.

For the cladding, the rate of steady-state creep de-

formation has the form:

˙εc

i=(B1+B2F)σ

(A3+D)G2e−Q4

RT +B3(1−D)+B4C+

+(B5+B6F)σ

(A3+D)G2e−Q4

RT +B7(1−D)+B4C,(30)

where Q4,Q5is activation energy (Cal/mol), Cis the

concentration of oxide mixture (weight percentage),

Gis the fuel pellet size (cm), B1, . . . , B7are the con-

stant coefﬁcients.

Graphs of the dependences of displacements and

stresses at the contact boundaries, as well as their two-

dimensional distributions in pellets and claddings are

presented in [12].

Fig. 2 shows two-dimensional radial, axial and

circumferential stress distributions for pellets, and

Fig. 3 shows similar distributions for the cladding

at the ﬁnal moment of time. Fig. 4 shows two-

dimensional distributions of radial, axial and circum-

ferential creep deformations for pellets, and Fig. 5

shows similar distributions for the cladding. Frag-

ments of distributions corresponding to the section

from the 5th to the 9th pellets are shown, and when

constructing deformed bodies, the applied displace-

ments are increased 10 times for pellets and 50 times

for the cladding for greater clarity. It can be seen from

the graphs presented that radial stresses are compres-

sive in almost the entire pellet, the greatest compres-

sive stresses are achieved in the center of the pellets.

The greatest axial and circumferential stresses are ob-

served at the outer surface of the pellets, the smallest

stresses occur at the inner surface.

(a) (b) (c)

Figure 2: Two-dimensional distributions in the nodes

of pellet elements for the moment of time

t= 3,2·108s: a — radial stresses; b — axial

stresses; c — circumferential stresses

(a) (b) (c)

Figure 3: Two-dimensional distributions in the nodes

of the cladding elements for the moment of time

t= 3,2·108s: a — radial stresses; b — axial

stresses; c — circumferential stresses

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

(a) (b) (c)

Figure 4: Two-dimensional distributions in the nodes

of pellet elements for the moment of time

t= 3,2·108s: a — radial creep deformations; b —

axial creep deformations; c — circumferential creep

deformations

(a) (b) (c)

Figure 5: Two-dimensional distributions in the nodes

of the cladding elements for the moment of time

t = 3, 2 · 108 s: a — radial creep deformations; b —

axial creep deformations; c — circumferential creep

deformations

Fig. 6, Fig. 7 shows graphs of the dependences

of the average values of radial stresses on the outer

surface of fuel pellets for the case of 100 pellets

(one point corresponds to one pellet) at various

points in time. The heat dissipation power is

constant (Fig. 6) or varies according to the

sinusoidal law according to the height of the fuel

column (Fig. 7). It can be seen from the above

graphs that radial stresses decrease 5–10 times over

10 years compared to the initial stresses in the

structure. Thus, taking into account creep defor-

mations leads to a noticeable decrease in stress values

in the structure.

In Fig. 8, Fig. 9 graphs of the average values of

radial and axial displacements on the outer surface

of pel-lets are presented (one point corresponds to

one pel-let, the total number is 100 pellets) for cases

of con-stant and sinusoidal heat dissipation power at

various

Figure 6: Dependences of the average values of the

radial stress σr(z)at r= 0.0038 m: constant heat

dissipation power

Figure 7: Dependences of the average values of the

radial stress σr(z)at r= 0.0038 m: sinusoidal heat

dissipation power

points in time. As can be seen from the ﬁgures, there

is a signiﬁcant displacement of the column of fuel pel-

lets relative to the initial position (up to 1.5 cm, which

is comparable to the height of the pellet).

Figure 8: Dependences of the average values of the

axial displacement uz(z)at r= 0.0038 m: constant

heat dissipation power

The Table 1 shows constant and maximum time

steps for a different number of pellets and tol = 10−4,

the Table 2 shows calculation time and number of

steps, and the Table 3 shows errors in calculations

with a variable step relative to calculations with a con-

stant step.

Fig. 10 and Fig.11 show graphs of the dependencies

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

Figure 9: Dependences of the average values of the

axial displacement uz(z)at r= 0.0038 m: sinusoidal

heat dissipation power

Table 1: Constant and maximum time steps

Number of

pellets

Const.

time step

Max.

time step

1 pellet 2·105s2,37 ·107s

2 pellets 5·104s1,09 ·107s

5 pellets 2·104s7,56 ·106s

10 pellets 1,5·103s3,69 ·106s

25 pellets 750 s1,24 ·106s

50 pellets 250 s 5,31 ·105s

100 pellets — 2,25 ·105s

the τstep and estimates of the local relative error err

from time to time for the case of 25 pellets and a given

accuracy tol = 10−4. The ﬁgures show how the time

step τchanges and that err does not exceed the value

of tol throughout the calculation.

Figure 10: Dependencies of the τstep from time for

25 pellets

From the results given in the Table 2, it can be

seen that when using a variable time step, it is pos-

sible to reduce the calculation time by 10–15 times

compared to the calculation with a constant step (10

Figure 11: Dependencies of err from time

for 25 pellets

Table 2: Calculation time and number of time steps

Number of

pellets

Const.

time step

Max.

time step

1 pellet 91,4 s

1600 steps

10,9 s

30 steps

2 pellets 122,6 s

6400 steps

23,1 s

52 steps

5 pellets 1623,9 s

16000 steps

134,7 s

85 steps

10 pellets 10075,5 s

213333 steps

664,7 s

176 steps

25 pellets — 3327,6 s

450 steps

50 pellets — 9213,2 s

911 steps

100 pellets — 25377,3 s

1882 steps

pellets). For 100 pellets (125 thousand unknowns) and

an accuracy of 10−4, the calculation time is approxi-

mately 7 hours, 1882 time steps were performed.

The displacement error was calculated as follows:

ε=v

(uII

−uI

ri)2+(uII

z−uI

zi)2

(uI

ri)2+(uI

zi)2

,(31)

where the number Idenotes a calculation with a con-

stant step, and the number II denotes a calculation

with a variable step.

For the cases of 25, 50 and 100 pellets, instead

of calculating with a constant step (they were not

completed due to high computational costs and the

inability to estimate the maximum possible step in

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

Table 3: Errors with respect to the calculation with a

constant step

Number of pellets tol = 10−4tol = 10−5

1 pellet 6,72 ·10−32,02 ·10−3

2 pellets 9,21 ·10−33,59 ·10−3

5 pellets 1,97 ·10−31,22 ·10−3

10 pellets 4,90 ·10−22,51 ·10−2

25 pellets 5,15 ·10−22,86 ·10−2

50 pellets 6,42 ·10−23,53 ·10−2

100 pellets 8,36 ·10−24,18 ·10−2

time), a calculation with a variable step was used at

tol = 10−6. When using a variable step, the number

of time steps decreases sharply, and with the increase

in the number of pellets, the errors regarding calcula-

tions with a constant step increase.

4 Conclusion

The formulation of the quasi-stationary problem of

multicontact interaction of a system of axisymmetric

thermoelastic bodies under thermomechanical loading

conditions, taking into account the creep process, is

presented. A numerical algorithm for solving such

problems based on the mortar method is described.

For the numerical solution of the problem simulat-

ing creep processes, an algorithm based on the use of

the implicit Euler method is considered, the Newton

method is used for linearization. The algorithm of au-

tomatic step selection based on Richardson extrapo-

lation was applied, which made it possible to signiﬁ-

cantly increase the time steps and reduce the calcula-

tion time. Other, more complicated algorithms based

on the use of artiﬁcial intelligence or neural networks

can also be applied. The results of applying the in-

troduced algorithm to solve a demonstration problem

simulating some processes in a fuel element section

for a mode with a constant heat dissipation power

are presented. In the future, it is planned to include

other signiﬁcant physical phenomena in the mathe-

matical model that must be taken into account when

fully modeling fuel elements, such as pellet cracking

and plasticity in the cladding.

Acknowledgements: The study was carried out

at the expense of a grant Russian Science Founda-

tion No. 22-21-00260, https://rscf.ru/project/22-21-

00260/.

References:

WSEAS TRANSACTIONS on APPLIED and THEORETICAL MECHANICS

DOI: 10.37394/232011.2022.17.10

Pavel Aronov, Mikhail Galanin, Alexandr Rodin

E-ISSN: 2224-3429

Volume 17, 2022

[1] I.V. Stankevich, M. E. Yakovlev and Si Tu Xtet,

Development of Contact Interaction Algorithm on the

Basis of Schwarz Alternating Method, Herald of the

Bauman Moscow State Technical University, Series

Natural Sciences, Special edition: Applied Mathematics,

2011, pp. 134–141.

[2] I. Babuska, The Finite Element Method with Penalty,

Mathematics of Computation, 27, 1973, pp. 221–228.

[3] B. I. Wohlmuth, A Mortar Finite Element Method

Using Dual Spaces for the Lagrange Multiplier, SIAM

Journal on Numerical analysis 38, 2000, pp. 989–1012.

[4] T. Gustafsson, R. Stenberg, J. Videman, Error

analysis of Nitsche’s mortar method, Numerische

Matematik 142, 2019, pp. 973–994.

[5] V. S. Zarubin, G. N. Kuvyrkin, Mathematical Models

of Continuum Mechanics and Electrodynamics, Bauman

Moscow State Technical University, Moscow, 2008, 512

[6] N. N. Malinin, Applied Theory of Plasticity and

Creep, Mechanical engineering, Moscow, 1975, 400 p.

[7] L. A. Rozin, Variational Problem Statements for

Elastic Systems, Leningrad University, Leningrad, 1978,

222 p.

[8] T. X. Duong, L. De Lorenzis and R. A. Sauer, A

segmentation-free isogeometric extended mortar contact

method, Computational Mechanics 63, 2019, pp. 383–

407.

[9] P. S. Aronov, M. P. Galanin and A. S. Rodin,

Application of the MSSOR for solving systems of linear

equations corresponding to problems with a changing

configuration of the contact surface, AIP Conference

Proceedings 2448, 2021, pp. 1–12.

[10] S. N. Korobejnikov, Nonlinear Deformation of

Solids, SB RAS, Novosibirsk, 2000, 262 p.

[11] E. Hairer, S. P. Norsett and G.Wanner, Solving

Ordinary Differential Equations. Nonstiff Problems, Mir,

Moscow, 1990, 512 p.

[12] P. S. Aronov, M. P. Galanin and A. S. Rodin,

Mathematical Modeling of the Contact Interaction of the

Fuel Element with Creep Using the Mortar-method,

KIAM Preprint.110, 2020, 24 p.

[13] D. L. Hagrman, A Library of Materials Properties

for Use in the analysis of Light Water Reactor Fuel Rod

Behavior, NUREG/CR-6150 TREE-1280, 1993, 445 p.

Contribution of Individual Authors to the

Creation of a Scientific Article (Ghostwriting

Policy)

The authors equally contributed in the present

research, at all stages from the formulation of the

problem to the final findings and solution.

Sources of Funding for Research Presented in a Scientific

Article or Scientific Article Itself

Funding was received for conducting this study by the Russian

Science Foundation No. 22-21-00260.

Creative Commons Attribution License 4.0 (Attri-bution

4.0 International, CC BY 4.0)

This article is published under the terms of the Creative

Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en _US

Conflicts of Interest

The authors have no conflicts of interest to

declare .