Contents
Introduction..................................................................................................................................... 4
1 Ice mechanics research review .................................................................................................... 7
1.1 Molecular structure.......................................................................................................... 7
1.2 Physical properties of ice .............................................................................................. 10
1.2.1 Ice thickness ...................................................................................................... 10
1.2.2 Phase composition ............................................................................................ 11
1.2.3 Salinity .............................................................................................................. 16
1.2.4 Ice density ......................................................................................................... 18
1.3 Mechanical properties of ice ......................................................................................... 19
1.3.1 Elastic properties ............................................................................................... 19
1.3.2 Compressive strength ........................................................................................ 24
1.3.3 Tensile and shear strengths ............................................................................... 30
1.4 Ice-structure interaction mechanism in case of vertical structures................................ 32
1.5 Researches in the field of numerical modeling of ice problems ................................... 34
1.6 Section conclusions ....................................................................................................... 52
2 Formulation of the problem and methodology of numerical modeling of ice impacts ............. 56
2.1 Formulation of the problem ........................................................................................... 56
2.2 Methodology of numerical simulation of ice impacts ................................................... 57
2.2.1 Numerical method selection ............................................................................. 57
2.2.2 Integration method selection ............................................................................. 61
2.2.3 General modeling techniques............................................................................ 64
2.2.4 Selection of constitutive models for ice material .............................................. 68
2.3 Section conclusions ....................................................................................................... 75
3 Numerical experiments on the impacts of ice formations on vertical structures and development
of the model ............................................................................................................................................ 77
3.1 Assessment of the influence of mesh pattern ................................................................ 77
3.2 Element size influence................................................................................................... 89
3.3 Section conclusions ....................................................................................................... 92
4 Verification of the methodology for modeling ice impacts on vertical offshore structures ...... 94
4.1 Comparison of load with analytical methods of national codes .................................... 94
4.2 Modeling of uniaxial compression tests ........................................................................ 97
4.3 Section conclusions ....................................................................................................... 99
Conclusion .................................................................................................................................. 101
References ................................................................................................................................... 103
Introduction
Timeliness of the topic. The fuel and energy complex is an important element of the Russian and
world economy. One of the main trends of its development at present is the development of hydrocarbon
deposits on the continental shelf. Based on various data, we can say that in the world today, up to 35%
of oil and gas is produced on the shelf and coastal waters.
The main hydrocarbon reserves of the Russian Federation are concentrated in the regions of the
Far North, which are characterized by difficult climatic conditions. At the same time, most of these
reserves belong to deposits located on the shelf of the Arctic Ocean. The development of these fields
requires the creation of special technical means, due to the difficult ice regime, which is common for the
northern seas.
Hydrocarbon production on the Arctic shelf currently involves the construction of unique marine
structures – offshore ice-resistant platforms, which are characterized by increased responsibility. These
platforms are exposed to various environmental factors. An important role in the complex of external
forces acting on the structure is played by loads from drifting ice fields. This requires a detailed
evaluation of the ice load in the design of such unique structures as offshore platforms. Currently, models
of projected platforms on a smaller scale, as well as special experimental devices are being created to
analyze the impact of ice formations, which affects the timing and cost of construction.
At the moment, the use of digital technologies is widespread in the world, including in the field of
construction. With the help of mathematical and computer modeling many problems are being solved,
such as creating of numerical models of structures, modeling of physical processes (heat transfer,
deformation, corrosion, cracking, crushing), describing of various material properties (plasticity,
viscoelasticity, fatigue, creep) and many others. One of the important questions is the modeling of loads
and impacts on buildings and structures.
Methods for evaluation of ice loads acting on hydraulic structures are described in modern
domestic and foreign standards. However, these methods have a several disadvantages, such as the lack
of a qualitative description of the anisotropic properties of ice, nature of the ice destruction in the contact
zone, influence of cracking in ice formation and so on. The influence of the geometric dimensions of the
structure on the process of its interaction with the ice cover is also little studied.
In this study will be considered the numerical simulation of the ice formations destruction process
during the impact on the vertical offshore structures. The use of modern methods of numerical modeling
can contribute to a more accurate evaluation of the impact of ice on offshore structures, which will
positively affect the analysis of the stress-strain state of structures and their design.
Nowadays state of problem. The question of the influence of the ice regime of freezing water
areas on structures has been actively discussed and considered since the second half of the 20th century.
4
Various conferences are regularly held and many articles on ice issues are being published. Increased
attention to this topic is largely associated with the development of fuel and energy complex in the world
and with the increasing rate of industrial development of oil and gas deposits on the Arctic shelf.
However, the question of exactly numerical simulation of processes of ice interaction and
structures on the shelf is appeared relatively recently, due to the development of digital technologies in
the world. Nowadays there are a number of studies on the impacts of ice formations on structures using
numerical simulation. But these studies do not sufficiently describe the properties of ice and ice
formations, the processes occurring at the contact of interaction, the influence of the geometric
dimensions of structures on the nature of the destruction of the ice and other important aspects that
should be considered in the simulation.
Goals and objectives of the study. The main goal of research – to develop a methodology for
numerical modeling of ice load formation processes on the vertical supports of offshore structures and
to evaluate its applicability to solving real engineering tasks.
In this paper, the objectives are as follows:
- to study the properties of ice as a material, to determine the degree of knowledge of the issue of
ice-structure interaction modeling;
- to study the mechanism of interaction between ice and offshore structures;
- select and apply constitutive material models, taking into account the various physical properties
of ice and fracture mechanism;
- to develop an algorithm for modeling the interaction of ice formation and vertical structures;
- to conduct numerical experiments using developed model in order to study and improve it;
- to verify the developed model, compare the results of numerical simulation with experimental
data.
The object of the research is the process of interaction process between ice formations and offshore
structures. The subject of research is numerical modeling of ice load formation processes.
Academic novelty. The academic novelty of the study is provided by the new formulation of the
problem of numerical simulation of the ice loads formation processes. Methodology of modeling of ice
impacts is developed in detail in order to take into account the many factors and features of ice-structure
interaction process that were not previously considered. The modeling of fracture of ice formations is
conducted considering various physical and mechanical properties of ice.
Theoretic and practical relevance. Theoretical significance is to consider the features of
numerical modeling of ice impacts, describing the approaches, advantages and disadvantages of the
methods used.
The practical significance is that numerical model being developed here as well as methods for the
numerical simulation of ice-structure interaction can contribute to an accurate assessment of the
5
influence of the ice regime in the analysis of the stress-strain state and the design of offshore ice-resistant
platforms and other structures, that will positively affect the cost and speed of design. The methods
presented in this thesis, will contribute to improvement of methods for calculating ice loads given in
modern codes.
Global problem of this study is the production of hydrocarbons in the Arctic region and its
development. The local problem of the research is the assessment of the impact of ice formations on the
shelf structures.
Methodology and research methods. In this study, numerical experiments are performed in the
SIMULIA Abaqus finite element software package. Final methodology and numerical model are verified
using by comparing the results with experimental data.
The following research methods were used in this work: review and analysis of existing researches
on the topic, deductive decomposition of the problem of numerical simulation of interaction, synthesis
of the general methodology for performing numerical experiments. axiomatic method for creating initial
numerical models, qualitative and mathematical comparison of the results.
The main states for defense:
- statistical relation of ice properties;
- constitutive models used for modeling of ice material;
- methodology for creating of numerical model;
- numerical model parameters and features.
Approbation. The main results of the master thesis were published in two papers in School of
Engineering Bulletin of Far Eastern Federal University and presented at the seminars and meetings of
the Offshore and Coastal Engineering department.
Also, the thesis was presented at seminars and scientific meeting at Kindai University and Toyo
Construction Company (Osaka, Japan).
The thesis structure is 109 pages, 4 chapters, 92 pictures and 20 tables.
6
1 Ice mechanics research review
1.1 Molecular structure
There are several ice types with different microstructures. But almost all of ice in the world is
ice Ih. Ice Ih has a hexagonal structure. When water freezes hydrogen bonds are formed between its
oxygen atoms. In solid Ih ice crystal the oxygen atoms of the nearest molecules line up in a repeating
unit cell of a tetrahedral shape – four nearby molecules are located at the vertices of a three-sided
pyramid, in the center of which is the fifth water molecule, as shown in the Figure 1.1
Figure 1.1 – Ice lattice structure
Another feature of the molecular structure of ice is the fact that the molecules are arranged in the
order of hexagonal symmetry, in consequence of which the structure of ice is called hexagonal. At the
same time, the hexagonal structure of ice consists of layers parallel to each other, called basal planes.
All basal planes are oriented perpendicular to the main crystallographic axis, called the c-axis [77].
Figure 1.2 shows the projections of the three-dimensional crystal lattice of ice along two main directions.
Figure 1.2 – Idealized arrangement of atoms in ice, where oxygen atoms are shown black, and
hydrogen atoms white: a) view of the crystal lattice along the C-axis; b) view of the crystal lattice
along the basal plane [64]
7
To analyze the effects of ice formations on structures, it is important to take into account the
microstructure of ice formations. The formation of natural ice features is closely related to the physical
processes occurring during freezing. During ice solidification in nature, various natural factors influence
it, such as snow, wind, currents, temperature changes, etc. In addition, air bubbles, salts in sea water,
alternate freezing and melting affect the microstructure of ice formations. As a result, the ice body has a
non-uniform microstructure.
The process of solidification of floating ice sheets initiates on or near the surface, and then
continues unidirectionally through the downward movement of the ice–water interface as heat is
transferred from the relatively warm water below to the colder atmosphere above [66]. The thicker the
ice, the slower it grows, since ice prevents the penetration of cold temperatures from the air. The top
layer of ice usually consists of crystals of a variety of shapes. The most common are needle-shaped
crystals with c-axes inclined in many directions, but granular ice can be observed. This layer is a small
part of the total thickness of thick ice sheets. As the ice–water interface advances, favorably oriented
crystals expand at the expense of less favorably oriented ones which are wedged out as presented in
Figure 1.3. The plate-like structure of the crystals denotes the orientation of the basal planes. This
process is termed geometric selection and accounts for the development of growth texture.
Figure 1.3 – A schematic sketch showing the process of geometrical selection during the unidirectional
solidification of water [66]
Thus, when sheets grow from thin layers in which crystals of all orientations are present, they are
expected to possess predominantly the c-axis horizontal texture. Sheets of first-year sea ice, for instance,
are generally composed of columnar-shaped grains, elongated in the growth direction [66] as presented
in Figure 1.4 and Figure 1.5. In all sea ice formed by unidirectional freezing, c-axis-horizontal
orientations develop rapidly after an initial ice skim forms, and they dominate the rest of the ice growth
[80]. Thus, the c-axes of crystals that make up sea ice, for instance, are more or less horizontally oriented
(i.e., perpendicular to the growth direction), and are occasionally aligned within this plane [66].
Also, during freezing process, the brine movement can be observed. As ice grows and the ice–
water interface advances downwards into the melt, salt ions are rejected from the ice. The salt builds up
ahead of the advancing interface, increasing the salinity of a thin layer of a few millimeters to a few
8
centimeters in thickness. The resulting gradient in salt concentration leads to diffusion of salt away from
the interface towards the less saline ocean [72]. Within sea ice, a downward flow increases the local bulk
salinity and porosity and ultimately leads to the formation of brine channels: Locally, brine is replaced
by more saline brine from above. Since the more saline brine is superheated (i.e. above the freezing
point) with respect to the temperature of the surrounding ice, it will dissolve surrounding ice partially to
attain thermodynamic equilibrium. These channels can be thought of as tubular river systems in which
the tributaries are arranged with cylindrical symmetry around the main drainage channels [67].
Figure 1.4 – Schematic summarizing the main ice textures, growth conditions and timescales, and
typical winter temperature and salinity profiles for first-year sea ice [72]
Figure 1.5 – Schematic drawing showing the heterogeneity of the structure of FY sea ice [67]
9
1.2 Physical properties of ice
After the appearance of the task of developing offshore oil and gas fields in the Arctic zone, where
the ice-free period could be only a few months, active research began in the areas of offshore
construction and ice engineering. Today, there are a lot of research on various characteristics of the ice.
However, due to the fact that ice, especially sea ice, is a rather heterogeneous and anisotropic material,
the research results vary widely. The anisotropy of ice is primarily due to its molecular structure.
The determination of the characteristics of ice has been carried out by many researchers since the
second half of the 20th century. The main characteristics of ice, which have been studied most widely,
are the elastic modulus, Poisson's ratio, density, compressive strength, tensile strength and bending
strength. Friction coefficients of ice on various materials, its salinity, porosity, and various dependencies
between the listed characteristics are also studied.
The following authors made a great contribution to the study of ice properties: Weeks W. F.,
Assur A., Frankenstein G., Garner R., Cox G. F. N., Kovacs A., Timco G. W., Frederking R. M. W.,
Sodhi D. S., Saeki H. A review of the data and many other significant studies of sea ice properties has
been performed by Timco and Weeks (2010) [74].
Also, quite extensive investigations of the properties of both fresh and sea ice have been carried
out by many Russian authors. In this thesis, will be used the studies of Ryvlin, Nazintsev Yu. L.,
Panov V. V. Dmitrash J. A., Moiseev V. I. and Doronin Yu. P.
1.2.1 Ice thickness
Ice thickness is one of the most important parameters of ice fields. Proper assessment of ice
thickness directly affects the accuracy of the ice load. Often, the assessment of the maximum thickness
of ice is made from the condition of the energy balance at the ice-water interface [74]:
𝜑𝑖 𝑑𝑡 = 𝜌𝑖 𝐿𝑑ℎ𝑖 ,
(1.1)
where 𝜑𝑖 – heat flow from ice to air, W/m2;
𝜌𝑖 – ice density, kg/m3;
𝐿 – latent heat of fusion of ice, J/kg
Thus, if we make the assumption that the temperature of the upper surface of the ice is equal to
the air temperature, and the temperature gradient along the vertical is equal to the temperature difference
between the lower and upper surfaces, the ice thickness ℎ𝑖 , m, can be determined by the equation [19,
20, 74]
ℎ𝑖 =
𝜆𝑖 (𝑇𝑏 − 𝑇𝑎 )
,
𝜑𝑖
10
(1.2)
where 𝜆𝑖 – thermal conductivity of ice, W/(m∙°С)
𝑇𝑏 and 𝑇𝑎 – temperature on the lower surface of the ice and air temperature, respectively, °С.
Integrating equation (1.1) with expressing 𝜑𝑖 from equation (1.2), gives the following equation for
ice thickness ℎ𝑖 [19, 20]:
1
2
𝑡
ℎ𝑖 (𝑡) = ℎ0 (1 +
(𝑇𝑏 − 𝑇𝑎 )
2
∫
𝜆
𝑑𝑡) ,
𝑖
𝐿𝜌𝑖
ℎ02
(1.3)
0
where ℎ0 – ice thickness at 𝑡 = 0.
Assuming that ℎ0 = 0 and ice properties (𝜆𝑖 , 𝐿, 𝜌𝑖 ) are unchanged, we can get the Stefan equation
to determine the thickness of the ice:
1
2
2𝜆𝑖
(𝑇𝑏 − 𝑇𝑎 )𝑡)
ℎ𝑖 = (
𝐿𝜌𝑖
(1.4)
In this case, the integral of the temperature difference in the equation (1.3) is often called the sum
of the number of freezing degree days, the approximate value of which is used in practice in
equation (1.4) and is determined by the following expression
𝑡
(1.5)
∫(𝑇𝑏 − 𝑇𝑎 )𝑑𝑡 ≈ ∑(𝑇𝑏 − 𝑇𝑎 ) 𝑡
0
The strong influence on the growth of ice thickness has snow. Assuming that the temperature of
the upper surface of the snow cover is equal to the air temperature 𝑇𝑎 , the formula (1.3), when accounting
for snow, takes the following form [19, 20]:
𝑡
ℎ𝑖 (𝑡) = ℎ0 (1 +
1
2
2
𝜆𝑖 (𝑇𝑏 − 𝑇𝑎 )
𝑑𝑡) ,
2∫
ℎ0 𝐿𝜌 (1 + 𝜆𝑖 ℎ𝑠 )
𝑖
0
𝜆𝑖 ℎ
(1.6)
The above equations usually show an overestimated value of ice thickness, since they do not take
into account the effect of wind on heat transfer between the ice surface and air, as well as the effect of
heat flow at the ice-water interface.
1.2.2 Phase composition
The phase composition of sea ice has a very large impact on almost all its properties. In sea ice, in
addition to crystals, there are a liquid phase, air and salts in a dissolved and solid state. The mass 𝑚, g,
of a certain volume of ice is determined by the following equation [19, 20, 51]
𝑚 = 𝑚𝑖 + 𝑚𝑤 + 𝑚𝑠 ,
where 𝑚𝑖 – mass of freshwater ice, g;
11
(1.7)
𝑚𝑤 – the mass of water in the liquid state, g;
𝑚𝑠 – mass of salts, g:
𝑚𝑠 = 𝑚𝑏𝑠 + 𝑚𝑐𝑠 ,
(1.8)
𝑚𝑏𝑠 – the mass of salts in the dissolved state (in brine), g;
𝑚𝑐𝑠 – mass of salts in the crystalline state, g;
Salinity of ice is the ratio of the mass of salts to the total mass of ice. Thus, the salinity of ice 𝑠𝑖 ,
fraction by weight, is determined by the equation
𝑠𝑖 =
𝑚𝑠
,
𝑚
Brine salinity 𝑠𝑏 , fraction by weight, is the mass of dissolved salts per unit mass of brine
𝑚𝑏𝑠
𝑠𝑏 =
𝑚𝑏
(1.9)
(1.10)
where 𝑚𝑏 – brine weight, g:
𝑚𝑏 = 𝑚𝑏𝑠 + 𝑚𝑤
(1.11)
Doronin Yu. P. on the basis of research Schwerdtfeger (1965) received the following dependence
of brine salinity on temperature:
𝑑𝑠𝑏
𝛽
=
,
𝑑𝑇 1 + 𝑚𝑏𝑠
𝑚𝑤
(1.12)
where 𝛽 – empirical coefficient of proportionality between brine concentration and temperature, °C-1.
Mass of brine 𝑚𝑏 , g, Doronin Yu. P. defines by the following expression
𝑚 ∙ 𝑠𝑖 − 𝑚𝑐𝑠
𝑚𝑏 = 𝑚𝑏𝑠 + 𝑚𝑤 =
,
𝑠𝑏
(1.13)
Nazintsev and Panov in their book [51] on the basis of the experiments of Nelson and Thompson
(1954) and Gitterman (1937) determined the values of 𝛽 and 𝑚𝑏𝑠 /𝑚𝑤 , used in the equation (1.12):
−1.848 ∙ 10−2 , at 0 ≥ 𝑇 ≥ −7.5 °𝐶
𝛽 = {−1.077 ∙ 10−2 , at − 7.5 ≥ 𝑇 ≥ −22.4 °𝐶
−0.532 ∙ 10−2 , at − 22.4 ≥ 𝑇 ≥ −30.0 °𝐶
(1.14)
𝛽𝑇 ∙ 10−2 , at 0 ≥ 𝑇 ≥ −7.5 °𝐶
𝑚𝑏𝑠
= { (5.670 − 1.077 ∙ 𝑇) ∙ 10−2 , at − 7.5 ≥ 𝑇 ≥ −22.4 °𝐶
𝑚𝑤
(16.770 − 1.532 ∙ 𝑇) ∙ 10−2 , at − 22.4 ≥ 𝑇 ≥ −30.0 °𝐶
(1.15)
Thus, the dependence of the ratio of the mass of salts to the mass of water in brine is represented
by a piecewise linear function of temperature as shown in Figure 1.6. In general, Nazintsev Yu. L.
identified the following equations [51]:
1 + (𝛽𝑇 + 𝑙)
,
𝛼𝑇 + 𝑘
𝛽𝑇 + 𝑙
𝑚𝑠𝑡 = 𝑠𝑖 (1 + 𝛾) (1 −
),
𝛼𝑇 + 𝑘
𝑚𝑏 = 𝑠𝑖
12
(1.16)
(1.17)
𝑠𝑏 =
𝛽𝑇 + 𝑙
,
1 + 𝛽𝑇 + 𝑙
(1.18)
where 𝑚𝑠𝑡 – mass of salts in brine (taking into account salts, including crystallized water), g;
𝛼, 𝛽, 𝑘 и 𝑙 – empirical coefficients presented in the Table 1.1.
Figure 1.6 – Chart of equilibrium salt concentrations in sea ice [51]
Table 1.1 – The values of empirical coefficients according to [51]
Ice temperature 𝑇𝑖 , °C
Coefficient
𝛼, °K-1
𝛽, °K-1
𝑘
𝑙
0 ≥ 𝑇 ≥ −7.5
-0.01848
-0.01848
0
0
−7.5 ≥ 𝑇 ≥ −22.4
-0.01367
-0.01077
0.0332
0.0567
−22.4 ≥ 𝑇 ≥ −30.0
-0.10567
-0.00532
-1.9920
0.1677
Notation: the boundaries of temperature ranges are determined by the points of deposition of salts
Na2SO4∙10H2O (at -7,5 °С) and NaCl∙2H20 (at -22,4 °С) [51].
Brine volume 𝜈𝑏 , ‰, is also used in practice. Nazintsev Yu. L. suggests the following formula for
determining the brine volume
𝜈𝑏 =
𝑚𝑏 𝜌𝑠𝑖
,
𝜌𝑏
(1.19)
where 𝜌𝑠𝑖 – sea ice density g/cm3;
𝜌𝑏 – brine density g/cm3.
USA researchers Frankenstein and Garner (1967) [21] propose to determine 𝜈𝑏 by the equation
13
𝜈𝑏 = 𝑠𝑖 (
𝑎
+ 𝑏) ,
|𝑇𝑖 |
(1.20)
where 𝑎 and 𝑏 – empirical coefficients presented in the Table 1.2.
Table 1.2 – Empirical coefficients 𝑎 and 𝑏 according to Frankenstein and Garner [21]
Coefficient
Ice temperature 𝑇𝑖 , °C
𝑎
𝑏
−0.5 ≥ 𝑇 ≥ −2.06
52.56
-2.28
−2.06 > 𝑇 ≥ −8.2
45.917
0.93
−8.2 > 𝑇 ≥ −22.9
43.795
1.189
The equation (1.20) is widely used by many authors, as well as regulatory documents. For cases
where a simple but less accurate equation is required Frankenstein and Garner determine the brine
volume at temperatures from 0.5 to 20 °C, taking 𝑎 = 49.185 and 𝑏 = 0.532.
Equations for determining the total phase composition of ice have been also described by Cox and
Weeks (1983) [14]. Based on previous studies by different authors, Cox and Weeks suggested using
some approximating functions to determine the volume fractions of brine, air, and solid salts in sea ice.
The volume fractions of brine 𝜈𝑏 , air 𝜈𝑎 and solid salts 𝜈𝑠𝑠 according to Cox and Weeks (1983) are
respectively defined as follows:
𝜈𝑏 =
𝜌𝑠𝑖 𝑠𝑖
,
𝐹1 (𝑇)
𝜌𝑠𝑖
𝐹2 (𝑇)
+ 𝜌𝑠𝑖 𝑠𝑖
,
𝜌𝑖
𝐹1 (𝑇)
𝜌𝑏
𝜈𝑠𝑠 = 𝐶
+ 𝜈𝑏 ,
𝜌𝑠𝑠
𝜈𝑎 = 1 −
(1.21)
(1.22)
(1.23)
where 𝐹1 (𝑇), 𝐹2 (𝑇), 𝐶 – empirical functions of temperature:
𝐹1,2 (𝑇) = 𝛼0 + 𝛼1 𝑇 + 𝛼2 𝑇 2 + 𝛼3 𝑇 3 ,
𝜌𝑏
+1
𝜌𝑖
,
𝜌𝑏 𝜌𝑏
−
𝜌𝑖 𝜌𝑠𝑠
𝐹2 (𝑇) −
𝐶=
(1.24)
(1.25)
where 𝛼0 , 𝛼1 , 𝛼2 , 𝛼3 – empirical coefficients presented in Table 1.3;
𝜌𝑖 – density of pure freshwater ice, t/m3;
𝜌𝑠𝑠 – the average solid salt density assumed to be constant at 1.5 t/m3 [14].
Nazintsev Yu. L and Panov V. V determine the density of brine 𝜌𝑏 , t/m3, using the expression
obtained by Cox and Weeks (1975) [51]:
𝜌𝑏 = 1 + 0.0008𝑠𝑏
14
(1.26)
In equation (1.26) the brine salinity 𝑠𝑏 is taken in ‰.
Table 1.3 – Empirical coefficients 𝛼𝑖
Coefficient
Ice temperature 𝑇𝑖 , °C
−2 ≥ 𝑇 ≥ −22.9
−22.9 > 𝑇 ≥ −30
𝛼0
𝛼1
𝛼2
𝛼3
𝐹1
-4.732
-22.45
-0.6397
-0.01074
𝐹2
8.903 10-2
-1.763 10-2
-5.33 10-4
-8.801 10-6
𝐹1
9899
1309
55.27
0.716
𝐹2
8.547
1.089
4.518 10-2
5.819 10-4
Kovacs (1997) [39] derived following equation based on ice core measurements made at 10- cm
increments through 44 Beaufort Sea winter ice floes shown in the Figure 1.7:
𝜈𝑏 = 41,64𝑠𝑖0,88 |𝑇|−0,67 ,
(1.27)
Figure 1.7 – Sea ice bulk brine volume versus average ice floe temperature and bulk salinity according
to Kovacs (1997) [39]
Since in the sea ice there are not only inclusions of brine but also air inclusions, the total porosity
of sea ice 𝜈𝑏 , volume fraction, is determined by the equation [74]
𝜈𝑡 = 𝜈𝑏 + 𝜈𝑎 ,
(1.28)
where 𝜈𝑏 – relative brine volume;
𝜈𝑎 – the relative air volume;
The volume of air in sea ice as well as the methods for determining it are less studied than the
volume of brine. In practice, usually the volume fractions of air and solid salts are neglected due to their
small value. The equation for determining the total porosity of sea ice 𝜈𝑡 , ‰, was obtained by
Kovacs (1997) [39]:
𝜈𝑡 = 19.37 + 36.18𝑠𝑖0.91 |𝑇|−0.69 ,
15
(1.29)
1.2.3 Salinity
Sea ice salinity has also been studied for quite a long time. Typically for measuring salinity of sea
ice cores are taken. Then the cores are cut into separate plates, placed in separate sealed containers and
left until complete melting. After special device (salimeter) is used to measure the electrical conductivity
of melt water. On the basis of electrical conductivity and temperature values the salinity is determined
[74].
Ryvlin (1974) [61] was one of the first who empirically expressed the salinity of ice 𝑠𝑖 , ‰, vs. ice
floe thickness. His equation takes into account the seawater salinity and the ice growth rate:
0.5
𝑠𝑖 = 𝑠𝑤 (1 − 𝑠𝑟 )𝑒 −𝑎∙ℎ𝑖 + 𝑠𝑤 𝑠𝑟 ,
(1.30)
where 𝑠𝑤 – salinity of water, ‰;
ℎ𝑖 – ice thickness, m;
𝑎 – growth rate coefficient;
𝑠𝑟 – salinity ratio:
𝑠𝑟 =
𝑠𝑖𝑓
,
𝑠𝑤
(1.31)
here 𝑠𝑖𝑓 – final bulk salinity at end of growth season, ‰.
Ryvlin A. Ya proposed to use the following values of the parameter 𝑎:
𝑎={
0.35 at 𝐺𝑅 ≥ 4 cm/day
0.5 at 𝐺𝑅 ≤ 0.5 cm/day
(1.32)
where 𝐺𝑅 – ice sheet growth rate, cm/day.
Ryvlin also suggests that where 𝐺𝑅 is not known, one may assume the value of 0.5 and 𝑎 ≈ 0.13.
At the same time, the salinity of the ice was investigated by American researchers Cox and Weeks
(1974) [15]. The authors analyzed the results of measuring the salinity of ice accumulated over several
years by many authors as it shown at Figure 1.8 and obtained the following dependences of the salinity
𝑠𝑖 , ‰, of ice on its thickness:
14.237 − 19.39ℎ𝑖 at ℎ𝑖 < 0.4 m
𝑠𝑖 = {
7.88 − 1.59ℎ𝑖 at ℎ𝑖 ≥ 0.4 m
(1.33)
Large and significant ice studies were performed by Timco and Frederking (1990) [75]. In this
paper, the authors analyzed more than 400 tests of small samples and obtained empirical equations for
some physical and mechanical properties of ice. As for the salinity of the ice, they used additional data
together with the data shown in Figure 1.8 and obtained modified equations:
𝑠𝑖 = {
13.4 − 17.4ℎ𝑖 at ℎ𝑖 ≤ 0.34 m
8 − 1.62ℎ𝑖 at ℎ𝑖 > 0.34 m
16
(1.34)
Figure 1.8 – Sea-ice bulk salinity vs. floe thickness according to Cox and Weeks (1974) [15]
Some empirical dependences of ice salinity were obtained by Kovacs (1996) [40] based on the
experiments of Fedotov (1973), Cox and Weeks (1974), Ryvlin (1974) and others. He combined all
Arctic and Antarctic first-year sea-ice bulk salinity vs. floe thickness data compiled from numerous
sources as it shown on Figure 1.9.
Figure 1.9 – All Arctic and Antarctic sea-ice bulk salinity vs. floe thickness data compiled from
numerous sources according to Kovacs (1996) [40]: a) for first-year ice; b) for multi-year ice
For one-year ice, when ice thickness is less than 2 m, the salinity 𝑠𝑖 , ‰, according to
Kovacs (1996) [40] is determined by equation
𝑠𝑖 = 4.606 + (
91.603
),
ℎ𝑖
(1.35)
where ℎ𝑖 – ice thickness, sm.
For multi-year ice the salinity 𝑠𝑖 , ‰, when ice thickness is from 2 to 9 m, is determined by equation
𝑠𝑖 = 1.85 + (
17
80217.9
)
ℎ𝑖2
(1.36)
It should be noted that the square value of the correlation coefficient 𝑟 2 obtained by Kovacs for
the regression equation (1.35) is 0.730, and for regression equation (1.36) it is equal to 0.222.
1.2.4 Ice density
Ice density is an important property of ice. The density of ice directly affects the buoyancy of ice
formations, as well as the load from its own ice weight. The density of ice was studied widely by many
authors. Reported values vary over a wide range from 0.72 t/m3 to 0,94 t/m3 [74].
Nazintsev Yu. L and Panov V. V [51] determine the density of pure freshwater ice 𝜌𝑖 (𝑇), t/m3, by
the equation
𝜌𝑖 =
𝜌0
,
1 + 𝛾𝑇
(1.37)
where 𝜌0 – pure freshwater ice density at 0 °C equal to 916.8 kg/m3;
𝛾 – bulk ice thermal expansion coefficient equal to 1.58∙10-4 K-1.
Nazintsev Yu. L [51] also performed calculations of the phase composition of sea ice using
equations (1.16)-(1.19). The density of sea ice 𝜌𝑠𝑖 , t/m3, he determined by the equation
𝜌0
𝜌𝑠𝑖 =
1 + 1.65 ∙ 10−4 ∙ 𝑇
(1.38)
Cox and Weeks (1983) [14] used the following equation, that was calculated by Pounder (1965),
to determine the density of pure freshwater ice 𝜌𝑖 , t/m3, which is directly or indirectly used in equations
(1.21)-(1.25):
𝜌𝑖 = 0.917 − 1.403 ∙ 10−4 𝑇
(1.39)
Authors also suggests the expression for determining the density of sea ice 𝜌𝑠𝑖 , t/m3, taking into
account its gas and brine volumes [14]:
𝜌𝑠𝑖 = (1 − 𝜈𝑎 )
𝜌𝑖 ∙ 𝐹1 (𝑇)
,
𝐹1 (𝑇) − 𝜌𝑖 ∙ 𝑠𝑖 ∙ 𝐹2 (𝑇)
(1.40)
where 𝜈𝑎 – relative air volume, volume fraction;
𝐹1 (𝑇), 𝐹2 (𝑇) – empirical functions of temperature determined by equation (1.24).
If we accept that there are no air inclusions in the ice, i.e. 𝜈𝑎 = 0, the density of sea ice can be
easily determined.
Some empirical equation for determining the density of ice monocrystals was used by Gammon et
al. (1983) [25] in series of experiments for obtaining elastic properties of ice:
𝜌0
𝜌𝑖 =
1 + 1.576 ∙ 10−4 𝑇 − 2.778 ∙ 10−7 𝑇 2 + 8.850 ∙ 10−9 𝑇 3 − 1.778 ∙ 10−10 𝑇 4
According to authors the equation (1.41) gives the maximum standard deviation ±1.5%.
18
(1.41)
1.3 Mechanical properties of ice
1.3.1 Elastic properties
As in the case of most materials, under small loads, deformations are recoverable and proportional
to stresses according to Hooke's law [66]:
𝜎𝑖𝑗 = 𝐶𝑖𝑗𝑘𝑙 𝜀𝑘𝑙 ,
(1.42)
𝐶𝑖𝑗𝑘𝑙 – elastic stiffness constants, Pa;
𝜎𝑖𝑗 and 𝜀𝑘𝑙 – the components of stress and strain tensors 𝑇𝜎 and 𝑇𝜀 respectively:
𝜎11 𝜎12 𝜎13
𝜎1
𝑇𝜎 = [𝜎12 𝜎22 𝜎23 ] = [𝜎6
𝜎13 𝜎23 𝜎33
𝜎5
𝜀11 𝜀12 𝜀13
𝜀1
𝑇𝜀 = [𝜀12 𝜀22 𝜀23 ] = [𝜀6
𝜀13 𝜀23 𝜀33
𝜀5
𝜎6 𝜎5
𝜎2 𝜎4 ],
𝜎4 𝜎3
𝜀6 𝜀5
𝜀2 𝜀4 ]
𝜀4 𝜀3
Ih ice possesses only five independent stiffness constants: 𝐶11 , 𝐶12 , 𝐶13 , 𝐶33 , 𝐶44 which are
parameters of the material that characterize the ability of a material to deform elastically. The stiffness
matrix for monocrystals of ice Ih become [66]:
𝐶11
𝐶12
𝐶13
𝐶𝑖𝑗 =
0
0
( 0
𝐶12
𝐶11
𝐶13
0
0
0
𝐶13
𝐶13
𝐶33
0
0
0
0
0
0
𝐶44
0
0
0
0
0
0
𝐶44
0
0
0
0
0
0
1/2(𝐶11 − 𝐶12 ))
These constants are expressed with respect to a rectangular coordinate system whose axes are
specified in terms of the unit cell of the monocrystal. For ice Ih, axis 𝑋3 is parallel to the c-axis and 𝑋1
and 𝑋2 may be taken as any pair of axes within the basal plane.
The values have been obtained by Gammon et al. (1983) [25]. They applied the method of
Brillouin spectroscopy in which incident laser light was scattered from thermally induced acoustic waves
and Doppler shifted. Table1.4 lists the weighted mean elastic constants calculated from the four sets of
Brillouin measurements, along with several quantities derived from these elastic constants.
The values in Table 1.4 were established using the value of the ice density calculated by the
equation (1.41). The authors also derived an equation for determining the values of stiffness coefficients
at other temperatures. According to value Gammon et al. (1983) the value 𝐶(𝑇) of any of the
fundamental stuffiness constants at temperature 𝑇 may be given by the relationship:
𝐶(𝑇) = 𝐶(𝑇𝑟 ) ∙ (1 − 𝑎(𝑇 − 𝑇𝑟 )),
where 𝑇𝑟 – reference temperature at which the constant was measured;
𝑎 – empirical factor equal to 1.42 ∙ 10−3 K−1.
19
(1.43)
Table 1.4 – Elastic parameters for ice monocrystals at -16 °C [25, 66]
Property
Dynamic elastic moduli
Bulk modulus
Poisson’s ratio
Units
GPa
GPa
-
Symbol
Value
𝐶11
13.929 ± 0.041
𝐶12
7.082 ± 0.039
𝐶13
5.765 ± 0.023
𝐶33
15.010 ± 0.046
𝐶44
3.014 ± 0.011
𝐾
8.899 ± 0.017
𝜇12
0.415
𝜇13
0.224
𝜇31
0.274
Thus, the authors determined the elastic properties of the ice monocrystals. In addition, the desired
values are the so-called "true" values. Since the elastic modulus has a strong dependence on the strain
rate and the time of their action, there is a true elastic modulus and effective elastic modulus.
The true elastic modulus of ice characterizes the deformation of a purely elastic, time-independent
character [33]. Its value is determined using dynamic measurements and therefore it is often called the
dynamic modulus of elasticity. An effective modulus characterizes time-dependent recoverable and
inelastic non-recoverable deformations. The effective modulus is determined by static methods, and
therefore it is sometimes called the static modulus of elasticity. The effective modulus is lower than the
true modulus of elasticity.
Since the problems solved in this work have a scale much larger than the scale of the microstructure
of one ice crystal, the best solution would be to consider the properties of polycrystals of ice Ih. In
principle the elastic properties of polycrystals free from porosity, inclusions and other defects can be
calculated from the fundamental elastic constants and from the orientations, sizes and shapes of the
grains [66]. In this case, two approaches can be considered to determining the elastic characteristics of
ice formation:
- consider that the ice feature has isotropic properties if we accept the fact that ice crystals or grains
are randomly oriented in space throughout the ice body;
- consider the ice growth process described earlier. Since the ice texture develops during growth
and in most cases has more or less established orientation of crystals and c-axes, in this case we can say
that ice formations will have orthotropic properties.
Isotropic elasticity. First, let’s consider the elastic properties of ice features, in the case of isotropy.
Elastic properties in that case are completely described by only two independent constants, chosen from
20
Young’s modulus 𝐸, the shear modulus 𝐺, Poisson’s ratio 𝜈 and the bulk modulus 𝐾 [66]. Any
combination of the two parameters described above can be represented, since they are interdependent
and the equations for their determination are well known from the strength of materials, for example:
𝐸
,
2(1 + 𝜇)
(1.44)
𝐸
,
3(1 − 2𝜇)
(1.45)
𝐺=
𝐾=
It is more common to talk about the most widely used elastic parameters: Poisson's ratio and
Young's modulus. The modulus of ice is a strong function of the loading rate, temperature, grain size,
and ice type. A large number of studies have been performed to determine the modulus of elasticity of
ice.
The true values of the parameters of ice elasticity were also calculated by Gammon et al.
(1983) [25] based on acoustic measurements. These values are presented In Table 1.5.
Table 1.5 – Polycrystalline (isotropic) averaged elastic parameters at temperature -16 °C [25]
Property
Units
Symbol
Value
𝐾
8.899
𝐸
9.332
𝐺
3.521
𝜇
0.325
𝜆
6.551
𝜇
3.521
Bulk modulus
Young's modulus
GPa
Shear modulus
Poisson's ratio
-
Lame constants
GPa
The dynamic modulus of elasticity 𝐸𝑖,𝑑 and static modulus of elasticity 𝐸𝑖,𝑠 , GPa, of freshwater ice
can be determined by the following equations [41]:
𝐸𝑖,𝑑 = 9.21(1 − 0.00146𝑇),
(1.46)
𝐸𝑖,𝑠 = 5.69 − 0.0648𝑇,
(1.47)
where 𝑇 – ice temperature, °C.
For sea ice, the elastic modulus also depends on the brine volume. The dependence of the modulus
of elasticity of sea ice 𝐸𝑠𝑖 , GPa, on the volume of brine was proposed by Weeks and Assur (1967) [81]:
𝐸𝑠𝑖 = 𝐸𝑖 (1 − 𝜈𝑏 )4 ,
(1.48)
where 𝐸𝑖 – modulus of elasticity of freshwater ice, MPa;
𝜈𝑏 – brine volume fraction.
An equation that directly defines the effective modulus of elasticity of sea ice proposed by
21
Vaundrey (1997) [78]. This equation is used in international ISO/FDIS 19906:2019(E) standard:
𝐸𝑠𝑖,𝑠 = 5.316 − 0.436𝜈𝑏0.5 ,
(1.49)
Langleben and Pounder (1963) [42] processed the results of over 300 measurements of elastic
modulus of sea ice from several locations in the Canadian Arctic and Greenland. They obtained a linear
relation between elastic modulus and brine volume with a small scatter of data as presented in Figure
1.10. True elastic modulus of sea ice 𝐸𝑠𝑖,𝑑 , GPa, as function of brine volume according to Langleben and
Pounder (1963) [42] is defined as follows
𝐸𝑠𝑖,𝑑 = 10 − 35.1𝜈𝑏 ,
(1.50)
Figure 1.10 – Elastic modulus as a function of brine volume for first-year sea ice according to [42]
The equations for determining the dynamic elastic modulus according to the Abele and
Frankenstein (1966) and Tabata (1958), respectively, are presented in [8] as follows
𝐸𝑠𝑖,𝑑 = 9.01 − 46.4𝜈𝑏 ,
(1.51)
𝐸𝑠𝑖,𝑑 = 9.27 − 28.2𝜈,
(1.52)
In equations (1.46)-(1.52) 𝜈𝑏 is the brine volume fraction and 𝜈 is total porosity.
Poisson's ratio is one of the main mechanical properties that characterize the behavior of a material
under load. Poisson's ratio 𝜇 is defined as the ratio of the lateral strain to the longitudinal strain in a
homogeneous material for a uniaxial loading condition. The Poisson’s ratio, just like the elastic modulus,
can be dynamic and static (true and effective).
One of the earliest proceedings in which an equation was proposed for determining the Poisson's
ratio is a monograph by Weeks and Assur (1967) [81]. Authors obtained the following formula of
dynamic Poisson's ratio 𝜇𝑖,𝑑 based on in-situ seismic observations of Lin’kov (1958) at Kap Shmit,
Siberia:
𝜇𝑖,𝑑 = 0.333 + 0.06105 exp (
22
𝑇
)
5.48
(1.53)
Figure 1.11 – (𝜇 − 1/3) vs ice temperature, based on in-situ seismic determinations [81]
Murat and Lainey (1982) [49] measured the longitudinal and transverse strains on simply
supported beams loaded in flexure. The tests were performed at different temperatures and loading rates
on columnar ice. As a result, authors proposed the following equations for the effective (static) Poisson's
ratio 𝜇𝑖,𝑠 :
𝜎̇ −0.29
𝜇𝑖,𝑠 = 0.24 ( )
+ 𝜇𝑖,𝑑 ,
𝜎̇1
(1.54)
𝜀̇ −0.29
𝜇𝑖,𝑠 = 0.0024 ( )
+ 𝜇𝑖,𝑑 ,
𝜀̇1
(1.55)
where 𝜎̇ and 𝜎̇1 – stress rate and unit stress rate respectively, kPa∙s-1; 𝜎̇1 = 1 kPa ∙ s−1 ;
𝜀̇ and 𝜀̇1 – strain rate and unit strain rate respectively, s-1; 𝜀̇1 = 1 s−1 .
Orthotropic elasticity. If to take into account the ice growth process of ice feature, it can be said
that the structure of the sheet has three mutually perpendicular mirror planes (two in the vertical plane)
and so possesses orthotropic symmetry [66]. In this case elastic properties of ice formations can be
described by three Young’s moduli 𝐸1 , 𝐸2 , 𝐸3 , by the three shear moduli 𝐺12 , 𝐺13 , 𝐺23 , and to six
Poisson’s ratios 𝜇12 , 𝜇21 , 𝜇13 , 𝜇31 , 𝜇23 , 𝜇32 (where 𝜇𝑖𝑗 = −𝜀𝑗 /𝜀𝑖 and 𝑖 – direction of uniaxial loading).
Schulson E. M., Duval P. (2009) calculated the values of these parameters for 3 texture types of
ice sheets, denoted S1, S2 and S3 ice. These values are presented in Table1.6 and related with rectangular
coordinate system where 𝑋1 and 𝑋2 lie within the horizontal plane of the sheet and 𝑋3 is parallel to the
vertical direction. In ice S1 c-axes are parallel to 𝑋3 , in ice S2 c-axes are randomly oriented in 𝑋1 –𝑋2
plane and in ice S3 c-axes are parallel to 𝑋1 .
The most convenient approach for modeling ice properties is to consider both isotropic and
orthotropic elasticity regarding to ice structure. As it was mentioned in subsection 1.2 ice in general
there are two different types of ice structure in ice floe. Top part of ice field is granular ice with more or
less randomly oriented c-axes. For this type of ice isotropic elastic properties can be considered. On the
23
other hand, bottom layers of ice floe are composed of columnar ice. This type of ice has different
parameters of elasticity and strength in vertical and horizontal directions, so it should be considered as
orthotropic ice.
Table 1.6 – Elastic properties of homogeneous orthotropic sheets of ice Ih at -16 °C [66]
Texture
Young’s modulus 𝐸,
GPa
Shear modulus 𝐺,
GPa
Poisson’s ratio
𝐸1
𝐸2
𝐸3
𝐺12
𝐺13
𝐺23
𝜇12
𝜇21
𝜇13
𝜇31
𝜇23
𝜇32
S1
9.71
9.71
11.8
3.42
3.01
3.01
0.415 0.415 0.224 0.274 0.224 0.274
S2
9.58
9.58
9.71
3.61
3.21
3.21
0.327 0.327 0.344 0.320 0.344 0.320
S3
11.8
9.71
9.71
3.01
3.01
3.42
0.274 0.224 0.274 0.224 0.415 0.415
1.3.2 Compressive strength
The compressive strength is another fundamental property of sea ice. The destruction of ice during
its impact on vertical structures occurs mainly due to the stresses reaching the compressive strength.
Thus, the compressive strength of ice plays the largest role in the magnitude of the load on the vertical
offshore structures.
Uniaxial compressive strength. The typical apparatus of compressive strength measurements is
shown in Figure 1.12 [74]. To conduct the compressive strength test, cores of sea ice usually are cut into
cylinders. The samples are placed in a test apparatus that can generate high loads. A load cell is included
in-line to measure the applied load. Among the important factors affecting the results of the experiment
are the stiffness of the testing machine, the parallelism of the faces of the samples, as well as the aspect
ratio and shape of the samples [66].
Figure 1.12 – Typical apparatus for a uniaxial compression test on ice [74]
24
Many investigators measured the compressive strength of sea ice. Based on their research, it can
be said that the compressive strength of sea ice depends on many factors. Among such factors, one can
especially distinguish the strain rate, the texture of ice, temperature and porosity.
If we talk about the strain rate, then it plays a very large role in the behavior of ice during
compression. It is well known that ice in compression can exhibit both plastic and brittle behavior
depending on the strain rate as presented in Figure 1.13. This property of ice has been observed and
confirmed by many authors.
Figure 1.13 – Schematic compressive stress–strain curves [66]
At low strain rates, usually less than 10-4 s-1, the ice exhibits plastic behavior, which is
characterized by a smooth stress-strain curve and hardening. Also, at these speeds, the compressive
strength of ice is little dependent on grain size, but slowly increases with decreasing grain size [13] as
presented in Figure 1.14 (a).
Figure 1.14 – Peak stress vs grain-size: a) strain-rates between 10-7 and 10-3 s-l; b) strain-rates between
3 ∙ 10-3 and 10-1 s-l
25
The brittle behavior of ice can be observed at strain rates of 10-3 and higher. Stress-strain curves
in this case are close to linear in shape as shown in the Figure 1.13. Also, these curves are characterized
by small "drops" associated with the appearance and development of microcracks. The effect of grain
size in this case is stronger than with plastic behavior as presented in the Figure 1.14 (b).
A consequence of both strain-rate hardening (ductile regime) and strain-rate softening (brittle
regime) is that the compressive strength reaches a maximum at the ductile-to-brittle transition and at
strain rate from 10-4 to 10-3 s-1 [13, 28, 65, 66, 67]. Thus, the transition from plastic to brittle behavior
does not occur instantly. The reason is that the overall inelastic deformation comprises a mixture of
creep via dislocation slip and cracking, the proportion of which varies in the vicinity of the transition
[66]. In Figure 1.13, the transition point is characterized by a strain rate equal to 𝜀̇𝑡𝑐 . The magnitude of
these strain rates, i.e., the rates of ductile-brittle transition, has been studied many times and ranges from
about 10-4 to 10-3 s-1. Since the ice strength is maximum at these rates, today the ice load acting on the
structures is calculated precisely at the strength corresponding to this transition zone. This behavior is
exhibited by both single crystals and polycrystals. Typical results of ice compression experiments
showing the dependence of strength on the strain rate are shown in Figure 1.15, Figure 1.16 and Figure
1.17.
Figure 1.15 – Peak stress vs strain-rate in experiments of Cole D. M. [13]
Figure 1.16 – Strain rate dependence of compressive strength of ice according to Gold L. W. [28]
26
Figure 1.17 – Unconfined compressive strength of single crystals of ice vs. strain rate [66]
It is also worth considering the effect of load direction on ice strength. The direction of the load
on the strength was investigated to a much lesser extent. Nevertheless, today the fact is known that the
compressive strength of ice reaches its minimum when the load direction is oriented at an angle of 45
degrees to the c-axes [66]. This feature is confirmed by the studies of Payton (1966) [67] and
Carter (1971) [66] as shown in figures below. In Figure 1.18 the loading-angle notation is as follows:
the first number gives the angle between the axis of the test cylinder and the vertical, while the second
number gives the angle between the sample and the c-axis of the single ice crystal being tested.
Figure 1.18 – Average failure strength in compression (solid circles) and in direct tension (open
circles) versus sample orientation [67]
Figure 1.19 – Unconfined brittle compressive strength (and tensile strength) of single crystals of ice at
−10 °C vs. orientation of the basal plane [66]
27
Schulson E. M. and Duval P (2009) collected a data from uniaxial compression experiments
performed by many authors such as: Jones (1997), Schulson et al. (2005) and Shazly et al. (2006) as
presented in the Figure 1.20. Taken collectively, the measurements indicate that dynamic strength 𝜎𝑐 ,
MPa, at -10 °C may be described by the following relationship with variance r2 equal to 0.71 [66]:
𝜎𝑐 = 9.8(𝜀̇)0.14 ,
(1.56)
Figure 1.20 – The compressive strength of ice at -10 °C according to Schulson and Duval (2009)
Timco and Frederking (1990) [75] developed some equations for calculating the compressive
strength of sea ice sheets. They compiled 283 compressive strength test results on first-year sea ice. They
derived the equations for the uniaxial compressive strength of first-year sea ice for several different grain
structures. For horizontally loaded columnar ice, the uniaxial compressive strength 𝜎𝑐 , MPa, is defined
as follows
𝜈𝑡
),
270
𝜎𝑐 = 37(𝜀̇)0,22 (1 − √
(1.57)
where 𝜀̇ – strain rate, s-1;
𝜈𝑡 – total porosity, ‰.
For vertically loaded columnar ice, the uniaxial compressive strength 𝜎𝑐 , MPa is defined as follows
𝜈𝑡
)
200
𝜎𝑐 = 160(𝜀̇)0.22 (1 − √
(1.58)
For granular ice, the uniaxial compressive strength 𝜎𝑐 , MPa is defined as follows
𝜎𝑐 = 49(𝜀̇)0.22 (1 − √
𝜈𝑡
)
280
(1.59)
The applicable range of strain rates for these equations is from 10−7 s −1 to about 2 ∙ 10−4 s −1.
Kovacs (1997) [39] obtained the following equation for the unconfined compressive strength of
first-year sea ice floes
𝜎𝑐 = 𝐵2 (𝜀̇)1/𝑛 𝜈𝑡𝑚 ,
28
(1.60)
where 𝐵2 , 𝑛 and 𝑚 – empirical parameters, values of which can be taken equal 𝐵2 = 2700, 𝑛 =
3, 𝑚 = −1.
Compressive strength of confined ice. The above features and the presented equations are obtained
on the basis of uniaxial compression strength tests of ice. However, when ice fields interact with
structures, a complex stress-strain state arises in them. In this case, it would be incorrect to characterize
this state by the ice strength characteristics obtained from uniaxial tests. Therefore, in recent years,
studies on the strength of ice under “confined” conditions corresponding to biaxial and triaxial stress
states of the ice sheet have begun to appear more often. The confinement induces biaxial and triaxial
stress states, which have a large effect on the strength of the ice and on its mode of failure [66].
There are several ways to perform biaxial and triaxial compression experiments on ice samples,
such as: triaxial loading using pressurized fluid, passive proportional loading, true triaxial loading and
indentation. The most convenient is the true triaxial loading method, because the three principal stresses
𝜎1 , 𝜎2 and 𝜎3 (𝜎1 is taken as the most compressive stress) can be varied independently and all loading
paths in principal stress space can be explored, including biaxial ones. The loading paths in this method
of testing are defined by the ratio (𝜎1 : 𝜎2 : 𝜎3 ) = (1: 𝑅21 : 𝑅31 ), where 𝑅21 = 𝜎2 /𝜎1 ≤ 1 and
𝑅31 = 𝜎3 /𝜎1 ≤ 1.
To consider the behavior of confined ice during compression, it is better to consider granular and
columnar ice separately. As mentioned earlier, polycrystalline granular ice has a random orientation of
single crystals and therefore can be considered as isotropic.
Some features of the behavior of granular ice during compression can be made based on the results
of experiments performed by Weiss and Schulson (1994) [82]. Authors used cubes (l55 mm on edges)
of fresh-water ice to investigate the failure of granular ice under multiaxial compressive loading using
solids platens at temperatures -10 °C, -20°C and -40 °C at a strain-rate of 10-3 s-l. Some particular triaxial
tests with 𝜎2 = 𝜎1 (and 𝜎3 = 𝑅𝜎1 ), as well as biaxial tests (𝜎3 = 0, 𝜎2 = 𝑅𝜎1 ) showed that the
intermediate principal stress 𝜎2 bas no significant effect on the compressive strength in either the lowconfinement or the high-confinement regime. Indeed within the scatter of the data, failure stress 𝜎1𝑓
biaxial tests is roughly equal to the uniaxial failure stress [82].
With triaxial compression of granular ice, an increase in compressive strength is observed. Under
triaxial loading the confining stress acts upon all inclined primary cracks, regardless of orientation.
Consequently, it raises the maximum principal stress to effect sliding and crack growth, thereby
increasing strength [66]. As in case of biaxial tests the triaxial strength appears to be independent of the
intermediate principal stress. In this regard, many authors determine the strength of ice by the stress
difference (𝜎1 − 𝜎3 ) known as the differential stress at failure. Generally speaking, under lower
confinement, the stress difference increases with increasing pressure and under higher confinement, the
stress difference is independent of pressure as presented in Figure 1.21 and Figure 1.22. This feature was
29
investigated by many authors [34, 59, 66, 82]. In the case of columnar ice, a different picture is observed
due to anisotropic properties.
Figure 1.21 – Maximum differential stress vs. confining pressure for fresh-water granular ice of 1 mm
grain size [66]
Figure 1.22 – Failure stress of ice vs. confining ratio according to Weiss and Schulson (1994) [82]
1.3.3 Tensile and shear strengths
The tensile and shear strengths of ice are less studied than compressive one, because of difficulties
of conducting experiments. However, some authors presented results of their investigations described
below.
The tensile strength of ice increases with decreasing temperature, but the values of ice tensile
strength can be several times less than the compressive strength. There are several methods to investigate
30
the tensile strength of ice, but the most usable is direct tensile test. If talk about single crystals of ice,
sufficient study of tensile behavior was performed by Carter (1971). His study was described by
Schulson and Duval (2009) [66]. Carter performed over 600 tests on cylindrical-shaped Ø50x150 mm
crystals of fresh-water ice loaded under either tension or compression at strain rates from 4 ∙ 10-7 s-1 to
2.5 ∙ 10-1 s-1 at temperatures from -30 °C to 0 °C. He also varied the orientation of the crystallographic
c-axis, from angle equal from 0° to 90° with respect to the tensile axis. He examined the following facts
which was described earlier:
- when deformed at strain rates below about 10-4 s-1 ice single crystals showed ductile behavior
and while at rates above about 10-3 s-1 they showed brittle behavior. Between these limits they exhibit
transitional behavior.
- the brittle tensile strength is essentially independent of strain rate and temperature but depends
strongly upon the orientation of the c-axis as presented in Figure 1.23.
In the case of ice polycrystals, the situation is more complicated due to the influence of the
structure and relative position of the crystals in the ice texture. Many factors influence tensile strength
such as: temperature, grain size, growth texture and others. It should be mentioned that under
“quasistatic” loading polycrystals fail under a tensile stress that is essentially independent of strain rate.
Figure 1.23 – Tensile strength of single crystals of ice at -10 °C vs. orientation 𝛼 of the basal plane,
where 𝛼 denotes the angle between the c-axis and the tensile load [66]
The most famous equations for determining the tensile strength of ice were obtained by Dykins
(1967) [21] and Weiss and Meyssonnier (2001) [83]. These equations are used in international standards
such as ISO/FDIS 19906:2019(E) [33].
For the vertical direction of columnar ice tensile strength 𝜎𝑡 , MPa is defined as follows
𝜈𝑏
)
310
𝜎𝑡 = 2.2 (1 − √
31
(1.61)
For the horizontal direction of columnar ice tensile strength 𝜎𝑡 , MPa is defined as follows
𝜈𝑏
)
140
𝜎𝑡 = 1 ∙ (1 − √
(1.62)
Significant measurements of tensile strength of columnar ice were performed by Richter-Menge
and Jones (1993) [58]. Tensile load was applied along the c-axes of the test specimens, which were
perpendicular to the growth direction of the ice. The data presented in this paper provide information on
the tensile behavior of columnar, first-year sea ice at four different temperatures (-20°, -10°, -5° and
-3°C) and two strain rates.
Shear strength of ice was measured by Timco and Frederking (1986) [76]:
𝜈𝑏
)
390
𝜎𝜏 = 1.5 (1 − √
(1.63)
1.4 Ice-structure interaction mechanism in case of vertical structures
The interaction of ice with vertical structures can occur under different scenarios and limiting
processes. There are 3 limiting processes that are associated with factors limiting the ice load on
structures [33, 54]: limit stress process, limit energy (or momentum) process, limit force process
which are presented in Figure 1.24.
In the case of limit stress process, the ice formation has enough driving force to begin its
destruction during the interaction and its cutting through the support to its full width.
In the case of limit energy (or momentum) process, the ice load is limited by the kinetic energy of
the ice formation, depending on its speed and mass. When an ice field collides with a structure, it stops
with partial cutting and dissipation of energy.
In the case of limit force process, the load is limited by natural forces acting on the ice formation.
After the ice formation stops, it can continue to act on the structure, transferring the load to it under the
action of external forces (wind, current, wave or other ice formations).
When assessing the impact of ice on the designed structure, several of the above scenarios are
considered, as well as their combinations. Moreover, as a rule, the smallest of the values obtained is
taken as the resulting load. However, due to the fact that the choice of ice regime parameter values to
the present is a difficult task, the calculated load values should be carefully justified.
Ice load during ice failure has the most important role in assessing the critical condition of a
structure. The mode of ice failure against the structure has a significant effect on the magnitude of the
ice action. The failure mode for sea ice (e.g. crushing, shear, flexure, creep) depends on parameters such
as ice thickness, presence of ridges, ice velocity, ice temperature, structure profile and plan shape, etc.
Also, different modes can occur on the same structure type depending on ice conditions and
32
interaction velocity, even during the same event as presented in Figure 1.25.
Figure 1.24 – Limit factors of ice loads [54]
Figure 1.25 – Failure modes of ice during interaction: a) crushing failure; b) bending failure;
c) splitting failure; d) buckling failure; e) pushing floes [42]
The process of ice impact in case of crushing failure mode is shown schematically in the Figure
1.26. One of the main features of this type of interaction is presence of “high pressure zones”. These are
small areas that quickly appear and disappear randomly at the ice-to-structure contact, on which the main
interaction force is concentrated as presented in Figure 1.26. At any one time, there are several HPZs,
and they occupy only a small fraction of the nominal contact area, generally but not always towards the
middle and away from the edges [54]. The presence of high-pressure zones is justified by the
heterogeneous structure and properties of ice and is confirmed by many authors. As a rule, currently
used ice load calculation methods, including in normative documents, do not take into account the
presence of HPZs. Instead, the so-called nominal contact area is use, which is the projection of the
structure onto the ice formation.
Figure 1.26 – Schematic picture of ice action during ice-structure interaction in case of ice sheet
compression (1 – ice sheet; 2 – structure; 3 – spalls and extrusion; 4 – high pressure zones; 5 –
pressure distribution over the contact surface) [35]
33
If circumferential crack is observed the bending type of flexural failure is occurred. Also, many
radial cracks can be observed behind this crack. Splitting failures are usually observed when the
interacting ice sheet has low lateral confinement or small size. The buckling type of flexural failure is
governed by a buildup of curvature in the ice sheet or occurs under special wave conditions. In this
regard, this process can be observed very rarely. The last type of interaction corresponds to the limit
force process explained before. The floes split and pass around the structure, while wind and wave
actions contribute significantly to the total force [53]. Conditions that induce ice failure by flexure
generally result in smaller ice actions than for crushing [33]. In addition, the interaction of ice with the
structure can cause dynamic loads, which is associated with the periodic destruction of ice during
compression.
1.5 Researches in the field of numerical modeling of ice problems
The problem of assessing the impacts of ice fields regarding the destruction of ice is a highly
specialized task. The issue of numerical modeling of the destruction of ice is little studied. Some studies
with different objectives in the field of numerical modeling of ice have been performed by several
researchers: Derradji-Aouat (2005) [18], Bjørnar Sand (2008) [63], J. Pernas-Sánchez et al (2012) [55],
Rüdiger von Bock und Polach and Sören Ehlers (2013) [60], Aleksandrov A.A. et al (2014, 2017) [2, 3],
Li Liang (2014) [44], Politko V.A. (2018) [56] and others.
J. Pernas-Sánchez et al (2012) performed numerical simulations of experiments conducted by
Pereira et al. (2006) and described by Carney et al. (2006) [12]. In these experiments, cylindrical ice
projectiles were launched onto a steel plate bounded to a strain gauge using a helium gas gun. The
experiment equipment is presented in Figure 1.27. The size of the ice cylinders was 17.5 mm in diameter
and 42.2 mm in length. The tests were carried out at different impact speeds.
Figure 1.27 – Equipment in the experiments described by Carney et al. (2006)
To perform the numerical simulation, three different numerical methods were used: the finite
element method in the Lagrange formulation (Lagrange), the arbitrary Lagrange-Eulerian method (ALE)
34
and smoothed particle hydrodynamics (SPH) method. Figure 1.28 shows three models of an ice cylinder
made by J. Pernas-Sánchez et al.
Figure 1.28 – Numerical models of the ice cylinder in the study of J. Pernas-Sánchez [55] (left –
Lagrange, in the center – ALE, right – SPH)
To model the ice material at high strain rates, the Drucker-Prager model was chosen. The process
of destruction of the ice cylinder when striking a steel plate at a speed of 152 m/s, as well as the
dependence of the impact force on the time of the experiment are shown in Figure 1.29.
The Lagrangian finite element solver showed the greatest convergence with experimental data,
despite the strong distortion of elements. Based on the results of this study, it can conclude that the
Drucker-Prager model shows high convergence with experimental data at high strain rates. However,
the real drift velocities of ice fields have values ten times less than in this experiment. In this connection,
additional computational studies of the applicability of the Drucker-Prager model are required for
modeling the impact of the ice formations on the offshore facilities.
Figure 1.29 – The process of destruction of the ice cylinder on impact (left) and the dependence of the
impact force on time (right) for the three numerical methods (from left to right: Lagrange, ALE, SPH)
35
Aleksandrov A.A. et al (2014, 2017) in their investigations used various approaches to modeling
the ice material.
In [2], the authors carried out a computational study in which ice parameters were determined
using the ANSYS software (the inverse problem was solved). Prior to the beginning of numerical
simulation, the authors conducted experiments on three-point bending of ice beams (700x100x100 mm)
and uniaxial compression of ice cubes (100x100x100 mm). Samples were obtained from the bottom
layer of ice, which formed naturally on the surface of the freshwater pond. Loading was carried out using
a hydraulic cylinder. The indenter had a cylindrical shape when testing beams and a plate shape when
testing cubes. A grid of dots, shown in Figure 1.30, was applied to the central part of the samples before
each test using a water-based paint. The load force, the geometric dimensions of the sample, and the
deformations at the grid nodes on the ice surface were measured.
Figure 1.30 – Grid of dots in the center of the ice sample [2]
The mesh motion during the tests was recorded using a three-dimensional computer vision system
and digital image correlation for contactless measurements of shape and deformation. The calculation
of the elastic modulus of ice was made using the formulas of the standard beam theory. An example of
the dependence of the strain values of the samples on the load is shown in Figure 1.31.
To determine the mechanical properties of ice, the authors used an optimization method using the
strain values obtained from the experiment. Two models of the material were used to model the ice: the
model of an ideally elastic body, which is deformed according to Hooke's law, and the viscoelastic
Maxwell model. Shear modulus 𝐺, MPa, in these material models was in the following dependence on
the modulus of elasticity 𝐸, MPa:
𝐺=
𝐸
,
2(1 + 3𝜈)
where 𝜈 – Poisson's ratio, which was taken equal to 0.33.
36
(1.64)
Figure 1.31 – Dependence of strain on the load at the center of the samples [2]
When determining the parameters of ice by the optimization method for the case of the Maxwell
model, the viscosity coefficient 𝜂 was also determined. The results of optimization are presented in
Tables 1.7-1.8.
Table 1.7 – Values of the properties of ice obtained using the standard beam theory and optimization by
the finite element method [2]
Sample
Standard beam theory
Numerical simulation
(elastic material
model)
№1
𝐸 = 390 MPa
𝐸 = 490 MPa
𝐸 = 500 MPa; 𝜂 = 250 GPa ∙ s
№2
𝐸 = 200 MPa
𝐸 = 220 MPa
𝐸 = 230 MPa; 𝜂 = 100 GPa ∙ s
№3
𝐸 = 310 MPa
𝐸 = 390 MPa
𝐸 = 400 MPa; 𝜂 = 420 GPa ∙ s
Numerical simulation
(viscoelastic material model)
Table 1.8 – Differences in displacements calculated by various methods relative to the experiments, mm
Sample
№1
№2
№3
The number of
compared points
Standard
beam theory
Numerical simulation
(elastic material
model)
Numerical simulation
(viscoelastic material
model)
1
0.031
0.023
0.021
25
-
0.029
0.026
1
0.091
0.119
0.088
25
-
0.121
0.107
1
0.029
0.022
0.018
25
-
0.029
0.025
37
Thus, there is little difference between the results of numerical modeling using elastic and
viscoelastic material models, but these results are quite different from the results of the calculating the
parameters according to the standard beam theory. According to the authors, the Maxwell viscoelastic
model shows the greatest convergence with experimental data. In this regard, adding the additional
parameters into the material model (in this case, the viscosity coefficient) allows to describe the behavior
of ice during deformation more accurately and to obtain more acceptable values of its characteristics.
Also, during bending tests at three points, various types of specimen failure were recorded. The
authors identified three main types of damage: fracture due to the formation of a straight crack, fracture
in the form of a zigzag crack, and fracture in the form of a curve crack, presented in Figure 1.32.
When searching for the relationship between the type of crack and the parameters of the external
load (loading rate, limit values of the load) or the physical and mechanical properties of ice, the authors
did not find a correlation.
Figure 1.32 – Typical types of cracks: a) straight, b) zigzag, c) curve [2]
In [3] Aleksandrov et al. considered the applicability of numerical modeling for solving the
problem of scaling ice loads in the study of ice impacts on vertical structures in the ice pool. According
to the authors, numerical modeling can help solve the problem of determining the physical and
mechanical properties of ice formations.
First, the authors conducted a series of tests in the ice pool. The block-conductor model was towed
in an ice field with low speeds of 0.04÷0.15 m/s, the thickness of the ice field was 20 and 30 mm. During
the tests, components of the ice load were recorded.
At the next research step, the authors performed a numerical simulation of the experiments. They
proposed the use of the isomorphism principle of the limiting surface to obtain the necessary
combination of material parameters. The principle of isomorphism of the limiting surface states that the
characteristic form of the limiting surface of the material is preserved. In this case, the limit surface
should be understood as the surface in the axes of the main stresses, upon reaching which the material
stops to obey the Hook’s law of elastic deformation (yield surface). The achievement of yield surface
by stresses means the beginning of plastic flow or -is a criterion of brittle fracture.
Based on the proposed principle, the authors solved the numerical problem of introducing a blockconductor model into an ice field in an ice pool. For a finite element modeling, the ANSYS software
38
was chosen. The block conductor used in the experiments, as well as the finite element model are
presented in Figure 1.33.
Figure 1.33 – Block-conductor during tests in the ice pool and FE model [3]
The surface of the block-conductor was modeled by membrane non-deformable finite elements.
The ice field was modeled using 4-node finite elements and fixed at the boundaries in all degrees of
freedom. The movement of the indenter is set using the pilot point and does not exceed 0.1 m.
The authors chose the Drucker-Prager Cap model as a model for the ice material. The shape of the
yield surface is approximated according to the data obtained as a result of tests for biaxial compression
of cubic samples of natural ice [65]. The data are presented in Figure 1.34. The approximation is
performed by solving an optimization problem for the control parameters of the model.
To restore the yield surface in the field of tensile stresses, the authors accepted the fact that the
compressive strength exceeds the tensile strength 5 times. When modeling the modulus of elasticity is
taken from the tests performed on the bending of cantilever beams 𝐸 = 23.5 𝑀𝑃𝑎, Poisson's ratio 𝜈 =
0.33.
Figure 1.34 – Chart of destruction of ice cubic samples under biaxial compression in the axes of main
stresses [65]
The results of the numerical simulation of experiments were compared with experimental data.
Comparison of data is presented in Table 1.9. The difference between the results of numerical simulation
and experimental research does not exceed 10%, which may indicate a close correspondence of the
Drucker-Prager Cap model to the actual behavior of ice under loads. The authors argue that the
39
coordination of the results of numerical modeling perfectly illustrates the application of the isomorphism
principle of the limiting surface to the solution of the problem of introducing a body into an ice field.
But looking at Figure 1.33, I can conclude. that the authors did not complete the qualitative creation
of the mesh. The ice field in thickness has only 2 finite elements and thus cannot adequately describe
the process of interaction with the block conductor.
Table 1.9 – Comparison of the results of experiments and numerical simulation [3]
Force of drag, N
Ice thickness, mm
Experiment Numerical model
20
36
39
30
64
59
Politko (2018) analyzed the applicability of the Mohr-Coulomb model for modeling the impact of
an ice field on vertical structures. He used the element erosion technique to model the ice fracture. Two
experiments were considered to verify the applicability of material model: interaction of a rectangular
horizontal stamp with the ice field [70] and laboratory model testing of the interaction of ice field with
a four-support structure in the basin of the Krylov Research Center [36].
When modeling the experiments [70], the physical and mechanical properties of ice presented in
Table 1.10 were determined. A numerical model of stamp-ice interaction process and boundary
conditions are shown in Figure 1.35. A graph comparing the ice load between an experiment and
numerical simulation is presented in Figure 1.36. Analyzing Figure 1.36, it can be said that the numerical
model showed similar results regarding the nature of the ice load in the contact zone of ice and a
rectangular stamp.
Table 1.10 – Physical and mechanical properties of ice in the numerical model created by Politko (2018)
Property
Units
Value
Elastic modulus 𝐸
MPa
1000
Poisson's ratio 𝜈
-
0.3
Angle of internal friction (to assign parameters of the Mohr-Coulomb
model)
degrees
35
Cohesion (to assign parameters of the Mohr-Coulomb model)
MPa
0.6
Ice salinity
‰
6.2
Ice temperature
°С
-2.7
40
Figure 1.35 – Numerical model for introducing a rectangular stamp into an ice field [56]
Figure 1.36 – Graph of ice load acting on stamp: 1 – field tests; 2 – numerical simulation [56]
When simulating laboratory tests of the interaction of an ice field with a four-support structure in
the ice pool of the Krylov Research Center, the same characteristics of ice are used, presented in
Table 1.10. The model that was used in the laboratory tests is shown in Figure 1.37. The numerical
model presented by the author is shown in Figure 1.38.
Figure 1.37 – Model of a four-support hydraulic structure used in field tests [36]
Figure 1.38 – Numerical model of the impact of the ice field on a four-support structure according to [56]
41
Numerical experiments, as well as laboratory experiments, were carried out at three different
angles of ice effect on the model: 𝛼 = 0°, 𝛼 = 30°, 𝛼 = 45°. The main purpose of the model tests in the
Krylov Research Center basin was to study the mutual influence of the supports at various angles of
action of the ice field. As a result of the studies, the influence coefficient of the supports was determined,
which shows how many times the load on the four-support structure is more than on one-support. The
picture of the destruction of ice during interaction in numerical simulation and during model tests, in
general, coincides and is presented in Figure 1.39. The results of the experiment and numerical
simulation are given in the Table 1.11.
Figure 1.39 – The destruction of ice in the interaction with the four-support structure at angle 𝛼 = 30°:
a) in numerical simulation; b) during the model tests [56]
Table 1.11 – The coefficient of influence of the supports of the four-support structures [56]
The angle of impact of
the ice field 𝛼
The results of the experiment, at a speed 𝑉, m/s
0.01
0.05
Numerical simulation
results
0°
1.73
2.34
2.1
30°
2.9
3.62
3.4
45°
2.53
2.51
2.5
Although the results of numerical simulation of both experiments show a good correlation with
field tests, the use of the Mohr-Coulomb model may be a wrong decision in case of assessing non-brittle
behavior, since this model does not take into account the average principle stress 𝜎2 when estimating the
complex stress-strain state of the ice field. Additional studies are needed to confirm the possibility of
application of the Mohr-Coulomb model for determination of the onset of plastic behavior. There is also
no discussion about the values of accepted parameters of the Mohr-Coulomb model.
As can be seen in Figure 1.35, the finite element mesh is incompatible in the region of transition
to larger finite elements, which indicates the use of contact pairs for modeling the interaction of these
two parts of the ice field. In my opinion, this approach is incorrect due to the relatively large errors that
occur during the deformation process. For adequate transmission of stresses throughout the body of the
simulated field, it is necessary to create a mesh with compatible nodes.
42
In addition, one of the main assumptions in these calculations is the fact that only brittle fracture
of ice without any hardening is considered at relatively fast strain rates. In this case, it is assumed that
for ice before fracture, only elastic deformations appear. Plastic properties of ice and strain-rate
dependence are not taken into account.
Further, in his work, the author investigated the influence of the shape of the structure and the
presence of zones of frozen ice on the ice load. Table 1.12 presents the results of numerical simulation,
which shows the values of the maximum ice load on the supports of various shapes in the plan and the
corresponding shape factors 𝑚.
Table 1.12 – The results of the study of the influence of the structure shape and the presence of frozen
ice on the ice load
Ice load, MN (support shape factor 𝑚)
Ice impact case
Lack of ice freezing with
support
The presence of ice freezing
with support
Front edge in the form of
rectangle
semi-circular
polyhedron
5.2 (1)
4.7 (0.9)
4.8 (0.92)
6.4 (1.23)
6.7 (1.29)
6.6 (1.27)
The author compared these results with the support form factor specified in Russian code
SP 38.13330.2018 [68] and got two main conclusions:
- in the case of the absence of freezing of ice for the front face of a semicircular outline and in the
form of a polygon, the simulation results show the support shape factors higher than in [68], namely
0.90 and 0.92, instead of 0.83;
- in case of freezing of ice with a rectangular support, the support shape factor becomes equal to
1.23 instead of 1.
Thus, the values of support shape factors obtained by the author are consistent with
SP 38.13330.2018 only for the case of a rectangular support in the absence of ice freezing and a semicircular outline (polyhedron) in the presence of freezing. The value 0.9 of shape factor for the case of
semicircular outline form in the absence of a freezing effect, obtained by the author, is consistent with
another code – ND 2-020201-015 [52].
When analyzing studies in the field of numerical modeling of ice impacts on vertical structures, it
is worth highlighting the following numerical experiments: Gürtner et al. (2009) [30], Gürtner et al.
(2010) [29], Hilding et al. (2011) [31], Hilding et al. (2012) [32]. These papers are united by the fact
that the authors performed a numerical simulation of the actions of level ice fields on the Norströmsgrund
43
lighthouse. This lighthouse was chosen by the authors, because it was involved in the Low Level Ice
Forces (LOLEIF) and Measurements on Structures in Ice (STRICE) projects in 1999-2003 years. In the
LOLEIF project, this lighthouse located in the subarctic Gulf of Bothnia was equipped with load panels
covered 162º of the lighthouse perimeter and special devices that monitor environmental conditions.
Thus, comprehensive measurements of ice loads were conducted, and full-scale ice load data was
obtained. The investigated interaction velocity in these papers is 0.15 m/s, which represents the relatively
low interaction velocities for real ice-structure interaction situations.
In these works, the authors used the principles of cohesive zone model (CZM) to simulate the
destruction of ice during interaction with the structure. Finite element models are presented in Figure
1.40, Figure 1.41 and Figure 1.42. The interpretation of this method in a finite element analysis is called
the cohesive element model (CEM) and will be described later in this thesis. The main idea of the method
is to model an ice body in the form of a set of volumetric finite elements having ice properties and
connected with cohesive elements of zero thickness, which are possible crack paths.
Figure 1.40 – Finite element model of ice interaction with Norströmsgrund lighthouse [30]
Figure 1.41 – Finite element analysis of ice interaction with Norströmsgrund lighthouse: a) Gürtner et
al. (2010) [29]; b) Hilding et al. (2011) [31]
44
Figure 1.42 – Simulation of ice-structure interaction process using CEM [32]
In these researches authors defined cohesive element traction-separation law also called
Tvergaard-Hutchinson law (Figure 1.43) for cohesive elements and elastic-plastic material with
softening for ice bulk elements to achieve the macroscopic ice crushing behavior. The Mises material
model was used for bulk elements. Cohesive element law is characterized through its piecewise
continuous curve with softening. Cohesive properties are different in vertical and horizontal direction to
accommodate anisotropy of fracture of sea ice. In all these studies, the von Mises yield criterion was
used to describe the plasticity of the ice material in bulk elements.
Figure 1.43 – Cohesive element traction-separation law defined in [29, 30, 31, 32]
For monotonic loading in fracture mode-I, the fracture energy, i.e. the area under the cohesion law
curve can be written as [30]:
𝐺𝐼𝐶 = 0,5𝑇1𝑚𝑎𝑥 𝑢1𝑐 (1 − 𝜆1 + 𝜆2 ),
(1.65)
where 𝑇1𝑚𝑎𝑥 – peak traction accommodated by the material in the fracture process zone;
𝜆1 , 𝜆2 – separation values on cohesive traction-separation curve (Figure 1.43);
𝑢1𝑐 – maximum allowable separation.
The issue of choosing values for the parameters of the materials used has been discussed many
times. Difficulties arise because today there are no qualitative test data that determine the values of
fracture parameters. The situation is complicated by the fact that the energy fracture parameters are
highly dependent on the mesh size, so in some cases they had to be scaled tens of times [29, 31]. The
values of material properties used in these papers are shown in Table 1.13.
45
Table 1.13 – Values of material properties used in [29, 30, 31, 32]
Paper title
Bulk element properties
Numerical simulation of ice
action to a lighthouse / Gürtner
et al. (2009)
von Mises yield criterion
𝐸 = 1000 GPa;
= 910 kg/m3 ;
𝜇 = 0.3;
𝑦 = 0.0007 MPa;
Numerical modelling of a full
scale ice event / Gürtner et al.
(2010)
von Mises yield criterion
𝐸 = 5 GPa;
= 910 kg/m3 ;
𝜇 = 0.3;
𝑦 = 1,5 MPa
Simulation of ice action loads
on offshore structures / Hilding
et al. (2011)
Simulation of Loads from
Drifting Ice Sheets on Offshore
Structures / Hilding et al.
(2012)
von Mises yield criterion
𝐸 = 5GPa;
= 910 kg/m3 ;
𝜇 = 0.3;
𝑦 = 2 MPa
von Mises yield criterion
𝐸 = 5GPa;
= 910 kg/m3 ;
𝜇 = 0.3;
𝑦 = 2 MPa
Cohesive element properties
𝐺𝐼𝐶 = 20 J/m2 (vert.);
𝐺𝐼𝐼𝐶 = 28 J/m2 (vert.);
𝐺𝐼𝐶 = 52 J/m2 (horiz.);
𝐺𝐼𝐼𝐶 = 52 J/m2 (horiz.);
𝜆1 = 0.08 (vert.);
𝜆2 = 0.45 (vert.);
𝜆1 = 0.10 (horiz.);
𝜆2 = 0.55 (horiz.)
𝑇𝑚𝑎𝑥 = 1 MPa (vert.);
𝑇𝑚𝑎𝑥 = 1.1 MPa (horiz.);
𝐺𝐶 = 5200 J/m2 ;
𝜆1 = 0.02;
𝜆2 = 0.55 (vert.);
𝜆2 = 0.45 (horiz.)
𝑢𝑐 = 6.8 mm (vert.);
𝑢𝑐 = 7.1 mm (vert.);
𝑇𝑚𝑎𝑥 = 1 MPa (vert.);
𝑇𝑚𝑎𝑥 = 1.1 MPa (horiz.);
𝐺𝐶 = 5200 J/m2
𝑇𝑚𝑎𝑥 = 1 MPa (vert.);
𝑇𝑚𝑎𝑥 = 1.1 MPa (horiz.);
𝐺𝐶 = 5200 J/m2
Figure 1.44 – Simulated horizontal global ice force on the lighthouse foundation and comparison to
measured loads according to Gürtner et al. (2009) [30]
This modeling method (CEM) shows close to real behavior and the destruction of the ice cover
when exposed to a vertical structure. Features such as the appearance of radial cracks, cyclic loading and
ice extrusion can be observed. However, in my opinion it is impossible to speak of a good agreement
46
between the results and field data, especially because results of modeling in these papers were often
compared with only a few field cases (most often with one). Numerical experiments in these papers have
several disadvantages. The authors did not perform a qualitative selection of constitutive model and
parameters of the ice material, which led to a large difference between the model and the real picture of
the interaction. As it usual for almost all simulations, the authors did not consider the dependence of
strength on the strain rate, temperature, as well as the ductile behavior of ice before fracture. The big
assumption is the use of von Mises plasticity for bulk elements. This plasticity model does not take into
account the influence of hydrostatic (mean) stress on total ice strength. This effect was described earlier
in subsection 1.3. So, we can say that there is no qualitative material modeling of ice in these papers.
One of the latest investigations in the field of ice-structure interactions was made by Wang et al.
(2019) [79]. In these paper authors also used the cohesive element method to model the ice actions on
Norströmsgrund lighthouse with new modeling techniques. First of all, some ice bulk elements with
random numbering are deleted to simulate natural ice sheet with initial defects. Then, the elastic modulus
was calculated by the empirical formula and assigned to the ice elements. Thirdly, cohesive elements
are inserted in ice sheet model which is meshed with tetrahedron elements. According to the authors the
use of hexahedral elements will cause “zig-zag” crack paths that are not applicable for the cohesive zone
method. A structured mesh with rectangular elements will create a 45° crack which will travel 2 times
longer than the actual one, leading to an extra energy consumption. The numerical model used in
calculations is presented in Figure 1.45. Ice element is assumed as elastic-perfect plastic material.
Figure 1.45 – Numerical model of ice sheet and lighthouse used in [79]
In general Wang et al. considered 3 different design cases: ice sheet without any initial defects,
ice sheet has only initial defects and ice sheet has initial defects and changed elastic modulus. Horizontal
force on the lighthouse and failure process of ice sheet are analyzed and compared with field data
measured on lighthouse. All of three cases showed different failure pattern as shown in Figure 1.46.
Values of simulated ice forces and comparison with field data is presented in Table 1.14.
47
Figure 1.46 – Ice failure patterns in three cases according to [79]
Table 1.14 – Comparison of the simulated mean horizontal forces with the measured data [79]
Case
Mean horizontal force, MN
Error
Case 1
3.21
15.5%
Case 2
2.59
6.8%
Case 3
2.92
5.7%
Measured results
2.78
-
Again, as in previously reviewed articles and theses, the same assumptions and disadvantages are
present. The Mises model was used as a model of plasticity, which in my opinion cannot be applied as
material model for such complex material as ice. Also, the authors did not take into account the influence
of temperature, strain-rate effects. As can be seen from Figure 1.45, a more or less good-quality ice field
mesh was made only in a small area near the interaction region. Away from this region, there is a rapid
enlargement of the finite element mesh up to 1 element in the ice field thickness, which in my opinion
leads to an incorrect distribution of stresses. Given all of the above, along with the fact that the
comparison is made only for 1 case of interaction, it cannot be said that the error values given in
Table 1.14 correspond to reality.
Nevertheless, the authors used some interesting methods. For example, they used a special
subroutine to account for the buoyance and gravity forces, acting only on spalled elements. Preliminary
removal of elements from the model to generate initial damage is also an interesting solution.
48
Good results in the field of numerical studies of the interaction of ice structures were obtained by
Sand (2008) [63]. In modeling, the author used the finite element method in conjunction with the
technique of removing elements and studied the effect of ice on both inclined and vertical structures.
Sand in his thesis considered the applicability of some failure surfaces for ice modeling. Among the
examined, the author identified the following:
- Hill yield criterion;
- Parisau parabolic criterion;
- Reinicke and Remer’s failure criterion
- Horrigmoe and Zeng’s failure criterion
According to Sand, criterion Horrigmoe and Zeng’s failure criterion is most appropriate, because
this orthotropic criterion contains not only a quadratic terms in stresses, but also linear terms in stresses
and has the dependence of the yield strength on the mean stress. The search for the values of the
parameters included in this criterion was performed by approximating the results of field experiments
data from papers of Schulson. Resulting failure surface used in numerical calculations is presented in
Figure 1.47. Moreover, the author noted that the obtained values can also be used in the Hill criterion.
Figure 1.47 – Fitted three-dimensional plots of Horrigmoe and Zeng’s failure criterion [63]
The behavior of the crushed or cracked ice was taken into account by replacing the “undamaged”
finite elements with the “damaged” ones with changed properties. This is so-called transformation of
state method that is shown in Figure 1.48. As in paper presented higher only elastic behavior of ice
before crushing was considered, and author did not take into account the dependence of strength on the
strain rate and temperature.
In terms of modeling the ice-vertical structure interaction problem, good results were obtained.
Author conducted several simulations of ice actions of Norströmsgrund lighthouse which was mentioned
before. Sand compared global ice forces with measured field data and considered failure patterns. One
can notice the occurrence of radial and longitudinal cracks during simulation as presented in Figure 1.49.
Comparison of the ice pressure distribution obtained by finite element computations with the measured
ice pressure distribution according to Sand showed similar trends.
49
Figure 1.48 – Transformation of state method used by Sand B. [63]
Figure 1.49 – Calculated failure pattern and distribution of ice pressure obtained by Sand (2008) [63]
Despite the widespread use of the finite element method, there is experience using other numerical
methods. An example of numerical experiments of ice-structure interaction using the discrete element
method is article of Liu et al. (2016) [45]. Authors presented ice load simulations for fixed and floating
structures, i.e. jack-up legs and the Kulluk floating drilling platform. They also compared the
computation time of a model using a computer processor Performing calculations using the GPU
significantly reduces the calculation time (up to 36 times).
Ice sheet in simulations was modeled using contacted spherical particles bonded together using a
50
parallel bond model, as the bonding disk (Figure 1.50). The maximum tensile and shear stresses acting
on the bonding disk are calculated based on the beam theory. In order to calibrate the relationship
between the bonding strength parameters assumed in the DEM model and realistic ice strength terms,
general ice sample tests of uniaxial compression test and 3-point bending test were simulated.
Figure 1.50 – Modeling the ice sheet using spherical bonded elements in Liu et al. (2016) [45]
Regarding vertical structures in this article the authors investigated the influence of the shape of
the structure and the mutual influence of several supports on the ice load. Example of ice-structure
interaction modeling with circular shape in this paper is presented in Figure 1.51. Results were compared
with calculations of the ice load using ISO 19906 recommendations.
Figure 1.51 – Snapshots of DEM simulations of an ice sheet moving through circular support
performed by Liu et al. (2016) [45]
The simulated ice force was less than 60% of the ISO 19906 formula value for the broadside
direction case. The average and maximum peak ice forces simulated by DEM are both larger for the
broadside case compared with the diagonal and round geometries but ISO 19906, which is based on
projected area, suggests otherwise. The DEM simulation results presented in this paper are promising,
but the accuracy and reliability of the tool still requires more validation studies and developments.
We can say that the simulation results in this paper showed poor convergence with analytical
methods for determining ice load. Along with the lack of proper consideration for the non-uniform
properties of ice and the high computational cost of calculations, the question of the applicability of
using the discrete element method in modeling ice effects remains open.
51
1.6 Section conclusions
The problem of assessing ice impacts on the offshore structures is becoming more global and
relevant in connection with the development of the energy complex and an increase in the rate of
hydrocarbon production in the Arctic zone. The impact of ice on structures is a multifactorial
phenomenon. The magnitude of ice loads depends not only on the climatic conditions, but also on the
characteristics of the ice itself, the interaction conditions at the contact and many other factors.
Sea ice is a complex material consisting of hard ice, brine, gas, the ratio of which depends on
temperature, texture, growth conditions, etc. Ways to describe the properties of ice have been developed
and improved since the second half of the 20th century. Thus, the studies of Russian and foreign authors
are based on a large number of field and laboratory experiments. Processing the results of these
experiments allowed to obtain analytical dependencies of the components of sea ice (fresh ice, water and
salts) on various properties. These properties vary greatly and depend on the area of research, methods
of conducting experiments and the age of the ice. Nevertheless, as a result of their work, researchers
managed to obtain many dependencies of ice parameters on various factors. As a rule, all parameters
show a strong dependence on ice temperature. The salinity, porosity and texture of ice also play a
significant role.
Comparison of methods for calculating the salinity of ice, brine volume, elastic modulus and
compressive strength is shown in Figure 1.52-Figure 1.55. In Figure 1.54 curve named “Gammon P. H
et al. (1983) + Weeks and Assur (1967)” means that values of elastic modulus of pure ice were selected
from Table 1.5, corrected using equation (1.43) and values for sea ice was calculated using equation
(1.48).
Figure 1.52 – The dependence of the salinity of ice on its thickness by various methods
52
Figure 1.53 – The dependence of the brine volume on ice temperature by various methods
Figure 1.54 – The dependence of the dynamic elastic modulus on ice temperature by various methods
53
Figure 1.55 – The dependence of the compressive strength of ice on strain rate
A significant scatter of data appears between different methods. The selection of appropriate
method should be done considering the testing approach that was used during experiments, age of
researches and applicability range of equations. The most important features of ice as material that must
be considering in numerical analysis are as follows:
- strong dependence of ice properties on temperature, salinity and its texture. Since during icestructure interaction process different conditions can appear at the same time, especially on contact with
structure, applied modeling methodology should have possibility to automatic recalculation of ice
properties based on current conditions (temperature, mean stress, strain rate, etc.);
- growth and texture forming conditions. Anisotropy of ice polycrystals has strong dependence on
ice texture. Thus, ice with random orientation of c-axes, for instance, granular ice, can be considered as
isotropic material. But in case of columnar ice, which has more or less ordered directions of c-axes,
polycrystal show clear anisotropic strength and elastic properties. So, the best way to describe this
feature is to consider columnar type of ice as transversely isotropic or orthotropic material;
- plastic-brittle transition of ice behavior based on interaction strain rates. Since this feature of ice
has not yet been taken into account by the authors earlier, this is an unresolved problem in modeling of
ice actions, which will ensure the scientific novelty of this thesis.
In case of mechanism of ice-structure interaction it can be said that li limit stress mode of
interaction is the most dangerous since has the highest forces. Crushing of ice is rapidly changing process
54
with occurrence of so-called high-pressure zones. To qualitatively take into account the influence of
these zones on the total ice effect, it is necessary to study the issue of modeling of ice destruction process
in details.
As for the worldwide experience of numerical modeling of ice impacts, a large number of papers
have been performed using various numerical methods and approaches. The most promising in my
opinion is the finite element method using the cohesive element method. In most cases, the authors do
not pay due attention to the selection of constitutive models suitable for describing the plastic behavior
and failure of ice. The considered papers of various authors do not take into account such features of ice,
as the dependence of its properties on temperature, strain rate, pressure-hardening (except Sand [63])
and texture. In this regard, it is difficult to talk about the correspondence of the results of numerical
models to real conditions.
Of great practical importance are the tools for optimizing the parameters of a material model for
specific conditions (for example, for experimental data). Solving the inverse problem, pre-setting some
restrictions, it is possible to numerically determine some empirical parameters of the material model.
55
2 Formulation of the problem and methodology of numerical modeling of ice impacts
2.1 Formulation of the problem
Offshore facilities on the continental shelf are unique facilities of increased responsibility, which
requires appropriate careful assessment of the effects of external factors. Proper design of such objects
ensures the reliability of structures and helps to prevent emergencies that may arise as a result of the
destruction of the structure.
Ice impacts are one of the most complex factors evaluated during design. This is directly related
to the spatio-temporal variability of the ice cover, the great variability of the physicomechanical
properties of ice, the influence of other environmental factors on the ice formations, such as waves,
currents, wind, etc. Therefore, it can be said that ice loads have stochastic nature, and load estimation
methods used today in the design of offshore structures use a number of significant assumptions that can
distort the real picture that arises in reality.
Numerical modeling methods can to some extent eliminate the difficulties that arise when taking
into account the random nature of the values of ice parameters, which significantly increases the
reliability of the magnitude of the ice load. In addition, as already mentioned, the destruction of ice
formations during their interaction with vertical structures depends on many factors (Section 1) and has
a strong influence on the formation of ice load. Thus, the assessment of the interaction of ice formations
with vertical structures using numerical modeling considering the process of ice destruction will allow
us to avoid several assumptions that affect the result.
In view of the above, to carry out a numerical simulation of the interaction of ice formations with
vertical structures, taking into account the destruction of ice, it is necessary to:
- analyze the methods for calculating ice parameters and determine their dependence on
temperature, salinity, strain rate, etc. (section 1);
- perform analysis and selection of constitutive models of plasticity and fracture, taking into
account their applicability for modeling ice impacts;
- to establish a single principle for creating a numerical model, that includes certain considered
environmental parameters, boundary conditions, geometry of ice formations, etc.;
- to make an algorithm for creating a numerical model that takes into account the above features.
- to establish parameters used in constitutive model to obtain the nature and picture of interaction,
close to realistic.
The best way to verify the model would be to reproduce (simulate) field and laboratory
experiments on ice impacts with identical conditions, such as dimensions of ice formation, its speed,
temperature gradient, structure width and shape etc.
If the results of numerical simulation are consistent with the experimental results to an acceptable
56
extent, the next step will be to simulate the effects of ice with the structure under different scenarios.
The most interesting and investigated feature of the interaction process is the effect of the complex
stress-strain state of ice on the magnitude of the ice load, as well as the effect of ice destruction on it.
Taking into account the destruction of ice upon contact will also allow to consider the process of
occurrence of high-pressure zones described earlier.
The solution of the above problems will ultimately allow to develop a methodology and describe
an algorithm for assessing ice impacts in real engineering tasks using numerical modeling.
2.2 Methodology of numerical simulation of ice impacts
2.2.1 Numerical method selection
Today, there is a wide selection of numerical simulation methods used to solve many engineering
problems. Among the most famous, the following methods will be discussed as possible methods for
simulation of ice-structure interaction.
Discrete element method (DEM). This numerical method was developed by Cundall and
Strack (1978) [16] and mainly used for modeling granular media (such as soils), consisting of a large
number of particles. The simulated medium is represented in the form of a set of individual particles that
have their own properties (density, elastic modulus, velocity, etc.). Particles can be presented in various
shapes and sizes. DEM allows to simulate the interaction between all particles, but usually it has
significant computational cost.
In ice problems, the method was mainly used to simulate the interaction of ice formations with
each other or with structures on a large scale as presented in Figure 2.1. In this case, ice formations
(fields, glaciers, icebergs) are modeled as separate particles interacting with each other. Moreover, in
most cases, it is assumed that the formations are indestructible. For example, the studies of Kim J-H.
and Kim Y. (2019) [37] and Richard and McKenna (2013) [57] can be distinguished.
Figure 2.1 – Simulation of the movement of vessel through discontinuous ice [57]
57
DEM has also been used to model the interaction of ice and structures on a more detailed scale. In
this case, the ice field is represented as a set of particles located in several layers as presented in Figure
1.51. The interaction between particles allows one to take into account the destruction, ridging and
moving of ice on the structure, which is very useful in modeling the actions of ice on inclined structures.
For example, the studies of Lu et al. (2012) [46] and Liu et al. (2016) [45] that was described earlier
(subsection 1.5) can be distinguished. However, this method has several significant disadvantages:
- the particles themselves are almost non-deformable and cannot be destroyed;
- the destruction of the ice field can only occur between particles, which imposes restrictions on
the path of conditional cracks;
- this method requires greater computing power compared with other considered numerical
methods;
- DEM is not widely used among software systems if compare with FEM and is more suitable for
modeling granular media than for solving a continuous medium problem.
Finite element method (FEM). To model the actions of ice formations, the finite element method
is most often used. This method has received the greatest distribution among methods of numerical
modeling and is actively developing today, because it has a wide range of advantages. The method is
based on the concept of partition of unity and widely used to solve problems of solid mechanics, heat
transfer, hydrodynamics and electrodynamics.
Using this method, all bodies are represented as they are, close to real geometry. Then the bodies
are discretized into finite elements – a mesh of each body is created. Finding a solution at finite element
nodes leads to solving a system of algebraic equations. The smaller the size of the finite elements (i.e.,
the more elements in the model), leads to a greater number of algebraic equations solved by the solver
and the greater the computational cost.
When modeling the ice-structure interaction using the finite element method, the ice formation is
often represented as a solid body. Usually the following two methods are used to simulate the destruction
of ice: cohesive element method and element erosion technique.
Cohesive element method. In this method, the simulated body is discretized using conventional
bulk elements, as well as special cohesive elements of zero thickness as presented in Figure 2.2. Bulk
elements are used to model the process of ice deformation, and destruction occurs due to the loss of
strength of the cohesive elements. Thus, it is assumed that the cohesive elements are possible paths for
the formation and development of cracks. The behavior of the cohesive elements is subject to predefined
traction-deformation rules. This method was widely used to model ice impacts on structures with taking
into account the ice destruction process. Among the studies considered before, the following can be
distinguished: Gürtner et al. (2010) [29], Li Liang (2014) [44], Salganik (2014) [62]. Among the
disadvantages of this method, one can distinguish the high computational cost, as well as the difficulty
58
of using in the case of complex geometry of ice formation.
Figure 2.2 – Mesh when using the cohesive element method
Element erosion technique. This method is also widely used to model ice-structure interaction
problems but is easier to use than the previous cohesive element method. In this case, the ice formation
is modeled as a completely continuous body using solid elements. To describe the behavior of ice under
loads, special material models are used, which describe:
- the yield criteria (yield surface) and flow rule;
- the damage initiation criteria;
- the damage evolution law.
The life cycle of the finite element in this method is as follows. First, the behavior of the material
completely depends on a predetermined constitutive model, then when a certain failure criterion is
reached, the destruction of the material begins according to the damage evolution law. When the degree
of damage reaches 100%, the material is considered destroyed and the finite element is removed from
the model. The following studies of several authors described earlier use this method: Sand (2008) [63],
Politko V.A. (2018) [56]. This method has high efficiency from the point of view of computation time,
however, a significant disadvantage is the strong dependence on the mesh size (the size of the cracks is
essentially connected with the size of the finite elements).
Moreover, when solving ice problems, two different formulations of finite elements can be used.
The Lagrangian-based finite element simulations is one of the most frequently used continuum
mechanics. This method relates the mesh and the material, which might cause excessive distortion and,
hence, numerical instability problems in which the material is deformed greatly. The ALE method
combines the methods of Lagrange and Euler (in the latter case, the mesh is fixed in space, and the
material flows through it). It takes advantage of both integration methods: the computational efficiency
of the Lagrange method and the ability to allow a significant deformation of the Euler method.
Extended Finite Element Method (XFEM). This method is an extension of the conventional
finite element method. It can model cracks and other discontinuities by enriching the degrees of freedom
in the model with additional displacement functions that account for the jump in displacements across
the discontinuity. In this way, cracks modeled by this method are independent of the mesh (mesh-free),
so the method is widely used in the fracture analysis. Nevertheless, it is practically not used to simulate
59
the interaction of ice with vertical structures, due to the difficulties in determining the position of the
initial crack and a small degree of integration into CAE software.
As an example of the application of the extended finite element method to solve ice problems, the
study of Lu et al. (2018) [47] in which authors performed the numerical modeling of ice destruction
during interaction of ice floe with icebreaker as presented in Figure 2.3.
Figure 2.3 – Crack path within the selected ice floe before and after the fracture in both the simulation
and reality according to Lu et al. (2018) [47]
After analyzing the considered methods of numerical modeling, it was decided to use the finite
element method in combination with cohesion element method. This method has the following
advantages compared to other methods:
- wide distribution and integration into a large number of software systems;
- relative ease of use and setting model parameters;
- clear and understandable interpretation of calculation results;
- good applicability for solving problems of solid mechanics;
- good agreement of results with field data (according to several authors).
A study of the effect of various material parameters on ice load using CEM was performed by
Feng et al. [22]. According to this study it was found that the structural response is highly sensitive to
the change of fracture energy rather than the cohesive strength used for defining the curve of cohesive
elements. The change of the shape or the initial slope of the curve has little effect on results. Also, the
stress-strain curve of the bulk elements has remarkable influence on simulations.
60
The solution comes down to solving the main equation of motion of the finite element method,
which has the following form:
[𝑀]{𝑥̈ } + [𝐶]{𝑥̇ } + [𝐾]{𝑥} = {𝐹},
(2.1)
where [𝑀], [𝐶] and [𝐾] – global mass, damping and stiffness matrixes respectively;
{𝑥̈ }, {𝑥̇ } and {𝑥} – nodal acceleration, velocity and displacement vectors respectively;
{𝐹} – applied load vector.
To perform numerical experiments, the SIMULIA Abaqus software package for finite element
analysis will be used. SIMULIA Abaqus – universal software system of finite element (FEM) analysis,
existing and developing over the past 40 years, is quite popular with specialists in the field of computeraided engineering (CAE) and linear and non-linear, stationary and non-stationary problems of mechanics
of a deformable solid and structural mechanics (including non-stationary geometrically and physically
non-linear contact problems interaction of structural elements). Also, it can be used for solving the
problems of fluid dynamics, heat transfer, electrodynamics, topology optimization and others.
2.2.2 Integration method selection
To simulate dynamic processes, including the impact of ice formations on a structure, two methods
for integration of motion are used: explicit and implicit methods.
In practice, implicit methods are usually used to simulate relatively long-running processes,
ranging from a few tenths of a second to several years (for example, the calculation of building structures
considering creep). Implicit methods are reduced to a series of solutions of quasistatic problems with
time-dependent loads. The time step may not be very small, since at each step a solution of the system
of equations is performed and balancing iterations associated with matrix operations are carried out.
Explicit methods are methods for solving equations of motion that are not related to solving
systems of equations, but using recurrence relations that express displacements, velocities, and
accelerations at current step through their values in the previous steps. To obtain a stable solution, the
time step in the calculations is very small compared to implicit methods. Such small steps allow to
calculate taking into account all the non-linearities and to thoroughly track the behavior of structures.
In general, the choice of integration method is based on the time duration of the simulated process,
as shown in the Figure 2.4.
Despite the fact that the process of ice interaction with structures is relatively lengthy, for the
correct consideration of the nonlinear properties of ice in the process of deformation, it is necessary to
use small time steps. In this regard, in this thesis the explicit integration method will be used for a
more detailed and quick assessment of the behavior of the ice during deformation.
61
Figure 2.4 – The application boundaries of implicit and explicit methods [6]
Explicit dynamics procedure solves every problem as a wave propagation problem. Out-of-balance
forces are propagated as stress waves between neighboring elements. As mentioned earlier, the time
increment in explicit dynamics procedure is very small to achieve system stability. A bounded solution
is obtained only when the time increment ∆𝑡 is less than the stable time increment ∆𝑡𝑚𝑖𝑛 . The stable
time increment is the minimum time that a dilatational (i.e., pressure) wave takes to move across any
element in the model. Its value for each element in the model can be calculated as follows:
∆𝑡𝑚𝑖𝑛 =
𝐿𝑒
,
𝑐𝑑
(2.2)
where 𝐿𝑒 – characteristic length of element;
𝑐𝑑 – the dilatational wave speed:
𝐸
𝑐𝑑 = √ ,
𝜌
(2.3)
where 𝐸 – is the Young’s modulus of material;
𝜌 – current material density.
Thus, decreasing element dimensions, increasing material stiffness, decreasing material
compressibility and decreasing material density can reduce the stable time increment and increase the
computational cost of solution.
Explicit method uses a central difference rule to integrate the equations of motion explicitly
through time, using the kinematic conditions at one increment to calculate the kinematic conditions at
the next increment. Initially, the solution comes down to defining accelerations at the beginning of the
current increment that can be expressed from equation (2.1) as follows
{𝑥̈ } = [𝑀]−1 ({𝐹} − ([𝐶]{𝑥̇ } + [𝐾]{𝑥}))
(2.4)
Thus, the acceleration of any node is determined completely by its mass and the net force acting
on it. Then the accelerations are integrated through time using the central difference rule, which
calculates the change in velocity assuming that the acceleration is constant. This change in velocity is
62
added to the velocity from the middle of the previous increment to determine the velocities at the middle
of the current increment:
𝑥̇
∆𝑡
(𝑡+ 2 )
= 𝑥̇
∆𝑡
(𝑡− 2 )
+
(∆𝑡(𝑡+∆𝑡) + ∆𝑡(𝑡) )
𝑥̈ (𝑡)
2
(2.5)
The velocities are integrated through time and added to the displacements at the beginning of the
increment to determine the displacements at the end of the increment:
𝑥(𝑡+∆𝑡) = 𝑥(𝑡) + ∆𝑡(𝑡+∆𝑡) 𝑥̇
(𝑡+
∆𝑡
)
2
(2.6)
The solution process of this method of integration of the equations of dynamics is shown in Figure
2.5 and briefly can be explained following steps:
a) the motion of the node points produces deformation in the elements of the mesh;
b) the deformation results in a change in volume (hence density) of the material in each element;
c) the rate of deformation is used to derive material strain rates using various element formulations;
d) constitutive laws take the material strain rates and derive resultant material stresses;
e)the material stresses are transformed back into nodal forces using various element formulations;
f) external nodal forces are computed from boundary conditions, loads and contact (body
interaction);
g) the nodal forces are divided by nodal mass to produce nodal accelerations;
h) the accelerations are integrated explicitly in time to produce new nodal velocities;
i) the nodal velocities are integrated explicitly in time to produce new nodal positions;
g) the solution process (cycle) is repeated until a user defined time is reached.
Figure 2.5 – Solution process of explicit method [7]
63
2.2.3 General modeling techniques
Geometry of ice formations and structures. Of greatest interest in terms of the shape of vertical
structures are single-cylinder cylindrical supports, due to their wide distribution in the construction of
offshore structures. Structures of other forms in the plan are also available, but much fewer are built. In
this regard, in this dissertation, the impacts of ice formations on circular structures will be investigated.
Among the various forms of ice formations, according to the methods of calculating ice load, they
distinguish level-ice sheets, hummocks and icebergs. The most common case of impact is the impact of
level ice. Its thickness can reach several meters and, in the case of multi-year ice, have a very
heterogeneous structure in thickness (Figure 1.4). The dependence of the ice load on the thickness in
real conditions is indirect, due to the variable properties through the thickness. In modern analytical
calculation methods, today, the dependence of the load on the ratio 𝑑/ℎ is considered, where d – width
(diameter) of the structure and ℎ is thickness of level-ice sheet. For offshore structures, this ratio usually
has a value of more than 5. Lower ratios (with small diameters) are more likely to correspond to various
port structures (for example, piers and berths). Thus, cases of ice-structure interactions with ratios of
𝒅/𝒉 greater than 5 will be considered.
As mentioned earlier, the case of interaction with limit stress mode is the most dangerous and
therefore of great interest. Thus, the numerical modeling will be performed under the assumption of the
infinite ice field conditions. Shape of ice sheet in plan is applied as rectangular. To reduce the influence
of boundary conditions, the dimensions of the ice field with respect to the structure diameter 𝒅 are
taken as follows:
- ice field width in the direction of movement – 6𝑑;
- width of front face – 12𝑑.
These dimensions are assigned based on the assumption that stresses propagate over a width of 4𝑑
(interaction zone). This fact has been investigated by researchers earlier [38] and is adopted in the current
national codes [68].
Applied physical and mechanical properties of ice. The importance of taking into account the
anisotropic properties of ice over the thickness of the ice field is described in Section 1. To consider the
heterogeneous structure of the ice field, the simulated ice body is divided into 2 connected bodies
with a thickness ratio of 1:3 according to SP 38.13330.2018 [68]. This division of the ice field is shown
in the Figure 2.6.
It is assumed that the upper body is granular ice layer with random c-axes orientations and
therefore with isotropic properties. The lower body is taken as a columnar layer with orthotropic
properties. For this case the most convenient way to describe ice properties is to accept the fact that the
ice field has S2 texture type (Table 1.6) at which c-axes are randomly oriented in 𝑋1 –𝑋2 plane (horizontal
64
plane) and X3 is parallel to the direction of growth of the ice field.
Figure 2.6 – Geometric dimensions and texture assignments of modeled ice sheet (vertical dimension
is scaled by 5 times)
The temperature gradient in ice is assumed linear. The initial temperatures of the ice field are
generated by a preliminary steady-state thermal analysis. The temperature calculation results are then
integrated into the main model. The boundary temperatures are: on the top face of the ice field – ambient
temperature; on the bottom face of the ice field – the freezing temperature of sea water, determined
according to SP 38.13330.2018 as follows [68]:
𝑡𝑏 = −0,057𝑠𝑤 ,
(2.7)
where 𝑠𝑤 – water salinity, ‰.
All equations of physical and mechanical properties of ice, which will be used in this thesis, are
presented in Table 2.1. Thus, the initial independent physical properties are temperature at the upper
boundary of the ice field, water salinity and ice thickness. Other physical and mechanical properties
depend on each other.
Meshing and boundary conditions. Since the cohesive element method shows significant
dependence on mesh size, the appropriate mesh settings must be applied. The main work in the
interaction of ice on the structures is done by the material near the contact. This region is discretized into
the smallest finite elements and has dimensions of about 2𝑑. Further, the finite elements of the ice field
will gradually increase in size up to the boundaries of the simulated ice body. To optimize the time of
the model creation process, as well as reduce the distortion of finite elements, their number in thickness
is taken equal for the whole body.
65
Table 2.1 – Equations for determining the physical and mechanical properties of ice
Ice properties
Ice salinity
Equation
𝑠𝑖 = 4,606 + (
91,603
)
ℎ𝑖
𝜌𝑖 = 0,917 − 1,403 ∙ 10−4 𝑇
Ice density
𝜌𝑖 ∙ 𝐹1 (𝑇)
𝜌𝑠𝑖 = (1 − 𝜈𝑎 )
𝐹1 (𝑇) − 𝜌𝑖 ∙ 𝑠𝑖 ∙ 𝐹2 (𝑇)
1 + (𝛽𝑇 + 𝑙)
𝛼𝑇 + 𝑘
𝛽𝑇 + 𝑙
𝑠𝑏 =
1 + 𝛽𝑇 + 𝑙
𝜌𝑠𝑖 𝑠𝑖
𝜈𝑏 =
𝐹1 (𝑇)
𝑚𝑏 = 𝑠𝑖
Brine salinity,
mass and
volume
𝐸𝑠𝑖 = 𝐸𝑖 (1 − 𝜈𝑏 )4
Elastic
properties
Equation number
Reference
Additional information
(1.35)
Kovacs (1996) [40]
Used in ISO/DIS 19906:2019
Will be used for ice thickness up to 2 m
Cox and Weeks
(1983) [14]
-
(1.39)
(1.40)
(1.16)
Nazintsev and Panov
(2000) [51]
(1.18)
Cox and Weeks
(1983) [14]
Weeks and Assur
(1967) [81]
(1.21)
(1.48)
Used in SP 38.13330.2018
-
Used in ISO/DIS 19906:2010
𝐸(𝑇) = 𝐸(𝑇𝑟 ) ∙ (1 − 𝑎(𝑇 − 𝑇𝑟 ))
𝐺(𝑇) = 𝐺(𝑇𝑟 ) ∙ (1 − 𝑎(𝑇 − 𝑇𝑟 ))
(1.43)
Gammon et al. (1983)
[25]
Based on measurements using Brillouin
spectroscopy
Timco and Frederking
(1990) [75]
Used in ISO/DIS 19906, ANSI API RP 2N2015, CSA/S471-04
𝜇𝑖,𝑠 (𝑇) = 𝜇𝑖,𝑠 (𝑇𝑟 ) ∙ (1 − 𝑎(𝑇 − 𝑇𝑟 ))
𝜈𝑡
)
270
(1.57)
𝜈𝑡
)
200
(1.58)
𝜈𝑡
)
280
(1.59)
𝜎𝑐 = 37(𝜀̇)0,22 (1 − √
Compressive
strength
𝜎𝑐 = 160(𝜀̇)0,22 (1 − √
𝜎𝑐 = 49(𝜀̇)0,22 (1 − √
66
The effect of the mesh pattern on the result will be later explored in this thesis. To simulate the
conditions of an infinite ice field, all the faces of the body except the frontal one are constrained from
movements in the direction orthogonal to the direction of motion. The initial velocity is also assumed
constant on these faces. Destroyed ice in the process of interaction continues to affect the structure.
Crushed ice increases the contact area, and also adds some vertical component acting on the ice field. In
this regard, the forces of gravity and buoyance cannot be neglected. Gravity is taken into account very
simply, by adding the appropriate load object to the bulk elements through the interface of the Abaqus
software. Buoyancy can be modeled in several ways.
The first is a Coupled Eulerian-Lagrangian (CEL) approach. In this case the two formulations of
finite elements (Eulerian and Lagrangian) are used simultaneously. Ice field and structure are modeled
by finite elements in Lagrangian formulation (material and mesh are interconnected). Water medium is
modeled by Eulerian domain in which the mesh doesn’t move during the analysis and water material
“flow” through the mesh. For this approach water material is usually modeled using equations of state
(EOS) that represents the hydrodynamic behaviour. Since this approach significantly increases the
computational cost of calculations, the second method of modeling buoyancy will be used.
The main idea of the second method is to reproduce the water pressure acting on the external faces
of bulk elements. If we assume that the vertical z-axis in the model is directed up, then buoyancy pressure
𝑝, Pa, in the node of submerged finite element can be calculated as follows:
𝑝 = 𝜌𝑤 𝑔(𝑧0 − 𝑧),
(2.8)
where 𝜌𝑤 – density of sea water, kg/m3;
𝑔 – gravitational acceleration, m/s2; 𝑔 ≈ 9,81 𝑚/𝑠 2 ;
𝑧0 – coordinate of water surface in vertical direction (z-axis), m;
𝑧 – coordinate of node in vertical direction, m;
In the toolkit of the Abaqus software package, there is no way to add a pressure depending on the
distance, the value of which will be updated every time increment. To add such pressure load the
following subroutine was created in the Fortran programming language, which was then initialized
at the beginning of each calculation time increment:
subroutine vdload (nBlock, ndim, stepTime, totalTime, amplitude,
& curCoords, velocity, dirCos, jltyp, sname, value)
! General subroutine settings
include 'vaba_param.inc'
dimension curCoords(nBlock,ndim), velocity(nBlock,ndim),
1 dirCos(nBlock,ndim,ndim), value(nBlock)
character*80 sname
! User-defined parameters
parameter (WaterDens = 1020, Gravity = 9.81, Zmax = 0, IceH = 1,
*
IceDens = 930)
! WaterDens - approximate water density
! Gravity - gravity acceleration
! Zmax - z-coord (in vertical direction) of ice bottom
67
! Ice dens - approximate ice density (average)
parameter (iX = 1,
*
iY = 2,
*
iZ = 3)
! Subroutine main
! Repeat every time increment
do i = 1, nBlock
if(curCoords(i,iZ) .Le. Zmax) then
Pressure = WaterDens*Gravity*(IceH*IceDens/WaterDens+Zmax
*
- curCoords(i,iZ)) ! Calculated buoyancy pressure
value(i) = Pressure
else
value(i) = 0.d0
endif
end do
return
end
2.2.4 Selection of constitutive models for ice material
Sea ice, as described earlier, is an extremely complex material. It can exhibit both plastic and
brittle properties. Moreover, depending on the structure, ice formations can be considered as isotropic
or orthotropic bodies. In this thesis, granular ice located mainly in the upper layers are assumed to be
isotropic, since, as mentioned above (subsection 1.3), it has a random orientation of single crystals. In
the case of columnar ice, the principle of isotropy cannot be applied, since the properties of such ice are
very different in different directions. A good solution would be to consider columnar ice polycrystals as
an orthotropic material.
To describe complex material behavior constitutive model are used. Constitutive model of material
is a set of mathematical relationships which describes the material behavior under different loading
conditions. Today there are a huge number of different models of materials that mathematically describe
their various features. To model the actions of ice on vertical offshore structures, it is necessary to
describe the behavior of the ice at all its stages up to destruction. For this, elasticity models, plasticity
models, as well as the damage initiation (failure) criteria and damage evolution laws, will be considered.
Elasticity model. As basic parameters for modeling, elastic properties are used. These properties
can be defined via the corresponding elasticity models. In this thesis, isotropic and orthotropic
elasticity is used to model ice material behavior where the stress versus strain relationship is linear and
the loading is kept within the elastic range. The models follow Hooke's law (equation 1.42), so assumes
that the stretch or compression remains in the elastic range of the material and the model will return to
its original shape after unloading. The elastic properties of polycrystals of isotropic and orthotropic ice
are presented in Table 1.5 and Table 1.6, respectively. The temperature dependence of the characteristics
is determined by the equation (1.43).
68
In case of cohesive elements, the constitutive response is defined in terms of traction versusseparation law. It assumes a linear elastic traction-separation law prior to damage. The elastic behavior
is written in terms of an elastic constitutive matrix that relates the nominal stresses to the nominal strains
across the interface. The nominal stresses are the force components divided by the original area at each
integration point, while the nominal strains are the separations divided by the original thickness at each
integration point. The default value of the original constitutive thickness equal to 1.0 will be used, which
ensures that the nominal strain is equal to the separation. The elastic behavior can then be written as [71]
𝑡𝑛
𝐸𝑛𝑛
𝑡 = { 𝑡𝑠 } = [ 𝐸𝑛𝑠
𝑡𝑡
𝐸𝑛𝑡
𝐸𝑛𝑠
𝐸𝑠𝑠
𝐸𝑠𝑡
𝐸𝑛𝑡 𝜀𝑛
𝐸𝑠𝑡 ] { 𝜀𝑠 } ,
𝐸𝑡𝑡 𝜀𝑡
(2.9)
where 𝑡𝑛 – normal component of nominal traction stress vector, Pa;
𝑡𝑠 and 𝑡𝑡 – shear components of nominal traction stress vector, Pa;
𝐸𝑖𝑗 – corresponding elastic (or shear) modulus, Pa;
𝜀𝑛 , 𝜀𝑠 , 𝜀𝑡 – normal and two shear components of strain vector respectively:
𝜀𝑖 =
𝛿𝑖
,
𝑇0
(2.10)
where 𝛿𝑖 – corresponding separations, m;
𝑇0 – original constitutive thickness of cohesive element; since 𝑇0 = 1 the strains and separations
are equal.
However, a simplified formulation can be used in which the behavior of the cohesive elements in
the normal and tangent directions is uncoupled. In this case each traction component 𝑡𝑖 depends only on
its conjugate nominal strain 𝜀𝑖 :
𝑡𝑛
𝐸𝑛𝑛
𝑡 = { 𝑡𝑠 } = [ 0
𝑡𝑡
0
0
𝐸𝑠𝑠
0
0 𝜀𝑛
0 ] { 𝜀𝑠 } ,
𝐸𝑡𝑡 𝜀𝑡
(2.11)
Plasticity model. To describe the plastic properties of ice, it is necessary to choose a plasticity
model. Typically, plasticity models consist of three main components [4]:
- the yield criterion that defines the material state at the transition from elastic to plastic behavior;
- the flow rule that determines the increment in plastic strain from the increment in load;
- the hardening (or softening) rule that gives the evolution in the yield criterion during plastic
deformation.
The following yield theories with corresponding materials models are most widely used in
structural analysis: von Mises yield theory, Tresca theory, Mohr-Coulomb theory and Drucker-Prager
theory. High usability of these models is associated with the ease of use – the parameters needed to
determine the yield surface can be obtained from standard material tests (yield strength, internal friction
angle, adhesion).
69
Von Mises yield criterion is isotropic and independent of hydrostatic pressure, which can limit its
applicability to microstructured materials and materials that exhibit plastic dilatation such as steels. This
material model as mentioned before was used by such researchers as Gürtner et al. (2009) [30], Gürtner
et al. (2010) [29], Hilding et al. (2011) [31], Hilding et al. (2012) [32] and Wang et al. (2019) [79]. It
includes an associated flow rule. Yield surface in axes of principal stresses is presented in Figure 2.7, a.
The von Mises yield criterion is [4]:
𝐹 = 𝜎𝑒 − 𝜎𝑦 = 0,
(2.12)
where 𝜎𝑦 – yield strength of material;
𝜎𝑒 is the von Mises equivalent stress, also known as the von Mises equivalent stress:
𝜎𝑒 = √3𝐽2 ,
(2.13)
(𝜎11 − 𝜎22 )2 + (𝜎22 − 𝜎33 )2 + (𝜎33 − 𝜎11 )2
2
2
2
𝐽2 =
+ 𝜎12
+ 𝜎23
+ 𝜎31
,
6
(2.14)
𝐽2 – second deviatoric stress invariant:
where 𝜎11 , 𝜎22 , 𝜎33 , 𝜎12 , 𝜎23 , 𝜎31 are stress components of Cauchy stress tensor 𝝈:
𝜎11 𝜎12 𝜎13
𝜎
𝝈 = [ 21 𝜎22 𝜎23 ].
𝜎31 𝜎32 𝜎33
The Tresca theory suggests that the material goes into a plastic state when the maximum shear
stress 𝜏𝑚𝑎𝑥 reaches the shear yield strength 𝜎𝑠ℎ𝑒𝑎𝑟 . Yield surface in axes of principal stresses is presented
in Figure 2.7, a. The Tresca yield criterion is [4]:
𝐹 = 𝜏𝑚𝑎𝑥 − 𝜎𝑠ℎ𝑒𝑎𝑟 = 0,
(2.15)
The above criteria do not consider the effect of hydrostatic stress on the condition of the transition
of the material into plastic behavior. Some materials show the influence of the hydrostatic stress
component on yielding onset: the higher the hydrostatic stress (confinement pressure) the higher the
yield strength. The two simplest models that take this feature into account are the Mohr-Coulomb and
the Drucker-Prager models.
The Mohr-Coulomb material model is used to represent the behavior of aggregate materials.
Aggregate materials such as soil, rock and concrete begin to plastically deform when the shear stress
exceeds the internal friction resistance between the material particles. The friction resistance is a function
of the normal force between the particles. The model defines yielding when the combination of pressure
and shear stress reaches the cohesion of the material particles. Mohr-Coulomb yield criterion was used
as fracture criterion by Politko (2018) [56]. Yield surface in axes of principal stresses is presented in
Figure 2.7, b. Yielding occurs when the shear stress 𝜏 on any plane in the material reaches this criterion:
𝜏 = 𝑐 − 𝜎𝑚 tan 𝜑,
where 𝑐 – material cohesion;
𝜑 – angle of internal friction;
70
(2.16)
𝜎𝑚 – hydrostatic stress:
𝜎𝑚 =
𝜎11 + 𝜎22 + 𝜎33
3
(2.17)
The classic Drucker-Prager model is used for modeling the granular materials such as soils, rock,
concrete and uses the outer cone approximation to the Mohr-Coulomb law. This yield surface is a
circular cone with the material parameters chosen such that it corresponds to the outer apices of the
hexagonal Mohr-Coulomb yield surface as presented in Figure 2.7, b. Yield criterion of classic DruckerPrager model is [5]:
𝐹 = 3𝛽𝜎𝑚 + √𝐽2 − 𝜎𝑦 = 0,
(2.18)
where 𝛽 and 𝜎𝑦 – material model parameters;
𝐽2 – second deviatoric stress invariant expressed in equation (2.14);
𝜎𝑚 – hydrostatic stress expressed in equation (2.17).
Figure 2.7 – Yield surfaces in axes of principal stresses: a) von Mises and Tresca; b) Mohr-Coulomb
and classic Drucker-Prager [56]
Among the other plasticity models that were used in different papers the following can be
distinguished: Extended Drucker-Prager models, Hill plasticity, Johnson-Cook, Orthotropic Yield
model, Horrigmoe and Zeng’s model and others. Some description of these models and discussion about
their applicability for modeling ice problems are described in [48].
For isotropic ice Linear Drucker-Prager criteria will be used for modeling the yielding of
pressure-dependent granular ice with different values of tension and compression yield stresses. Yield
surface shown in Figure 2.8 provides for a possibly noncircular form in the deviatoric plane (π-plane) to
match different yield values in triaxial tension and compression, associated inelastic flow in the
deviatoric plane, and separate dilation and friction angles. The linear Drucker-Prager criterion is written
as
𝐹 = 𝑡 − 𝑝 tan 𝜑 − 𝑐 = 0,
(2.19)
𝑝 = −𝜎𝑚
(2.20)
where 𝑝 – pressure stress:
71
𝜑 – the slope of the linear yield surface in the p–t stress plane and is commonly referred to as the
friction angle of the material;
𝑡 – value determined by the following equation:
1
1
1 𝐽3 3
𝑡 = 𝜎𝑒 (1 + − (1 − ) ( ) ) ,
2
𝐾
𝐾 𝜎𝑒
(2.21)
where 𝐾 – the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression and,
thus, controls the dependence of the yield surface on the value of the intermediate principal stress;
𝐽3 – third deviatoric stress invariant:
2 3 1
𝐼 − 𝐼 𝐼 + 𝐼3 ,
27 1 3 1 2
(2.22)
where 𝐼1 , 𝐼2 and 𝐼3 – first, second and third stress invariants, respectively:
𝜎22
𝐼2 = |𝜎
32
𝐼1 = 𝜎11 + 𝜎22 + 𝜎33 ,
𝜎23
𝜎11 𝜎13
𝜎11
|
+
|
|
+
|
𝜎33
𝜎31 𝜎33
𝜎21
𝜎11 𝜎12 𝜎13
𝜎
𝐼3 = | 21 𝜎22 𝜎23 |
𝜎31 𝜎32 𝜎33
(2.23)
𝜎12
𝜎22 | ,
(2.24)
(2.25)
Figure 2.8 – Yield surface of linear Drucker-Prager criterion in the meridional plane
The input parameters when using Simulia ABAQUS software are the angle of internal friction 𝜑,
as well as the stress 𝜎̅0 , which determine the initial shape of the yield surface (at zero values of plastic
strains 𝜀̅𝑝𝑙 ). The yield strengths described above also depend on the strain rate and temperature.
In Simulia ABAQUS, the stresses 𝜎̅0 , which determine the initial yield surface, as well as the
equivalent stresses 𝜎̅, with the corresponding plastic strain values 𝜀̅𝑝𝑙 , which determine the hardening,
can be interpreted by the user as the yield strength under uniaxial compression, tension or pure shear. It
is obvious that the most convenient is the determination of hardening through uniaxial compression
stresses, since in most cases the data on the nonlinear nature of ice deformation were obtained by the
authors through uniaxial compression tests. On the other hand, for brittle materials, including ice, the
cohesion 𝑐 can be interpreted as the yield stress on pure shear, but when determining the hardening
through compression stresses, the cohesion value 𝑐 will be automatically calculated in the model
72
depending on the angle of internal friction 𝛽, defined by the user, according to the following
equation [72]:
1
𝑑 = (1 − tan 𝛽) 𝜎̅
3
(2.26)
In addition, the angle of internal friction β for brittle materials is a conditional characteristic, and
for a material such as ice, the selection of the optimal value is a difficult task. In this regard, the following
algorithm for modeling ice using the linear Drucker-Prager model in the Simulia ABAQUS software
package will be used [48]:
1) Determination of values of the ultimate strength of ice for uniaxial compression 𝜎𝑐 and shear
𝜎𝜏 , corresponding to several values of temperature 𝑡 and strain rate 𝜀̅̇𝑝𝑙 .
2) Assuming that 𝜎̅ = 𝜎𝑐 , and 𝑐 = 𝜎𝜏 , the calculation of the angle of internal friction values 𝛽 by
the equation
𝛽 = tan−1 (
3 ∙ (𝜎𝑐 − 𝜎𝜏 )
)
𝜎𝑐
(2.27)
Thus, the shape of the yield surface will most accurately describe the strength characteristics of
the simulated ice. When the stress level in the volume element of the compressive strength is reached,
its deformation will continue as perfectly plastic.
Failure criterion and post-fracture behavior. To determine the conditions under which the
material is destroyed, the failure criteria are used. When using the damage evolution laws fracture criteria
are also called “damage initiation criteria”. Often, fracture criteria are represented as failure surfaces in
the space of principal stresses. In the structural analysis of structures, the behavior of materials after
destruction is usually not considered, since structures are designed from the condition of preventing their
destruction. However, in this thesis, the process of ice load formation is considered taking into account
ice destruction. In the simple case, it can be accepted that the ice is destroyed after failure. In this way,
instantaneous failure is modeled. Bur after the appearance of cracks and partial destruction, the ice is
still able to affect the structure up to complete destruction. The so-called damage evolution laws can be
used to describe the behavior of ice after failure. In this way, gradual failure (damage) is modeled.
Fracture energy-based damage models are commonly used for this purpose.
Many criteria determine the destruction of the material when the value of plastic deformation 𝜀̅𝑝𝑙
𝑝𝑙
reaches a critical value 𝜀̅𝑚𝑎𝑥 . In most cases, the fact is accepted that the destruction of the material occurs
when it reaches a yield surface (for example, when calculating steel structures). This is especially valid
when describing the fracture of brittle materials such as rock, concrete and ice. In this case, it can be
accepted that the condition of fracture onset is the appearance of plastic deformations (𝜀𝑝𝑙 > 0). These
criteria work well when describing the fracture of ductile materials such as metals. Estimation of optimal
value of critical strains is a very difficult task, since in experiments even among samples with similar
73
parameters, critical strains vary greatly However, in order to maintain the stability of the system and
prevent excessive distortion of the finite elements, the criterion of critical plastic strains will be applied
to bulk elements. A plastic strain value of 1 will be used as a criterion for the removal of elements for
the first modeling iterations. This value is several times higher than the possible real values for ice and,
accordingly, will not affect the ice-structure interaction process.
In subsequent numerical experiments, the main process of destruction of ice is modeled using
cohesive elements. So, the traction of cohesive elements is not infinite and have limits. In Abaqus
software there are several failure models that can be applied to cohesive elements. The initial response
of the cohesive element is assumed to be linear as discussed above. However, once a damage initiation
criterion is met, material damage occurs according to a damage evolution law also called tractionseparation law (Figure 2.9). The cohesive elements do not undergo damage under pure compression.
Figure 2.9 –Typical traction-separation response of cohesive elements [71]
Damage initiation refers to the beginning of degradation of the response of a material point. The
process of degradation begins when the stresses satisfy the following damage initiation criteria:
𝑚𝑎𝑥 {
〈𝑡𝑛 〉 𝑡𝑠 𝑡𝑡
, , } = 1,
𝑡𝑛0 𝑡𝑠0 𝑡𝑡0
(2.28)
where 𝑡𝑛0 , 𝑡𝑠0 and 𝑡𝑡0 – maximum values of corresponding tractions.
The damage evolution law describes the rate at which the material stiffness is degraded once the
corresponding initiation criterion is reached. The stress components of the traction-separation model are
affected by the damage according to:
(1 − 𝐷)𝑡𝑛̅ ,
𝑡𝑛 = {
𝑡𝑛̅ ,
𝑡𝑛̅ ≥ 0
𝑡𝑛̅ < 0
(2.29)
𝑡𝑠 = (1 − 𝐷)𝑡𝑠̅ ,
(2.30)
𝑡𝑡 = (1 − 𝐷)𝑡𝑡̅ ,
(2.31)
where 𝑡𝑛̅ , 𝑡𝑠̅ and 𝑡𝑡̅ are the stress components predicted by the elastic traction-separation behavior for
the current strains without damage;
74
𝐷 – scalar damage variable represents the overall damage in the material and captures the
combined effects of all the active mechanisms. It initially has a value of 0. Then the value of 𝐷
monotonically evolves from 0 to 1 upon further loading after the initiation of damage.
In the case of an orthotropic material, individual strength limits for compression, tension and shear
in all three directions can be considered.
2.3 Section conclusions
Thus, the advantages and disadvantages of numerical methods were considered. The finite element
method with an approach to modeling ice fracture by the method of cohesive elements was chosen as a
tool for solving the tasks described above. Numerical modeling will be performed using the Simulia
ABAQUS software and for integration of general equation of motion (2.1) explicit method will be used.
The basic principles of creation a numerical model were also identified, as well as constitutive
models for modeling ice material. To use the above principles in numerical modeling, it is necessary to
follow the certain algorithm, starting with a description of the simulated situation (the size and shape of
the ice formation and structure, the speed of ice formation, temperature) and ending with the
interpretation of the calculation results (history of ice load, cyclic movement of the structure). Typical
scheme for performing numerical experiments is presented in Figure 2.10.
Figure 2.10 – Typical calculation algorithm
75
In the numerical model there is a joint work of bulk elements modeling the deformability of ice,
as well as cohesive elements modeling its destruction. The parameters of a model with such a complex
interaction require additional justification. Thereby, to analyze the applicability of the developed model
it is necessary to perform numerical experiments. The main goal of experiments is to study the developed
model of the ice-structure interaction by varying its parameters and assessing the quality of the
interaction picture. Subsequent simulations will be performed according to following plan:
1. Studying the influence of a mesh pattern;
2. Studying the influence of the size of bulk elements;
3. Selection of the optimal value of critical strains at which bulk elements are removed;
4. The study of correspondence of the plasticity response and the pattern of deformation of bulk
elements to real conditions;
5. Examination of the formation of a complex stress-strain state and its effect on the ice load value.
76
3 Numerical experiments on the impacts of ice formations on vertical structures and
development of the model
3.1 Assessment of the influence of mesh pattern
Since the solution of the tasks is performed by the numerical finite element method, it is necessary
to evaluate the effect of the finite element mesh on the modeling results. This assessment is especially
important when using the cohesive element method since the faces of bulk elements are possible crack
paths. There are several types of finite elements in the library of the Abaqus/Explicit module, with which
the volumetric body of an ice field can be discretized. They differ in geometric order, degrees of freedom,
the presence of hourglass control and other functions. Since connected finite elements can only be of
low geometric order and the ice field mesh must be conformal, quadratic bulk elements are not
considered further. Available cohesive and linear bulk finite elements with temperature degree of
freedom in the Abaqus/Explicit module are presented in Table 3.1.
Table 3.1 – Available types of finite elements for modeling of ice field
Code name
Description
Bulk elements
C3D4T
A 4-node thermally coupled tetrahedron, linear
displacement and temperature
C3D6T
A 6-node thermally coupled triangular prism, linear
displacement and temperature
C3D8T
An 8-node thermally coupled brick, trilinear
displacement and temperature
77
Picture
Continuation of Table 3.1
Cohesive elements
COH3D6
A 6-node three-dimensional cohesive element
COH3D8
An 8-node three-dimensional cohesive element
Thus, four mesh patterns will be compared: unstructured tetrahedral mesh using C3D4T and
COH3D6 elements, unstructured prism mesh using C3D6T, COH3D6 and COH3D8 elements,
structured hex mesh using C3D8T and COH3D8 elements, unstructured hex mesh using C3D8T and
COH3D8 elements. Finite element models for all 4 cases are presented in Figure 3.1. Minimum and
maximum mesh sizes are the same for all cases and equal to 0.25 m and 1.5 m respectively. Only four
“layers” of elements in thickness are created to reduce calculation time.
Figure 3.1 – Finite element models with different mesh patterns: a) tetrahedral mesh; b) prism mesh; c)
structured hex mesh; d) unstructured hex mesh
To perform this study, a model is created with the following parameters:
- diameter of the structure 𝑑 equal to 5 m;
78
- ice sheet dimensions 60x30x1 m;
- drift velocity 𝑣 equal to 0.5 m/s;
- water salinity 𝑠𝑤 equal to 20 ‰;
- air temperature 𝑡𝑎 equal to -20 °С;
- only isotropic ice and linear behavior will be considered.
Other physical and mechanical properties are determined according to the accepted methodology
(section 2). The time of simulation is assumed to be 5 seconds, so that the structure crashes into the ice
field to full width. As an evaluation criterion, the general picture of ice destruction, as well as the history
of ice load during the simulation, are used.
In Table 3.2 the information is presented about the components of the computer on which all
calculations are performed.
Unfortunately, in the current version of the ABAQUS software package it is impossible to use a
graphic card to perform calculations with explicit time integration. All calculations are performed on
CPU. General modeling information is presented in Table 3.3. Ice load history in the direction parallel
(x-axis) and orthogonal to the direction of movement (y-axis) is presented in Figure 3.2 and Figure 3.3.
Table 3.2 – Computer specifications
Component
General information
Motherboard
GigabyteX570 Aorus Ultra (rev. 1.0) PCI Express 4.0
Central processing unit
(CPU)
Random-access memory
(RAM)
Storage device
AMD Ryzen 9 3900x, 12 cores 4.1 GHz, 24 threads
DDR4 DIMM 16x4 GB (64GB total) 3200 MHz
Sabrent 1TB Rocket NVMe 4.0 Gen4 PCIe M.2 Internal SSD Extreme
Performance Solid State Drive (SSD);
Read speed 5000 MB/s; Write speed 4400 MB/s
Table 3.3 – General modeling information
Case
Number of finite
elements (bulk /
cohesive)
Total calculation time,
hours
Peak total ice force
along direction of
motion 𝐹𝑥𝑡𝑜𝑡 , MN
Tetrahedral mesh
124 261 / 236 075
48.63
(4.01 s out of 5)
-
Prism mesh
46 272 / 103 672
12.61
7.358
Unstructured hex mesh
26 312 / 71 918
4.88
9.785
Structured hex mesh
29 808 / 81 264
2.25
6.224
79
Figure 3.2 – Load history with different mesh patterns (along x-axis)
Figure 3.3 – Load history with different mesh patterns (along y-axis)
As we can see mesh patterns showed quite different results. First, let us consider the structured
hex mesh case. In my opinion this type of mesh showed the worst result. Peak force with extremely large
value of 16.676 MN occurs at 1.34 s of calculation time, which corresponds to a depth of penetration of
the structure of 0.67 m, and then decreases to average values around 6.23 MN. Obviously, this does not
correspond to the real mechanics of the ice-structure interaction. Also, these values are several times
higher than other cases. The picture of ice destruction at the end of the calculation is presented in Figure
3.4.
80
Figure 3.4 – Ice crushing during impact of ice field in case of structured hex mesh: a) general view;
b) top view; c) side view
As we can see, the vertical cohesive elements destroyed along several parallel faces of the central
finite elements. To understand the reason, it is necessary to consider the very beginning of the
interaction. Figure 3.5 shows that, at first, the two closest finite elements to the structure begin to interact
with it. Then, shear stresses 𝜏𝑥𝑧 arise on the lateral faces of these elements. Since the mesh is structured,
the propagation of stresses passes deep into the field in a straight line without any bulk elements in the
pass. The cause of the unrealistic destruction of vertical cohesive elements is the rapid propagation of
shear stresses 𝜏𝑥𝑧 , while cohesive elements orthogonal to them are not working. This is clearly seen
when looking at the fracture pattern presented in Figure 3.6.
Figure 3.5 – Propagation of shear stress 𝜏𝑥𝑧 in structured hex mesh: a) at time 0.01 s; b) at time 0.02 s
Figure 3.6 – Fracture pattern at time 5 s (structured hex mesh)
81
Some similarities can be seen with the splitting failure mode described in section 1 and presented
in Figure 1.25, but this is not the right case, since longitudinal separating cracks form on the faces of
each row of finite elements. Also, the horizontal load along y-axis shows very unstable behavior. Some
rapidly changing peaks can be seen from time 1 s. Thus, for the correct description of the destruction
process, the use of a structured hex mesh is not allowed.
The second type of mesh to consider is unstructured hex mesh. The result of the calculation at the
last second of the interaction is presented in Figure 3.7. As for the general picture of the interaction it is
obvious that the asymmetric destruction occurs with the breaking away of large parts of the ice field.
Although the presence of an unstructured mesh made it possible to avoid the propagation of shear
stresses over the entire width of the modeled field, there are still small regions consisting of “direct”
parallel faces of elements form weak zones (Figure 3.8) along which the field is destroyed and which
are the cause of the asymmetric interaction picture. Reducing the mesh size in the case of a hexagonal
mesh will only worsen the situation, since a decrease in the size of the elements leads to a straightening
of the boundaries of elements and the mesh tends to become structured.
Figure 3.7 – Ice crushing during impact of ice field in case of unstructured hex mesh: a) general view;
b) top view; c) side view
Figure 3.8 – Shear faults at time 0.58 s in case of unstructured hex mesh
82
The ice load in this case rapidly increases to high values (about 9.785 MN), and then gradually
decreases to values of about 1,5 MN. The maximum value corresponds to the beginning of interaction
(time range 0.3-0.6 s), that is obviously not correct. The effect of asymmetric fracture can also be seen
by looking at the load history in the y direction (Figure 3.3). After splitting off a large fragment of the
field, there was no contact with solid ice on this side. In this regard, the horizontal load from ice on the
other side of the structure along the y-axis was significant throughout almost the entire calculation time.
Fracture pattern is shown in Figure 3.9. Given all the above, we can conclude that despite the efficiency
of computing, the use of a hexagonal mesh at all is unacceptable when modeling the impact of ice on
structures.
Figure 3.9 – Fracture pattern at time 5 s (unstructured hex mesh)
Next, consider the prism mesh case. The result of the calculation at the last second of the interaction
is presented in the Figure 3.10. This case shows more realistic picture. Spalls and slight hummocking
are present. An interesting fact is the periodic formation of “wing” fracture surfaces in an ice field. For
example, as can be seen in Figure 3.11 (a), at the beginning of the interaction, the mean (hydrostatic)
stresses have increased values in the region of future cracks. With further interaction, the ice fails mainly
on these surfaces, then the process repeats Figure 3.11 (b). For clarity, the fact of the formation of such
cracks is shown in the Figure 3.12.
As for the ice load, then, as in the past case, it quickly increases at the beginning of the interaction
to the value 4.398 MN at 0.15 s, and then decreases to average values of about 0.318 MN. However,
further behavior shows a different result then previous one. From about 1.8 s, the load begins to slowly
increase and reaches its peak value of 7.358 kN at a time of 4.49 seconds, which corresponds to almost
full-diameter penetration of the structure into the field. At a time of about 4.7 s, the horizontal load along
83
the y axis increased to 1.2 MN due to the destruction of ice on one side of the structure, but then quickly
returned to low values.
The presence of a peak in ice load at the beginning of the interaction is connected with a small
contact area and in my opinion is not correct. As can be seen from Figure 3.11, a, one finite element has
a significantly higher level of stress than the other surrounding elements. This may be due to many
factors, such as improper contact stiffness in the model, or the incorrect principle of a mesh creation in
the contact zone. Therefore, additional studies of the applicability of this pattern are required.
Figure 3.10 – Ice crushing during impact of ice field in case of prism mesh: a) general view; b) top view;
c) side view
Figure 3.11 – Formation of “wing” failure surfaces in ice field in case of prism mesh pattern: a) at time
0.2 s; b) at time 5 s
The disadvantages of this type of finite elements can be seen by considering the central section of
the model (Figure 3.13). Although the destruction pattern is acceptable in top view, in this section there
is a complete discrepancy with real observations of ice impacts. With the current drift velocity and the
diameter of the structure, the type of interaction should correspond to crushing failure mode of the ice
field at the contact. Bekker's monograph [9] describes the process of field destruction in thickness for
84
this mode. In nature, in the upper and lower layers of ice spalls are formed along inclined surfaces,
contributing to the concentration of stresses in the central part of the ice field. However, the Figure
3.13, c, shows that after partial crushing of the ice, the horizontal connected finite elements are destroyed
simultaneously over the entire thickness, due to which the ice field is more or less equally vertically
stressed. With an increase in the number of finite elements in thickness, a different picture may be
observed, however, a decrease in the size of the finite elements may not be practical in comparison with
the use of tetrahedral finite elements.
Figure 3.12 – Fracture pattern at time 5 s (prism mesh)
Figure 3.13 – Central vertical section of model during interaction (prism mesh): a) at 0.01 s; b) at 0.1 s;
c) at 0.25 s; d) at 0.55 s
85
Modeling of ice with a prismatic mesh requires significantly longer calculation time than in the
case of a hexagonal mesh, but such a mesh allows to take into account the uneven development of cracks
in the vertical plane to a good extent.
The last case to consider is a tetrahedral mesh. This case is considered the most expensive in terms
of computation cost. This is due to the much larger number of finite elements needed to discretize the
model with the same size settings. The calculation was not completed until the end, since at a time point
of 4.01 s in the simulation, the time increment became very small, possibly due to strongly deformed
elements. As a result, to calculate a model with the same initial conditions using the tetra mesh, almost
6 times more time is required than in the case of prism mesh. The simulation result for time in model
equal to 4.01 s is shown in Figure 3.14. The deformation pattern in this case is quite different from the
previous ones. The fracture pattern is shown in Figure 3.15.
Figure 3.14 – Ice crushing during impact of ice field in case of tetrahedral mesh: a) general view; b) top
view; c) side view
Figure 3.15 – Fracture pattern at time 4.01 s (tetra mesh)
86
In this case, there are no large chipped pieces as well as cracks propagation deep into the ice field.
Figure 3.14, c shows a hummocking process, which has a much greater degree than in previous cases.
This picture clearly corresponds to the real case of destruction, when a gradual and uniform destruction
of ice at the contact is realized. The reason for the strong difference between the tetrahedral mesh and
the prismatic one can be seen if we analyze the beginning of the interaction and the central vertical
section of ice field.
As can be seen in Figure 3.16, a, at the very beginning (up to 0.1 s), the main destruction begins
along the horizontal and vertical faces of the elements and the difference with the prismatic mesh is not
noticeable. However, with further interaction Figure 3.16, c, the inclined faces of the tetrahedral
elements are included in the work and serve as fracture surfaces.
Figure 3.16 – The development of ice destruction at the contact: a) at 0.1 s; b) at 0.2 s; c) at 0.35 s; d)
at 0.5 s
The qualitative picture of fracture, as in the case of a prismatic mesh, can be estimated by
considering the interaction in the central section of the model. This section is presented in Figure 3.17.
Unlike the case of a prismatic mesh, the pattern of destruction is more consistent with the process
described by Bekker (2004) [9]. At the beginning of interaction (Figure 3.17, b) the central part of the
ice field has the highest stress values, since spalls have occurred in the upper and lower parts. Uneven
destruction of the ice field in thickness is clearly visible in the Figure 3.17, d. The lower part of the ice
field is destroyed to a greater extent due to high ice temperatures and, accordingly, low ice strength.
Thus, the fracture process in this case is periodic in nature with the formation of inclined fracture
surfaces. This interaction physics will be better visible if the number of finite elements in the thickness
of the ice field is increased.
87
Figure 3.17 – Central vertical section of model during interaction: a) at 0 s; b) at 0.01 s; c) at 0.05 s;
d) at 0.25 s
As for the nature of the ice load, this pattern is the only one considered that has almost no peak
load at the beginning of the interaction. The presence of this peak can also be due to the influence of
inclined faces and spalls, that are destroyed and reduce the contact area. Also, load is increasing in slow
manner during analysis. Unfortunately, it was impossible to estimate the load when embedding the
structure at full width, because the calculation was interrupted. An additional study of the criterion for
the removal of bulk and cohesive elements can solve the problem of a strong decrease in the time step.
At the end of the study of the mesh pattern, the following conclusions can be drawn:
- the hexagonal mesh showed the highest value from the calculated loads;
- the presence of a structured mesh speeds up the calculation, but leads to an incorrect picture of
the interaction and stress distribution in an ice field;
- the presence of inclined faces of finite elements significantly changes the pattern of destruction;
- the tetrahedral mesh has significantly higher computational cost;
- the tetrahedral mesh shows the most realistic picture of the interaction.
In my opinion, using a tetrahedral finite element would be the best solution, since this case has an
acceptable picture of failure and the nature of the load. It has high computational cost, but this problem
can be solved by using more powerful computers. As for the process of interaction with structures in the
first seconds, a rapidly increasing load can be a problem of ensuring the correct initial contact.
88
3.2 Element size influence
At this stage of the model study, the influence of the size of finite elements on the final simulation
result was considered. As the first case in a series of experiments, the model from the previous subsection
(tetrahedral mesh pattern) was adopted. When creating a new model, all sizes and boundary conditions
were preserved, but the size of the finite elements of the ice field was changed. Moreover, to optimize
the calculation time, the field was divided into 4 zones as presented in Figure 3.18.
Figure 3.18 – Mesh zones with different element sizes
The first zone has a radius equal to the diameter of the structure and the minimum size of the finite
elements is approximately 0.15 m. By thickness, the ice field is discretized into 8 layers. The second
zone – two diameters with an approximate size of the finite elements equal to 0.3 m (4 elements in
thickness). The third and fourth zones have 2 and 1 end elements in thickness, respectively, with a
maximum size of about 2.2 m. This decision was made based on the destruction pattern presented in the
Figure 3.15, namely, it is clear that zone of destroyed cohesive elements does not extend much deep into
the field, when the field interacts almost over the full diameter. Thus, a more accurate description of
contact with the structure was provided, and the total number of finite elements was increased by only
1.5 times. The calculation results are presented in the Table 3.4. The history of the load along and across
the direction of movement is presented in the Figure 3.19 and Figure 3.20.
The nature of the load in the cases considered is, as expected, very similar, but for the case of the
fine mesh, the vibrations are smoother. Moreover, with a decrease in the size of the elements, the load
has higher values throughout almost the entire interaction time. The picture of destruction in the plan is
shown in the Figure 3.21 and visually not much different. The central section of the model is shown in
Figure 3.22.
89
Table 3.4 – General modeling information
Case
0.30 m (4 elements in
thickness)
0.15 m (8 elements in
thickness)
Number of finite
elements (bulk /
cohesive)
Total calculation time,
hours
Peak total ice force
along direction of
motion 𝐹𝑥𝑡𝑜𝑡 , MN
124 261 / 236 075
48.63
(4.01 s out of 5)
-
195 904 / 378 476
141.68
10,096
Figure 3.19 – Load history with different mesh sizes (along x-axis)
Figure 3.20 – Load history with different mesh sizes (along y-axis)
90
Figure 3.21 – Fracture pattern at time 5 s (model with fine mesh)
Figure 3.22 – Central vertical section of model with 8 elements in thickness during interaction: a) at
0.01 s; b) at 0.10 s; c) at 0.25 s; d) at 0.55 s
Noticeable improvements can be seen in the overall picture of the interaction. As in the case of a
coarse mesh, at the beginning of the interaction, high compressive stresses are concentrated in the center
of contact zone (Figure 3.22, a). A small peak of the load that occurs in the first fractions of a second is
associated with the process of deformation and growth of stresses in bulk elements. Then the first
spalling occurs, and the total ice force drops to low values. A further increase in the load is associated
91
with an increase in the contact area and, accordingly, the number of contacting elements. The inclined
surfaces of destruction are now visible even more clearly (Figure 3.22, d).
Also in the process of destruction, there were moments when the destruction of the ice field did
not occur along two, but along one inclined surface throughout the thickness as presented in Figure 3.23.
Such a picture can also develop and is associated with low strength of the lower layers, which are
destroyed earlier than the rest.
Figure 3.23 – Central vertical section of model with 8 elements in thickness at time 4.0 s
3.3 Section conclusions
The results of modeling the ice impacts on vertical structures according to the accepted method
showed good results. As expected from the finite element method, the influence of various mesh patterns
on the picture of ice field destruction is large. The presence of inclined cohesive elements leads to a
more realistic picture of the interaction. Therefore, an ice field of tetrahedral bulk elements is the most
acceptable option from the considered. The main problems when using tetrahedral elements are as
follows:
- assignment of ice properties to inclined cohesive elements, since there is little information on the
deformation and strength characteristics of ice along inclined axes;
- big time for calculating models, which can be solved by using a more powerful computer.
The dependence of the load on the size of the finite elements was weak for the considered case.
However, the pattern of ice destruction at the contact is much better when using a fine mesh and
coincides with the real observations of some authors.
The nature of the destruction of ice in general corresponds to the crushing failure mode in real
conditions. Unequal destruction of the ice field in height during modeling is a significant factor that
92
confirms the relevance and reliability of the applied method. Also, this nature of the destruction is
contrary to Gladkov’s assumptions [27], on which the current methodology for determining ice loads in
the regulatory document SP 38.13330.2018 [68] is based. In determining the values of the factor 𝑘𝑏 , he
admitted that all layers of the ice field are destroyed simultaneously.
To assess the possibility of applying the considered methodology for solving real engineering
problems, it is necessary to verify it.
93
4 Verification of the methodology for modeling ice impacts on vertical offshore structures
4.1 Comparison of load with analytical methods of national codes
To confirm the efficiency of the model for solving real engineering problems, it is necessary to
verify the used methodology of ice-structure interaction modeling. This can be done, for example, by
methods such as analyzing the correspondence of the qualitative picture of ice destruction to recorded
cases and comparing the magnitude of the impact with the measured in the field.
Another way is to compare the simulation results with the analytically calculated value of the load
according to various theories. Although the analytical methods used to date have several disadvantages,
they are still successfully used in the design of offshore structures. The first stage of verification will be
the comparison of magnitude of the ice load obtained as a result of modeling with the load calculated by
the methods described in the national codes below. Most of national standards require to take into
account several limiting factors of ice load described earlier (section 1). The methods for calculating the
corresponding scenarios of the interaction of ice formations and structures have large differences.
However, since the principle of infinite ice field is applied in the modeling, a comparison with the load
only for limit stress scenario will be performed.
SP 38.13330.2018 [68]. This regulatory document is a standard of the Russian Federation and
contains requirements for the calculation, design, justification of the reliability and safety of hydraulic
structures. In accordance with these standard, the ice load form drifting ice field on structures with a
vertical front face in case of limit stress scenario 𝐹𝑏,𝑝 , MN, is determined by the modified Korzhavin’s
formula proposed by Gladkov M. G. [27]:
𝐹𝑏,𝑝 = 𝑚𝑘𝑏 𝑘𝑉 𝑅𝑏ℎ,
(4.1)
where 𝑚 – structure shape factor;
𝑘𝑏 – indentation factor taking into account the influence of complex stress-strain state of ice field
on the ice strength limit;
𝑘𝑣 – factor taking into account the influence strain rate on the ice strength limit;
𝑅 – compressive strength of ice, MPa;
𝑏 – structure width, m;
ℎ - ice thickness, m.
STO Gazprom 2-3.7-29-2005 [69]. This document is a standard of the Russian company
Gazprom and contains methods for calculating the ice load, supplementing and clarifying the design
requirements, and the provisions of Russian national codes. The load from the impact of ice formations
is determined with the values of 𝑑/ℎ ≥ 10 according to following equation
𝐹𝑏,𝑝 = 𝑚𝑘𝑏 𝑅𝑏ℎ,
94
(4.2)
The parameters in equation (4.2) are the same as in equation (4.1).
ND 2-020201-015 [52]. These rules of the Russian Maritime Register of Shipping establish
requirements that are specific to floating drilling rigs and offshore stationary platforms, take into account
the recommendations of the IMO Code for the design and equipment of floating drilling rigs adopted by
the IMO Assembly on December 2, 2009. In this case ice load is determined by equation
𝐹𝑏,𝑝 = 𝑚𝑘𝐿 𝑘𝜐 𝑅𝑏 0.85 ℎ0.9 ,
(4.3)
where 𝑘𝐿 – factor that takes into account the effect of ice field area and structure diameter on the
load.
Other parameters in equation (4.3) are the same as in equation (4.1).
ISO/FDIS 19906:2019(E) [33]. This document was developed by the International Organization
for Standardization (ISO) and establishes requirements, provides recommendations for the design,
construction, transportation, installation and decommissioning of offshore structures associated with the
activities of the oil and gas industry in the Arctic and cold regions. Based on this document the global
ice load is determined as follows:
𝐹𝑏,𝑝
ℎ 𝑛 𝑏 𝑚
= 𝑏ℎ𝐶𝑅 (( ) ( ) + 𝑓𝐴𝑅 ) ,
ℎ1
ℎ
(4.4)
where 𝐶𝑅 – the ice strength coefficient, MPa;
𝑛 and 𝑚 – empirical coefficients;
𝑓𝐴𝑅 – empirical term given by
ℎ
𝑓𝐴𝑅 = 𝑒 −𝑏⁄3ℎ √1 + 5 ,
𝑏
(4.5)
CAN/CSA-S471-04 [11]. This document is Canada's national standard.
This standard use power-law dependencies of ice pressure on nominal contact area, which were
obtained as a result of processing a large amount of ice pressure data on the hulls of the Kigoriak, Polar
Sea, MV Arctic, Manhattan and Oden icebreakers in Canadian waters. The code recommends two
different formulas for calculating the ice load depending on the ratio of the width of the structure to the
thickness of the ice field.
If the global ice load 𝐹𝑏,𝑝 is determined as follows:
(𝐷 +1)
𝐶𝑝 𝐴𝑁 𝑝
, 𝑏/ℎ < 10
𝐹𝑏,𝑝 = {
𝐶𝑝 ℎ(𝐷𝑝−𝐸𝑝 +1) ∙ 𝑏 (𝐸𝑝 +1) , 𝑏/ℎ ≥ 10
(4.6)
where 𝐴𝑁 – nominal contact area, m2;
𝐶𝑝 , 𝐷𝑝 и 𝐸𝑝 – empirical parameters.
Elforsk rapport 09:55 [24]. This report by the Swedish energy research company contains
recommendations, a description of the ice load mechanisms and statistics from the Baltic Sea. In
95
accordance with the recommendations of this report, the ice load from the level ice fields acting on
structures can be evaluated by the original Korzhavin’s equation:
𝐹𝑏,𝑝 = 𝑚𝑘𝑐 𝑘𝑏 𝑅𝑏ℎ,
(4.7)
where 𝑘𝑐 – contact factor that covers the fact that ice under continuous crushing is not in contact with
the whole nominal area.
GL 2005. IV-Part 6 [26]. These rules are a German national standard and are designed to evaluate
the properties of sea ice and ice loads on offshore structures. Ice load is determined using the formula
obtained by several authors at Iowa University based on small-scale tests:
𝐹𝑏,𝑝 = 𝑘𝑏 0.5 ℎ1.1 𝑅,
(4.8)
where 𝑘 – empirical factor.
As a simulated case of ice impact on structures, the case of a tetrahedral mesh, described in
section 3.1 of this thesis, is accepted. The maximum load was considered in two ways: load surges were
considered and not considered. This decision was made because the cause of these surges is unknown
and their number is very small compared to other data points. Calculations using analytical methods
were carried out strictly in accordance with all clauses of national standards. The strength of ice was
calculated in accordance with the recommendations of each regulatory document. In the case when there
is no indication in the codes on the strength of ice, then its value was taken equal to that calculated by
the SP 38.13330.2018. A comparison of the results is presented in the Table 4.1 and in the Figure 4.1.
Table 4.1 – Comparison of ice load and difference with simulation results
Parameter
Modeling
(with
surges)
SP
Ice load, MN
10.10
Difference
with modeled
load (surges
are
considered), %
Difference
with modeled
load (surges
are not
considered), %
MRS
ISO
FDIS
CAN
CSA
Elforsk
GL
Modeling
(without
surges)
7.83
11.44
9.36
3.94
1.59
4.17
4.59
0.00
28.99
13.29
7.85
156.25
534.55
142.10
119.72
119.72
70.34
148.92
103.72
16.63
188.80
10.19
0.00
The magnitude of the ice load calculated by various analytical methods in some cases differs
several times. If we compare these values with the values of the load obtained as a result of modeling
we get a mixed result.
96
Figure 4.1 – Comparison of the ice load according to different codes with the result of numerical
simulation
If we take into account surges, then the load has value that is close to calculated according to more
modern standards (SP 38.13330.2018, MRS 2018 №2-020201-015, ISO/FDIS 19906:2019(E)). In this
case the methodology of the ISO normative document showed the greatest convergence, but when
compared with other standards, the difference reaches 534.55% for Elforsk rapport 09:55. This fact may
also be a coincidence, therefore, it is impossible to judge the correctness of the ice load obtained during
modeling. Analytical methods also have their own assumptions and disadvantages, so it is necessary to
perform several iterations of calculations with different conditions, and then perform a comparison.
4.2 Modeling of uniaxial compression tests
Finally, the last approach to verify the numerical model is to simulate uniaxial compression tests.
Such laboratory data are valuable, since the parameters of ice and load can be regulated during the
experiment, and different cases can be considered. The author took part in a series of such experiments
conducted in Laboratory of Ice Mechanics of Far Eastern Federal University (FEFU). One of the main
goals of tests was to determine the compressive strength of sea ice at different strain rates.
Ice samples for testing were prepared as follows. First of all, ice cores were cut from ice in Novik
Bay, located near the university using the core sampler (Figure 4.2, a). The ice thickness at the time of
core sampling was approximately 25 cm. Then the cores were transported to the laboratory and placed
in the refrigerator. Second, samples of 200 mm height were cut out from cores, their sizes and weights
were measured. For the compression tests, 20 cylindrical samples with a diameter of 90 mm and a height
97
of 200 mm were prepared (Figure 4.2, b). The tests were performed using the Shimadzu press (Figure
4.2, c). The maximum possible load on the press is 100 kN. Apparatus has the possibility to control
loading speed. Test were performed on following speeds: 0.1 mm/s, 0.5 mm/s, 1 mm/s, 1.5 mm/s, 2.5
mm/s using 5 samples in each case. After crushing the temperature of each sample was measured.
Figure 4.2 – Compression tests preparation: a) core sampling; b) sample setting; c) Shimadzu press
During the verification process, the reaction force magnitude of the upper plate, as well as the
overall pattern of fracture of the samples will be compared with experimental data. The finite element
model with boundary conditions is presented in Figure 4.3.
When creating a numerical model of experiments, the following assumptions were made:
- support and loading plates of the press are modeled using rigid shells;
- roughnesses in the surface of the samples are not taken into account;
- the size of the finite elements of the sample is taken equal to 12 mm;
- the coefficient of friction is taken equal to 0 to prevent the confining of displacements in the
vicinity of plates.
Figure 4.3 – Finite element model in numerical experiments on uniaxial compression of ice samples
A comparison of forces in modeling and laboratory experiment is presented in the Figure 4.4. The
results were mixed. On the one hand, the nature of the load is very different. for example, during
98
modeling, there are constant power surges throughout the entire time. The output frequency was assumed
to be 0.01 second, as in a laboratory experiment. On the other hand, the force graph obtained from the
results of a laboratory experiment runs right through the upper points almost until the sample is
destroyed. In my opinion, in general, the result is acceptable, because when evaluating the force along
the upper boundary, as when comparing the ice load according to regulatory documents, there is a good
agreement.
Figure 4.4 – Comparison of forces in modeling and laboratory experiment
4.3 Section conclusions
In this section, verification of the developed methodology for numerical modeling of the impact
of ice formations on the vertical supports of structures was performed.
This modeling methodology showed a good correspondence of the load values to the load
determined by the current documents. Even among themselves, the methods of various codes are very
different from each other. For example, the very important fact that almost all foreign standards consider
the regime of viscous-brittle fracture when the strain rate is within approximately 5∙10-4÷10-3 s-1. Thus,
such regulatory documents assume the maximum strength of ice during interaction. The only exceptions
is the code GL 2005. IV-Part 6 [26] who propose using the empirical formula obtained by Kovacs (1997)
[39] to calculate the strength of ice based on tests of ice samples from the Beaufort Sea. The Russian
standards have a 𝑘𝜐 factor which was first introduced in SNiP 2.06.04-82* at the proposal of M. G.
Gladkov [27], although Korzhavin, when developing the original formula, also took into account the
influence of the strain rate on the strength of ice. The proposed method takes into account the rate of ice
deformation during the entire simulation time, which is an undoubted advantage over analytical methods.
99
As for verification using laboratory test data, the result is mixed. On the one hand, the general
nature of the increase in load during modeling correlates well with the laboratory experiment. Also, the
maximum force does not differ much from that obtained in the experiment. On the other hand, when
modeling, there are multiple fluctuations in the force, although the frequency of the data corresponds to
the frequency of recording the force of the press machine. This can be due to many factors and is one of
the problems that can be solved with further research.
100
Conclusion
Summarizing the results obtained in this dissertation, the following conclusions can be drawn:
1. Based on the analysis of the current state of problem of assessing the impact of ice fields on
structures, it can be concluded that this problem has been solved for a rather long time. In terms of
determining the physicomechanical characteristics of ice, a significant amount of research has been
performed. Numerical modeling is a promising way to solve this problem. To date, there are a number
of works in this area. In almost all the reviewed articles and theses, the authors did not pay due attention
to the selection and substantiation of constitutive models of ice material. Such features of ice as the
strong dependence of the parameters on temperature, the complex stress-strain state, and also the strain
rate were not taken into account.
2. A comprehensive study of methods for determining the parameters of ice, depending on many
factors was made. An algorithm has been compiled to calculate the full set of necessary characteristics
of ice, depending on some initial conditions.
3. As a method for solving the problem, a numerical finite element method is chosen. The solution
tool is the SIMULIA Abaqus software package. As an integration method for the main equation of
motion, an explicit integration method was justified and adopted. Tasks were set and a plan was drawn
up for numerical experiments performed in the thesis. Initial numerical models of the interactions of the
ice field with structures were created on the basis of some assumptions based on the known results of
experimental studies of sea ice.
4. The statement of the problem based on the developed methodology of numerical modeling
includes the following studies and improvements:
- the tetrahedral mesh pattern is accepted as the most optimal from the point of view of a qualitative
picture of interaction;
- reducing the size of the finite elements does not lead to a significant change in the magnitude of
the ice load, but contributes to a more detailed description of the process of ice loads formation and the
formation of a more accurate picture of ice destruction;
- the accepted threshold in the criterion for the removal of bulk elements was found to be applicable
5. The numerical model of interaction with the introduced improvements was verified in two ways:
сomparison of load with analytical methods of national codes and modeling of uniaxial compression
tests. As a result of verification calculations, numerical models showed good convergence with the
reference data. The nature of the destruction of the ice field during modeling is in good agreement with
field observation, for example, the formation of piles of ice fragments (hummocking) and uneven
destruction of the ice field in thickness.
6. The developed methodology for numerical modeling of ice formations impacts can be used to
101
solve real engineering problems and will significantly improve the reliability and reduce the risk of
offshore structures. The method will allow to evaluate the complex stress-strain state of structure under
action of ice loads from drifting ice fields.
7. One of the main problems in using the developed modeling approach is the high computational
cost. When solving real engineering problems, it is recommended to use powerful computers with a lot
of CPU cores.
8. The model can be more deeply researched and improved in the following terms:
- improving methods for calculating the deformation and strength parameters of ice;
- taking into account the influence of deformations in mutually perpendicular directions in the
material model of cohesive elements, i.e., a coupled traction model can be used;
- to study the influence of the general damping parameters of the model (linear and quadratic
viscosity);
- searching for appropriate initial contact between finite elements of structure and ice field;
- assessment of the degree of influence of various water models on simulation result
- application of the Coupled-Eulerian-Lagrangian method in numerical model;
- verification of the model can be supplemented by modeling real recorded cases of ice impacts on
structures.
102
References
1 Albrektsen A. Prediction of the response from ice forces to a lighthouse structure: master's
thesis. - Trondheim : Department of Structural Engineering, NTNU, 2009.
2 Aleksandrov A.V., Platonov V.V., Matantsev R.A. Study of ice failure and deformation
processes [Conference] // Proceedings of the ASME 2014 33rd International Conference on Ocean,
Offshore and Arctic Engineering OMAE2014. - San Francisco : American Society of Mechanical
Engineers, 2014.
3 Alexandrov A.V., Matantsev R.A. Primeneniye metodov komp'yuternogo modelirovaniya
v reshenii zadach masshtabirovaniya ledovykh nagruzok [Application of computer simulation methods
in solving problems of scaling ice loads] [Journal] // Trudy Krylovskogo gosudarstvennogo tsentra Proceedings of the Krylov State Center. - St. Petersburg : FSUE Krylov State Scientific Center, 2017. 379 : Vol. I. - pp. 70-77.
4 ANSYS Documentation. Mechanical APDL. Material reference. Chapter 4: Nonlinear
material properties.
5 ANSYS Documentation. Mechanical APDL. Theory reference. Chapter 4: Structures with
Material Nonlinearities.
6 ANSYS Documentation. Mechanical Applications. Explicit Dynamics Analysis Guide.
Chapter 3: Transforming an Implicit Model to run in Explicit Dynamics.
7 ANSYS Documentation. Mechanical Applications. Explicit Dynamics Analysis Guide.
Chapter 6: Explicit Dynamics Theory Guide.
8 Assur A. Flexural and Other Properties of Sea Ice Sheets: Condensed Version [Conference] //
Physics of Snow and Ice : proceedings. - 1967. - Vol. 1. - pp. 557-567.
9 Bekker A.T. Veroyatnostnyye kharakteristiki ledovykh nagruzok na sooruzheniya
kontinental'nogo shel'fa [Probabilistic characteristics of ice loads on the continental shelf structures]
[Monograph]. - Vladivostok : Dalnauka, 2004.
10 Bjerkås M. Meese A., Alsos H. Ice Induced Vibrations – Observations of a Full Scale Lockin Event [Conference] // Proceedings of the 23rd International Offshore and Polar Engineering. Anchorage : International Society of Offshore and Polar Engineers (ISOPE), 2013. - pp. 1272-1279.
11 CAN/CSA S471-04. General requirements, design criteria, the environment and loads. Ontario : Canadion standards association, 2004. - 116 p.
12 Carney K., Benson D., DuBois P., Lee R. A phenomenological high strain rate model with
failure for ice [Journal] // International Journal of Solids and Structures. - 2006. - 43. - pp. 7820-7839.
13 Cole D. M. Strain-rate and grain-size effects in ice [Journal] // Journal of Glaciology. [s.l.] : Cambridge University Press, 1987. - 115 : Vol. 33. - pp. 274-280.
103
14 Cox G. F. N., Weeks W. F. Equations for determining the gas and brine volumes in sea ice
samples [Journal] // Journal of Glaciology. - 1983. - 102 : Vol. 29. - pp. 306-316
15 Cox G. F. N., Weeks W. F. Salinity variation in sea ice [Journal] // Journal of Glaciology. 1974. - 67 : Vol. 23. - pp. 109-120.
16 Cundall P. A. Strack O. D. L. The distinct element method as a tool for research in granular
media [Book]. - Minneapolis : Department of Civil and Mineral Engineering, University of Minnesota,
1978. - Vol. I.
17 Daiyan H., Sand B. Numerical Simulation of the Ice-Structure Interaction in LS-DYNA
[Conference] // Proceedings of 8th European LS-DYNA Users Conference. - Strasbourg : [s.n.], 2011.
18 Derradji-Aouat Ahmed Explicit FEA and constitutive modelling of damage and fracture
on polycrystalline ice - simulations of ice loads on offshore structures [Conference] // Proc., 18th Int.
Conf. on Port and Ocean Eng. Under Arctic Conditions. - Canada : , 2005. - Vol. I. - pp. 225-238.
19 Doronin Yu. P. Fizika okeana [Ocean Physics] [Book]. - St. Petersburg : Russian State
Hydrometeorological University, 2000. - 305 p. (in Russian).
20 Doronin Yu. P., Kheisin D. E. Morskoy lod [Sea ice] [Book]. - Leningrad :
Gidrometeoizdat, 1975. - 320 p. (in Russian).
21 Dykins J. E. Tensile Properties of sea ice grown in a confined system [Conference] //
Physics of Snow and Ice: proceedings. - Japan : Bunyeido Printing, 1967. - Vol. 1. - pp. 523-537.
22 Feng D., Pang D. S., Zhang J. Parameter Sensitivity in Numerical Modelling of IceStructure Interaction With Cohesive Element Method [Conference] // Proceedings of the ASME 2016
35th International Conference on Ocean, Offshore and Arctic Engineering (OMAE). - Busan : ASME,
2016.
23 Frankenstein G., Garner R. Equations for determining the brine volume of sea ice from 0.5 °C to -22.9 °C [Journal] // Journal of Glaciology. - 1967. - 48 : Vol. 6. - pp. 943-944.
24 Fransson L. Bergdahl L Elforsk rapport 09:55 // Recommendations for design of offshore
foundations exposed to ice loads. - 2009. - 43 p.
25 Gammon P. H., Kiefte H., Clouter M. J., Denner W. W. Elastic constants of artificial and
natural ice samples by Brillouin spectroscopy [Journal] // Journal of Glaciology. - [s.l.] : Cambridge
University Press, 1983. - 103 : Vol. 29. - pp. 433-460
26 GL 2005 IV. Part 6 - Offshore Installations. Chapter 7 - Guideline for the Construction of
Fixed Offshore Installations. - Hamburg : Germanischer Lloyd, 2005. - 29 p.
27 Gladkov M. G. Nagruzki i vozdeystviya l'da na morskiye gidrotekhnicheskiye
sooruzheniya: dis. … dokt. tekh. nauk [Loads and impacts of ice on marine hydraulic structures. Dr.
thesis (Eng.)]. - St. Petersburg : [s.n.], 1997. - 181 p.
28 Gold L. W. Engineering properties of fresh-water ice [Journal] // Journal of Glaciology. -
104
[s.l.] : Cambridge University Press, 1977. - 81 : Vol. 19. - pp. 197-212.
29 Gürtner A. Bjerkås M., Forsberg J., Hilding D. Numerical modelling of a full scale ice
event [Conference] // Proceedings of 20th IAHR International Symposium on Ice. - Lahti : [s.n.], 2010.
30 Gürtner A., Bjerkås M., Kühnlein W., Jochmann P., Konuk I. Numerical Simulation of
Ice Action to a Lighthouse [Conference] // Proceedings of the ASME 28th International Conference on
Ocean, Offshore and Arctic Engineering (OMAE). - Honolulu : ASME, 2009.
31 Hilding D., Forsberg J., Gürtner A. Simulation of ice action loads on offshore structures
[Conference] // Proceedings of 8th European LS-DYNA Users Conference. - Strasbourg : DYNAmore
GmbH, 2011.
32 Hilding D., Forsberg J., Gürtner A. Simulation of Loads from Drifting Ice Sheets on
Offshore Structures [Conference] // Proceedings of 12th International LS-DYNA Users Conference. Detroit : Livermore Software Technology Corporation, 2012.
33 ISO/FDIS 19906:2019(E). Petroleum and natural gas industries — Arctic offshore
structures. Final Draft. International Organisation for Standartization. ISO/TC 67/SC 7. - Vernier : ISO
copyright office, 2019. - 583 p.
34 Jones S. J. The confined compressive strength of polycrystalline ice [Journal] // Journal of
Glaciology. - [s.l.] : Cambridge University Press, 1998. - 98 : Vol. 28. - pp. 171-177.
35 Jordaan I. J. Mechanics of ice-structure interaction [Journal] // Engineering Fracture
Mechanics. - [s.l.] : Elsevier, 2001. - 17 : Vol. 68. - pp. 1923-1960.
36 Karulina M., Shkhinek K., Thomas G Theoretical and experimental investigations of level
ice interaction with fourlegged structures [Conference] // Proc. of 21st Int. Conf. on Port and Ocean Eng.
under Arctic cond., POAC 11. - Montreal : Curran Associates, Inc, 2011. - Vol. 1. - pp. 235-246. - ISBN:
978-1-61839-225-1.
37 Kim J-H., Kim Y. Numerical simulation of concrete abrasion induced by unbreakable ice
floes [Journal] // International Journal of Naval Architecture and Ocean Engineering. - [s.l.] : Elsevier
B.V., 2019. - 1 : Vol. 11. - pp. 59-69.
38 Korzhavin N. K. Vozdeystviye l'da na inzhenernyye sooruzheniya [The effect of ice on
engineering structures] [Книга]. - Novosibirsk : Publishing House of the Siberian Branch of the
Academy of Sciences of the USSR, 1962. - стр. 204 p. (in Russian).
39 Kovacs A. Estimating the full-scale flexural and compressive strength of first-year sea ice
[Journal] // Journal of geophysical research. - 1997. - С4 : Vol. 102. - pp. 8681-8689.
40 Kovacs A Sea Ice. Part I: Bulk Salinity Versus Ice Floe Thickness // CRREL Report 96-7. USA : CREEL, 1996.
41 Kovalev S. M., Lebedev G. A., Nedoshivin O. A., Sukhorukov K. K. Mekhanicheskiye
svoystvo morskogo l'da [Mechanical property of sea ice] [Book]. - St. Petersburg : Gidrometeoizdat,
105
2001. - p. 75 p. (in Russian).
42 Langleben M. P., Pounder, E. R. Elastic parameters of sea ice / ed. Kingery W. D.. - USA :
MIT Press, 1963. - pp. 69-78.
43 Li H. Bjerkås M., Høyland K. V., Nord T. S. Panel loads and weather conditions at
Norströmsgrund lighthouse 2000-2003 [Conference] // Proceedings of the 23rd IAHR International
Symposium on Ice. - Ann Arbor : [s.n.], 2016.
44 Li Liang Ledovaya nagruzka na gidrotekhnicheskiye sooruzheniya s naklonnoy gran'yu:
dis. … kand. tekh. nauk [Ice load on hydraulic structures with an inclined face: Cand. dis. (Eng.)]. - St.
Petersburg : [s.n.], 2014. - 136 p. (in Russian).
45 Liu J., Liu X., Chen Y. Numerical Simulations of Ice Loads on Fixed and Floating Offshore
Structures using the Discrete Element Method [Conference] // Proceedings of Arctic Technology
Conference. - St. John’s : Offshore Technology Conference, 2016.
46 Lu W., Lubbad R., Løset S., Høyland K. Cohesive Zone Method Based Simulations of
Ice Wedge Bending: a Comparative Study of Element Erosion, CEM, DEM and XFEM [Conference] //
Proceedings of 21st IAHR International Symposium on Ice. - Dalian : Dalian University of Technology
Press, 2012. - pp. 920-938.
47 Lu W., Heyn H.-M., Lubbad R., Løset S. A large scale simulation of floe-ice fractures and
validation against full-scale scenario [Journal] // International Journal of Naval Architecture and Ocean
Engineering. - [s.l.] : Elsevier B.V., 2018. - 3 : Vol. 10. - pp. 393-402.
48 Makarov O., Bekker A., Gogoladze D. Analysis of constitutive plasticity models in
relation to numerical modelling of ice impacts [Journal] // FEFU: School of Engineering Bulletin. Vladivostok : FEFU, 2020. - 43 : Vol. 2. - pp. 141-154.
49 Murat J. R., Lainey, L. M. Some experimental observations on the Poisson's ratio of seaice [Journal] // Cold Regions Science and Technology. - 1982. - 2 : Vol. 6. - pp. 105-113.
50 Nazintsev Yu. L., Dmitrash Zh. A., Moiseev V. I. Teplofizicheskiye svoystva morskogo
l'da [Thermophysical properties of sea ice] [Book]. - Leningrad : Publishing house of Leningrad State
University, 1988. - 260 p. (in Russian).
51 Nazintsev Yu. L., Panov V. V. Fazovyy sostav i teplofizicheskiye kharakteristiki morskogo
l'da [Phase composition and thermophysical characteristics of sea ice] [Book]. - St. Petersburg :
Gidrometeoizdat, 2000. - 84 p. (in Russian).
52 ND 2-020201-015 Rules for the classification, construction and equipment of mobile
offshore drilling units and fixed offshore platforms. - St. Petersburg : Russian Maritime Register of
Shipping, 2018 г.. - 454 p. (in Russian).
53 Nord T. S., Kvåle K. A., Petersen Øyvind W., Bjerkås M., Lourens E.-M. Operational
modal analysis on a lighthouse structure subjected to ice actions [Journal] // Procedia Engineering. -
106
2017. - Vol. 199. - pp. 1014-1019.
54 Palmer A. Croasdale K. R. Arctic Offshore Engineering [Book]. - Singapore : World
Scientific Publishing Co. Pte. Ltd., 2012. - p. 372.
55 Pernas-Sánchez J., Pedroche D. A., Varas D., López-Puente J. Numerical modeling of
ice behavior under high velocity impacts [Journal] // International Journal of Solids and Structures. 2012. - 14 : Vol. 49. - pp. 1919–1927.
56 Politko V. A.
Vozdeystviye
ledovykh
poley
na
odno-
i
mnogoopornyye
gidrotekhnicheskiye sooruzheniya s vertikal'noy opornoy chast'yu: dis. … kand. tekh. nauk [Influence
of ice fields on single and multi-bearing hydraulic structures with a vertical supporting part. Cand. dis.
(Eng.)]. - Moscow : [s.n.], 2018. - 148 p. (in Russian).
57 Richard M., McKenna R. Factors influencing managed sea ice loads [Conference] //
Proceedings of the 22nd International Conference on Port and Ocean Engineering under Arctic
Conditions (POAC). - Espoo : [s.n.], 2013.
58 Richter-Menge J. A., Jones K. F. The tensile strength of first-year sea ice [Journal] //
Journal of Glaciology. - [s.l.] : Cambridge University Press, 1993. - 133 : Vol. 39. - pp. 609-618.
59 Rist M. A. Murrel S. A. F. Ice triaxial deforntation and fracture [Journal] // Journal of
Glaciology. - [s.l.] : Cambridge University Press, 1994. - 135 : Vol. 40. - pp. 305-318.
60 Rüdiger von Bock und Polach, Sören Ehlers Model scale ice — Part B: Numerical model
[Journal] // Cold Regions Science and Technology. - [s.l.] : Elsevier, 2013. - 94. - pp. 53-60.
61 Ryvlin A. Ya. Metod prognozirovaniya predela prochnosti ledyanogo pokrova na izgib
[Method for predicting the ultimate strength of ice cover in bending] [Journal] // Problemy Arktiki i
Antarktiki - Problems of the Arctic and Antarctic. - 1974. - 45. - pp. 79-86. - (in Russian).
62 Salganik E. A., Shkhinek K. N. Ice induced vibrations of offshore structures [Journal] //
Magazine of Civil Engineering. - St. Petersburg : Peter the Great St. Petersburg Polytechnic University. 4 : Vol. 48. - pp. 72-78.
63 Sand Bjørnar Nonlinear finite element simulations of ice forces on offshore structures:
doctoral thesis. – Luleå., 2008. - 264 p.
64 Sanderson T. J. O. Ice mechanics : risks to offshore structures [Book]. - London : Graham
& Trotman, 1988. - 253 p.
65 Schulson E. M. Brittle Failure of Ice [Journal] // Engineering Fracture Mechanics. - 2001. 17-18 : Vol. 68. - pp. 1839-1887.
66 Schulson E. M., Duval P. Creep and fracture of ice [Book]. - New York : Cambridge
University Press, 2009. - p. 417.
67 Schwarz J. Weeks W. F. Engineering properties of sea ice [Journal] // Journal of
Glaciology. - [s.l.] : Cambridge University Press, 1977. - 81 : Vol. 19. - pp. 499-531.
107
68 SP 38.13330.2018 Loads and impacts on hydraulic structures (from wave, ice and ships). Moscow : Ministry of Construction of Russia, 2018. - 127 p. (in Russian).
69 STO Gazprom 2-3.7-29-2005 Metodika rascheta ledovykh nagruzok na ledostoykuyu
statsionarnuyu platformu [Methodology for calculating ice loads on an ice-resistant stationary
platform]. - [s.l.] : Gazprom LLC, 2005. - 16 p. (in Russian).
70 Taylor R., Frederking R., Jordaan I. The nature of high pressure zones in Compressive
Ice Failure [Conference] // Proceedings of The 19th IAHR International Symposium on Ice. Vancouver : St. Joseph Communications, 2008. - Vol. 2. - pp. 1001-1010. - ISBN: 978-0-9810446-0-6.
71 The Abaqus documentation. Abaqus. Elements. Special-Purpose Elements. Cohesive
elements. Defining the constitutive response of cohesive elements using a traction-separation
description.
72 The Abaqus documentation. Abaqus. Materials. Inelastic Mechanical Properties. Other
plasticity models. Extended Drucker-Prager models.
73 Thomas D. N. Dieckmann G. S. Sea Ice. Second Edition [Book]. - Chichester : Blackwell
Publishing Ltd, 2010. - p. 640.
74 Timco G. W., Weeks W. F. A review of the engineering properties of sea ice [Journal] //
Cold Regions Science and Technology. - : Elsevier B.V, 2010. - 2 : Vol. 60. - pp. 107-129.
75 Timko G. W., Frederking R. M. W. Compressive Strength o f Sea Ice Sheets [Journal] //
Cold Regions Science and Technology. - 1990. - 3 : Vol. 17. - pp. 227-240.
76 Timko G. W., Frederking R. M. W. Field Measurements of the Shear Strength of
Columnar-Grained Sea Ice [Conference] // Proceedings of 8th IAHR International Symposium on Ice. Iowa City : [s.n.], 1986. - Vol. 1. - pp. 279-292.
77 Tyshko K. P., Cherepanov N. V., Fedotov V. I. Kristallicheskoye stroyeniye morskogo
ledyanogo pokrova [Crystal structure of sea ice cover] [Book]. - St. Petersburg : Gidrometeoizdat,
2000. - 66 p. (in Russian).
78 Vaudrey K. D. Ice engineering: Study of related properties of floating sea-ice sheets and //
US Navy Civil Engineering Laboratory, Technical Report R 860. - 1977. - p. 81.
79 Wang Y., Zou Z.-J., Wang F., Lu T.-C. A Simulation Study on the Ice Fracture Behaviors
in Ice-Lighthouse Interaction Considering Initial Defects and Change of Elastic Modulus [Conference] //
Proceedings of the 29th International Ocean and Polar Engineering Conference. - Honolulu :
International Society of Offshore and Polar Engineers (ISOPE), 2019. - Vol. 173. - pp. 433-449.
80 Weeks W. F. Ackley S. F. The growth, structure and properties of sea ice // CRREL
Monograph 82-1. - Hanover : [s.n.], 1982. - p. 136.
81 Weeks W. F., Assur A. The mechanical properties of sea ice // U.S. Army CRREL
Monograph. - Hanover : U.S. Army Materiel Command, 1967. - Vol. 2.
108
82 Weiss J. Schulson E. M. The failure of fresh-water granular ice under multiaxial
compressive loading [Conference] // Proceedings of the 12th International Symposium on Ice of IAHR. Trondheim : [s.n.], 1994. - Vol. 1. - pp. 495-504.
83 Weiss J., Meyssonnier, J. Micro Mechanics of Ice // LOLEIF-Report No. 7 EU-Contract
No. MAS3-CT-97-0098. - 2001.
109
Отзывы:
Авторизуйтесь, чтобы оставить отзыв