CC BY-NC-ND
0ELSEVIER
Contents lists available at ScienceDirect
Cold Regions Science and Technology
journal homepage: www.elsevier.com/locate/coldregions
Numerical analysis of the bending strength of model-scale ice
Rüdiger U. Franz von Bock und Polach
Aalto University, Finland
Norwegian University of Science and Technology, Norway
ARTICLE INFO ABSTRACT
Performance simulation tools are of high significance for the design and especially the optimization of ships and offshore structures. However, for ice covered waters such tools are hardly available and are either costly as ice model tests or have a limited range of validity, such as semi-empirical formulas. This arises from the complexity of ice as material and insufficient knowledge on its mechanics. This paper presents a numerical analysis for model-scale ice in which material parameters are developed that can represent: tension, compression and in-situ downward bending. Those parameters are incorporated into a material model following the Lemaitre damage law. The developed material characteristics for model-scale ice are intended to support the design process of ships and offshore structures. The key phenomenon joining the deformation processes in bending together with those in compression and tension, proved to be the through thickness dependency of properties. This analysis and development is a continuation of previously presented parameters for compression and tension and is developed in agreement with experimental evidence.
© 2015 The Author. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
CrossMark
Article history: Received 27 October 2014 Received in revised form 4 June 2015 Accepted 9 June 2015 Available online 25 June 2015
Keywords: Model-scale ice Bending strength simulation Temperature gradient Strength property distribution
1. Introduction
The simulation of ice-structure and ice-ship interactions became in the past decade of increasing significance. This development is driven by increasing activities in the high North and increasing computational capacities. Ice model tests are still the state of the art method to validate the performance of designs, but time consuming and of significant costs. Therefore ice model tests are not suitable for design optimization, which however needs to be performed to determine economically suitable designs (see von Bock und Polach et al., 2014). Therefore, this paper is a contribution towards the development of a numerical version of the Aalto model ice basin.
1.1. The significance of ice failure in bending
The bending failure of ice is of high significance for ships in ice due to the inclined contact interfaces with the ice (see e.g. Enkvist, 1972; Lindqvist, 1989 or Valanto, 2001). Fig. 1. illustrates the bending failure of level ice after contact with the Finnish icebreaker Urho. This is also the reason, why the so called bending strength is often found in semi-empirical resistance formulations such as those of Lindqvist (1989). Fig. 2 illustrates the state of the art theory on ice failure (see Varsta, 1983), where despite bending as main failure mechanism the initial contact causes local compressive failure (crushing). The more vertical the contact area or structure becomes, the more compressive features
E-mail address: ruediger.vonbock@aalto.fi.
are included in the failure process. This underlines the importance to represent the failure in bending and compression with the same model. The work in this paper is a further development of von Bock und Polach and Ehlers (2013). In consequence to model in future ice structure interactions the compliance with cantilever beam bending tests is introduced as criterion, to establish the physical connection between tensile and compressive strength properties for horizontal loading and downward bending.
12. Beam bending simulation with state of the art model
The beam bending experiments to be reproduced numerically are conducted in the same ice sheet for which in von Bock und Polach and Ehlers (2013) the compressive and tensile material parameters are determined. The initial stress-strain relationship is governed by the elastic strain modulus, E, until yielding occurs, after which the hardening modulus, H, determines the stress-strain relationship for compression, Hc, and tension, Ht. In the plastic regime of H the stresses additionally increase by the evolution of damage according to the material model of Lemaitre and Desmorat (2001). The damage based material model of Lemaitre and Desmorat (2001) proved to represent elasto-plastic behavior including softening in tension and compression well (see von Bock und Polach and Ehlers, 2013) and is in consequence also applied on the bending failure. The damage law represents the increasing damage in the inter-granular junctions, which fail once reaching the critical damage value, dc. Fig. 3 reflects the principle stress-strain diagram, with the hardening moduli of less than 2% of
http://dx.doi.org/10.1016/j.coldregions.2015.06.003
0165-232X/© 2015 The Author. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Fig. 1. Bending failure of ice (marked crack) in front of IB Urho.
the elastic-strain modulus. Furthermore, the yield stress is low and not higher than 3% of the ultimate failure stress. More detailed information are found in von Bock und Polach and Ehlers (2013).
The material parameters developed in von Bock und Polach and Ehlers (2013) assume constant material properties throughout the thickness. Furthermore, representative material parameters are determined for each loading case with the help of numerical simulations, with which corresponding force displacement curves are generated. Consequently, those material parameters for tension and compression should suffice to model the beam bending experiment as it was done
in Fig. 4. However, the obtained numerical result shows a poor agreement with the experiments with a too weak response stiffness. The result in Fig. 4 shows that the currently available model is not suitable to characterize the mechanical behavior sufficiently.
13. Thickness dependency of strength
The hardening modulus, H, is the parameter dominating the stressstrain relationship of the model-scale ice. In the uni-axial horizontal loading cases, as in von Bock und Polach and Ehlers (2013) the through thickness distribution of H cannot be captured, since any distribution satisfying the axial stiffness, S, results in compliance with the experiments. Eq. (1) shows the basic dependency of the stiffness, S, on the cross-sectional area, A, and H as a function of the thickness coordinate z. Fig. 5a illustrates example curves for H, which are all of the same global stiffness and would produce the same force-displacement curve.
S = J H(z)dA (1)
Fig. 5b illustrates the principle stress distribution in bending for a constant hardening modulus. The bending stiffness is a function of the second moment of area where the distribution of H and the distance of each curve point, H(z), to the neutral axis which is of significance. In consequence, the distributions in Fig. 5a deliver the same axial response, but different bending stiffnesses.
Kerr and Palmer (1972) stated that the through thickness distribution of the elastic strain modulus is of significance for sea ice, which the model of von Bock und Polach and Ehlers (2013) did not account for. Furthermore, based on the suggested property distribution functions of Kerr and Palmer (1972) a four point beam bending experiment of sea ice has successfully been simulated by Ehlers and Kujala (2014). This additionally indicates the significance of the through thickness dependency of the material properties.
The through-thickness dependency of the ice properties is commonly associated with the temperature gradient over thickness. The temperature boundary conditions of the ice are, the ambient air temperature at the top and the basin water temperature at the bottom, which is analogous for sea ice (see Timco and Weeks, 2010).
The exact impact of the temperature gradient on the distribution of the model-ice properties over thickness is not found in literature. Therefore, analogy between model-scale ice and sea ice is assumed. The existence of a property or strength gradient is supported by visual observation in tensile tests, where the upper layers appeared to bear most of the load.
Fig. 3. Principle stress-strain relationship of Aalto model-scale ice defined in von Bock und Polach and Ehlers (2013).
Displacement [mm]
Fig. 4. Measured and simulated force-displacement curves for in-situ cantilever beam bending tests.
The cohesion of the grains is the principal strength parameter of ice and follows a very non-linear and exponential temperature dependency for temperatures above — 20 °C (see Fish and Zaretsky, 1997a). Fish and Zaretsky (1997a) stated that experimental data with a grain size of less than 1 mm - as the Aalto model-scale ice (see von Bock und Polach et al., 2013) - follows the same trends as Labrador ice with « 8 mm grain size diameter.
Neither for sea ice nor for model-scale ice a link between tensile, compressive and cantilever beam bending tests has been established, which is the main objective of this paper.
1.4. Paper structure and development process of the material parameters
This section introduces the process of the material parameter development which is aligned with the course of the paper. After the cantilever beam simulation with the model of von Bock und Polach and Ehlers (2013) failed, a constant hardening modulus is determined, HB, that reassembles the correct bending stiffness. HB is set as tensile hardening modulus of the top-layers Ht top. Subsequently, the distribution of the tensile hardening modulus Ht (z) is defined by Ht, top on top and the global tensile stiffness (see Eq. (1)) as boundary conditions. The characteristics ofthe elastic strain modulus, E, are assumed analogous to Ht. Once E(z) is determined, the tensile yield limit, ayt (as in von Bock
und Polach et al., 2013) and the critical damage parameter, dc, are defined, which concludes the tensile properties.
The analysis of the compressive properties starts with an effective thickness analysis to investigate their through thickness dependency. The compressive yield strength, ayc, is unknown and cannot be determined directly from experimental evidence. Therefore, it is determined iteratively based on agreement with force-displacement curves of beam bending simulations. The compliance between simulations and experiments is based on the bending force associated with the average experimental ITTC bending strength (ITTC, 2002). The 1TTC bending strength is a worldwide standard index reflecting the inherent mechanical properties of sea ice and model-scale ice and hence considered as best parameter of reference.
2. Experiments
The in-situ cantilever beam bending tests (see Fig. 4) are conducted in the same ice sheet as the tensile and compressive tests in von Bock und Polach et al. (2013), for which the material parameters have been determined in von Bock und Polach and Ehlers (2013), and with the same testing speed of 6 mm/s. The tank water contains 0.3% ethanol as weakening dopant and a freezing point between — 0.2 °C and - 0.3 °C. In total eight cantilever beam tests are conducted, which
Hardening modulus, H Bending stress,
Fig. 5. a) Various hardening modulus distributions resulting in the same axial stiffness, b) bending stress distribution for a thickness independent/constant distribution.
Table 1
Average slopes of the experimental force-displacement curves.
Kim«««««««««« Ina h
Fig. 6. Principle dimensions of the cantilever beam test.
are cut out with a stencil to reduce possible geometric variations. The geometry of the numerical model is identical to the stencil (see Fig. 6). The beam dimension comply with the guidelines of ITTC (2002) and have the thickness, h, of the ice-sheet, 25 mm, the length, l, of 175 mm (7 x thickness, h), and a width, b, of 58 mm (see also Fig. 6). The radius at the interface of root and beam, r, is 5 mm and the distance between plunger and free beam end, d, is 15 mm. The eight tests deliver an average of ab = 59.1 kPa. The measured data and accuracies are found in Appendix A, Table 7.
The experimental force-displacement curves of the bending tests do not increase linearly from zero until failure, but are arced and have changing slopes. The force levels at which the slopes change are in all tests similar and their averages are visualized in Fig. 7, with the belonging gradient values in Table 1.
The through thickness distribution of the tensile and compressive properties is of significance in the cantilever beam tests, however the tests in von Bock und Polach et al. (2013) did not explicitly focus on this. Therefore, observations are made in other ice sheets experiments, which indicate a strong through thickness dependency of the tensile properties, which is however less significant for the compressive properties. The mentioned ice sheets had a thickness of 25 mm and an average ITTC bending strength of 31 kPa and 35 kPa.
Fig. 8 displays the temperature distribution through the model-scale ice thickness, which is considered to affect the property distribution strongly (see Fish and Zaretsky, 1997b). The temperatures on top and bottom originate from air and water temperature measurements. The measurement in the middle of the ice is subjected to an estimated uncertainty of ± 3 mm for locating the probe.
Experiment section Slope [N/mm]
Total average 0.51
Section 1 0.89
Section 2 0.58
Section 3 0.27
In bending the cantilever beam experiences additional buoyancy forces when submerged, which are acting against the load-cell and increase the measured force. However, under deflection also parts of the top surface are flooded. Fig. 9 shows the flooded areas 2a and 2b, which became evident after the crack (1) was visible. However, since the crack (1) opens and becomes visible a certain time after the failure, a conservative estimate is made, by considering area 2a as flooded in all numerical calculations. In consequence, it is assumed that all encountered effects refer to the material constitution of the model-scale ice, buoyancy forces or tip flooding. Other phenomena such as hydrody-namic effects are considered to be of negligible significance due to the low plunger speed of 6 mm/s.
3. Tensile and elastic property distribution
3.1. A hardening modulus for the correct bending stiffness
The average hardening moduli for tension and compression are defined in von Bock und Polach and Ehlers (2013) as Ht = 3.07 MPa and Hc = 1.14 MPa. However, in order to simulate the correct bending stiffness an average H « 20 MPa would be required (see Fig. 10). The required hardening modulus is about seven times higher than the tensile hardening modulus. The modification of the hardening modulus would provide the correct bending stiffness, but in turn noncompliance with the tensile or compressive experiments. Furthermore, Fig. 10 shows that a linear elastic simulation based on the elastic modulus determined in von Bock und Polach et al. (2013) would be too stiff. The stiffness of the experiments lies between the simulations with the current material model (Fig. 4) and the linear elastic runs.
3.2. The distribution functions and specimen stiffness
In bending the top (as well as the bottom) layer suffers the highest strain and is hence considered responding mainly in the plastic regime with Ht. In order to achieve compliance with the experiments, it is assumed that at least the top of the model-scale ice needs to be of
Displacement [mm]
Fig. 7. Simplified representation of beam bending force-displacement curves.
Fig. 8. Measured temperature distribution over the model-scale ice thickness.
Ht« 20 MPa, while Ht in lower layers is distributed to comply with the total tensile stiffness. The target bending stiffness is achieved by mixed responses of layers with tensile stiffness properties on top, Ht, compressive properties at the bottom Hc (bottom) and elastic properties, E, around neutral axis.
Experimental observations and laboratory experience indicate that the upper layers carry significantly higher tensile load than the lower layers, which applies especially on layers above the water line. Based on this and the results in Fig. 10 the top layers are assigned with a tensile hardening modulus of Ht « 20 MPa, while the lower layers will have a significantly lower hardening modulus to comply with the overall tensile stiffness. Furthermore, it is assumed that the distributions of the elastic strain modulus and the tensile hardening modulus are qualitatively similar. The high stiffness on top and the compliance with the overall stiffness leads to a strongly parabolic distributions. This complies with the suggested temperature-strength relationship of Fish and Zaretsky (1997b) and the principal idea of the parabolic nature of ice property distributions through thickness (e.g. Kerr and Palmer, 1972 for sea ice).
The average hardening modulus in tension (see von Bock und Polach and Ehlers, 2013 with Ht, avg = 3.07 MPa) is about 7 times smaller than the average hardening modulus required to achieve the correct bending stiffness. Therefore, the tensile hardening modulus on the top of the updated model is 7 Ht, av. The two unknowns a and b are determined by the boundary conditions, BC. NB: the units of f(z) are MPa and z are in meters. Furthermore, the coordinate z of the top layer refers to z = 0.00034 m = 0.00068 m/2, which refers to the middle of the top layer thickness. The same principle applies on the bottom layer (layer thickness equals 0.68 mm, von Bock und Polach et al., 2013).
f (z) = a + b ■ z
The unknowns a and b are determined with the boundary conditions, BC.
BC 1: HWop = 7 ■ 3.07 MPa.
BC 2: Ht,avg ■ h = fya + b ■ z-0 5dz.
The first boundary condition, BC, is the set value on top and the second BC refers to the tensile stiffness. Based on the continuous function in Eq. (2) the beam properties are divided and averaged over five segments, whereas the top-segment is divided into another two in order to allocate separate properties to the part above the waterline. Furthermore, for the actual model it is assumed that the properties in the top layers above the water line are constant and equal to the highest strength values of Ht, top.
The same principle is applied on the elastic strain modulus, where the second boundary condition maintains the same plate stiffness as the average elastic modulus of 110 MPa (see von Bock und Polach et al., 2013). Additional information on the elastic strain modulus is found in Appendix B. The distribution of the tensile hardening modulus is shown in Fig. 11 and the elastic strain modulus in Fig. 12. Both Figures show the continuous distribution function, the value for each single layer and the values incorporated into the FE model. The coefficients of the distribution functions for Eq. (2) are found in Table 2.
33. The tensile yield strength
The yield strength is represented as von Mises stress and marks the threshold between the elastic and plastic domain (see Fig. 3). As in von Bock und Polach et al. (2013) the yield stress is determined by a linear elastic simulation of the infinite plate loading experiment, for the highest experimental loading case causing elastic response.
Fig. 9. Failed cantilever beam with marked crack (1), flooded surface (2a, 2b) and the dry area (3).
5 4.5 4 3.5 3 2.5 2 1.5 1
- -.Ht = 3.07 MPa, Hc = 1.14 MPa • H = 20 MPa
Linear interpolation (H = 20 MPa) E = 110 MPa Experiments
Displacement [mm]
Fig. 10. Force-displacement curves simulated with the damage material (H), linear elastic (E) and the experiments.
The new FE-model uses the distribution of the elastic modulus (Fig. 12) and delivers 1.2 kPa von Mises stress as yield strength at the bottom,
Oyt, bottom = 1.2 kPa.
Weeks and Anderson (1958) stated that sea ice has an overlaying layer from which the skeleton layer grows, which is according to Weeks and Anderson (1958) of low strength. It is assumed that such overlaying layer also exists for model scale ice. The experiments to determine the elastic strain modulus (see von Bock und Polach et al., 2013) showed that after partly plastic deformation occurs the ice still has a significant amount of elasticity left. This indicates that the layers with highest strain, i.e. only top and/or bottom — are most likely to be deformed plastically, while the rest elastically. Furthermore, the weakest layer at the bottom is loaded in tension (the ice sheet is loaded with a weight on top), in which ice tends to fail at lower strains than in compression (see also Karr and Das, 1983). This leads to the assumption that in the infinite plate test the weakest layer at the bottom is loaded in its most vulnerable stress state and is the first to deform plastically, whereas is ambiguous, whether the top - loaded in compression - remains in the elastic regime or not.
As for the elastic strain modulus and the tensile hardening modulus the yield strength is considered to follow also the function in Eq. (2). The boundary conditions, BC, for deriving the function constants are listed below for the tensile yield strength on top and
bottom (ayt, top/bottom) and the values of the coefficients, a and b, are found in Table 3.
BC 1: ayUtoP = 7.5 kPa
BC 2: Oytbottom = 1.2 kPa
The value for ayt, top is an estimate, based on increasing it iteratively while maintaining compliance with experiments (see Fig. 13, yield related peak). Fig. 13 displays the tensile test with the here derived parameters of Fig. 11 for the Lemaitre material model, the simulation with the average properties from von Bock und Polach and Ehlers (2013) and a reference measurement, which is close to the average force-displacement curve. The red circle marks the peak related to the yield strength, which indicates good agreement with the reference simulation and the measurement. A higher yield strength on the top than the stated 7.5 kPa, leads also to a significantly higher peak which would reduce the compliance with the experiments.
3.4. The critical damage parameter and tensile specimen simulations
The critical damage parameter is the damage value at which the material fails. The determination of the critical damage parameter based on specimen-simulations is affected by the yield strength, ay,
<j> <n CD a
F 0.015
Continous property function Property distribution per layer ■FE-model property distribution
Tensile hardening modulus [MPa]
Fig. 11. Distribution of the tensile hardening modulus.
Strain modululus of elasticity [MPa]
Fig. 12. Distribution of the strain modulus of elasticity.
the hardening modulus, H, and the structural stiffness i.e. the thickness and the mutual impact of the layers on each other. Therefore, it is not possible to conduct numerical specimen tests for each of the layers in Fig. 11 separately, since the adding or removing of layers affects the stress development and point of failure, i.e. the dc value. Ultimately, the number of constraints is not sufficient to determine a unique dc for each layer. Therefore, the criterion for the dc value is that the numerical tensile test can be reproduced. Simulation results show that it suffices to assign one dc value to all layers, because once the first layer fails the stress increase in the remaining cross-section leads to immediate failure through thickness (see also Fig. 13). Table 4 summarizes the obtained material properties and opposes them against the parameters from von Bock und Polach and Ehlers (2013) for through-thickness constant properties.
4. The compressive properties
4.1. The effective thickness and hardening modulus
Numerical tests are conducted to investigate the through thickness dependency of the compressive properties. The specimen thickness (i.e. in consequence the cross-sectional area) is systematically reduced and the hardening modulus increased in order to maintain the force displacement curve in von Bock und Polach and Ehlers (2013). Fig. 14 shows the hardening modulus as function of the effective specimen thickness together with images of the encountered failure patterns. At a thickness of 2/3rd of the original specimen thickness the failure pattern is still identical to the one of the original specimen. The reduction of thickness changes the local flaw concentration, which in turn affects the failure pattern and the failure origin.
A hardening modulus of 6 MPa or higher, as in Section III of Fig. 14, causes a change in the failure pattern and creates a failure behind the actual specimen in the adjacent ice plate by pushing the specimen through it, which is not observed in experiments. Consequently, the element stiffness is too high and the element or grain stiffness cannot exceed 6 MPa for the hardening modulus in order to comply
with the experiments. The maximum stiffness of the harder top layer, Hc, max, - when modeling the through thickness properties - must be smaller or equal than 5.25 times the average hardening modulus of 1.14 MPa in compression, Hc, avg (see Eq. (3)) and is therewith smaller than in tension.
Hc,max — 5-25 ■ Hc,avg- (3)
The value of Hc = 6 MPa is the theoretical extreme value and furthermore the property distribution in compression appears to be more uniform than in tension. Furthermore, if the property function Hc was similar to Ht (see Fig. 11) the compressive loaded section at the bottom had a small gradient and a nearly linear progression. As a consequence and missing additional evidence that would allow specifying the gradient, Hc is considered constant through thickness.
4.2. The compressive yield strength
The numerical analysis of the infinite plate test to determine the yield strength deliver solely information on the tensile yield strength and the compressive yield strength cannot be determined directly. Karr and Das (1983) stated for sea ice that the yield strength is significantly higher in compression than in tension and hence it is assumed that also here the yield strength is to be well above the maximum tensile value of 7.5 kPa. The compressive yield strength is assumed to be constant over the thickness in analogy to the assumptions of the hardening modulus (Section 4.1). Fig. 15 shows the force-displacement curves of a numerical compressive test for the three yield strengths of 20 kPa, 30 kPa and 40 kPa and the corresponding responses of numerical in-situ cantilever beam bending tests in Fig. 16.
The yield strength in compression is affecting the stiffness and magnitude of force at failure. The number of experimental force-displacement curves is reduced to experiments, where the loading point is located 15 mm ± 2 mm from the free edge as in the simulations. Since, the loading point of all the equally sized specimens affects the load progression.
Table 2
Coefficients of the property distribution functions.
Strain modulus of elasticity -58.53 15.16
Tensile hardening modulus -2.14 0.43
Table 3
Coefficients of the property distribution function for tensile yielding.
-617 8617
Displacement [mm] Fig. 13. Tensile simulations with highlighted yield related peak.
Table 4
Tensile material properties.
Constant material properties Thickness dependent
(von Bock und Polach material properties
and Ehlers, 2013)
Elastic strain modulus, E [MPa] 148 770-43
Hardening modulus, Ht [MPa] 3.07 21.5-0.77
Yield strength ayt [kPa] 0.5 7.5-1.2
critical damage value, dc 0.001 0.004
The bending strength, ob, of both, the experiments and the simulations is calculated following ITTC, based on the maximum forces. Fig. 17 indicates a linear correlation between the compressive yield strength and the flexural strength. Therefore, the intersection of the linear fit with the line representing the experimental average ab points at the required ayc = 23 kPa to eventually simulate the target ITTC ab. The intersection between the two lines in Fig. 17 is between ayc = 22 kPa-23 kPa. However, since the amount of water on top of the beam is rather underestimated (see Fig. 9) ayc = 23 kPa is determined.
Fig. 18a displays this representative bending force form the experiments, including upper and lower accuracy bounds. The accuracy bounds refer to the inherent uncertainties in the geometric measurements described in Appendix A. The simulation with ayc = 23 kPa fails at the target load level and can represent the material properties for the average flexural strength of the experiments. Furthermore,
Fig. 18b shows the simulation for a dry top and indicates that the flooding of the beam tip has a small impact on the force-displacement progression.
A model where the elastic strain modulus in compression is equal to the elastic strain modulus in tension requires for each layer an individual dc. Table 5 is compiling such data for the yield strength of 23 kPa for the thickness sections shown in Fig. 12. However, it must be considered that due to the large degrees of freedom for determining these parameters such as the total structural stiffness, Table 5 presents a solution for which other values may be valid as well. Nevertheless, the presented solution reflects the order of magnitude of which any other solutions are expected.
5. Summary of results
The results of the experiments (see Table 1 and Fig. 7) indicate that the beam undergoes a change in stiffness over the deformation process. This points towards layered or through thickness distributed properties, where the transition from the domain of the elastic strain modulus, E, to the hardening modulus, H, in the different layers is changing the global beam response.
Through-thickness dependent ice properties are defined, based on a qualitative analogy to distributions for sea ice according to Kerr and Palmer (1972) and the observation that the upper layers are significantly stronger than the lower layers. Fig. 19 displays the through thickness distribution of the elastic strain modulus and the tensile hardening
Fig. 14. Dependency of the hardening modulus and the failure pattern on the effective thickness in compressive tests.
Fig. 15. Force-displacement curves in compression with varying yield strength.
modulus. Both fulfill the requirements complying with the required plate bending stiffness and the axial tensile stiffness, respectively (see also Section 3).
A possible thickness dependent distribution of the compressive hardening modulus, Hc, is found to have - even for extreme values (Section 4.1) - a less strong gradient than for E or Ht and is hence considered constant over thickness. As Fig. 19 indicates that the gradient in the lower layers, which are loaded in compression, is small and furthermore constant properties ease the application in engineering. The yield strength in tension ayt has qualitatively a similar distribution as the tensile hardening modulus with values from 1.2 kPa (bottom) to 7.5 kPa (top). The compressive yield stress ayc is analogous to the Hc constant over the thickness and found to be with 23 kPa significantly higher than ayt (see Fig. 20), which complies in principle with findings for sea ice (Karr and Das, 1983).
Fig. 18 presents the successful in-situ cantilever beam simulations, which include the flooding of the beam tip, due to submergence. As afore mentioned the compressive yield strength is determined for compliance with the ITTC bending strength (see Fig. 17), which is internationally recognized as (model) ice property index.
6. Discussion
Material parameters for the Lemaitre damage model are presented for mode-scale ice of a thickness of 25 mm and an ITTC bending strength
of 59.1 kPa. The parameters are a further development of those found in von Bock und Polach and Ehlers (2013) and can represent tensile, compressive and flexural behavior.
6.1. Compliance of simulation with observations and measurements
The assumption of through thickness constant material parameters as initially done in von Bock und Polach et al. (2013) and von Bock und Polach and Ehlers (2013) delivered satisfying results for tensile and compressive specimen tests.
Measurements of air and water temperature revealed a temperature gradient between both media. The interface-layer between them is the model-scale ice sheet over which the temperature changes (see Fig. 8). The temperature gradient in the ice affects the properties in thickness direction. In addition to this, the equations found in Cox and Weeks (1983) state that the strength properties of sea ice are depending not only on the brine volume, but also on the temperature, even though Cox and Weeks (1983) did not state whether the temperature is constant or a gradient over thickness. Eventually Fish and Zaretsky (1997b) stated how the temperature and its gradient is affecting the ice properties. The simulation of the beam bending test and the derived solution for an ice property variation over thickness shows its significance to join compressive, tensile and flexural behavior in one numerical analysis.
Displacement [mm]
Fig. 16. Force-displacement of experiments and simulations for varying yield strengths.
Fig. 17. The ITTC ab of the simulations and the average experimental ab.
In Valkonen et al. (2007) it is reported that in pressure measurements of the Aalto model-scale ice interacting with a structure the pressure distribution is not homogeneous over the ice thickness. It is found that most load is transferred from the top layer. The distribution in Fig. 12 shows that the elastic strain modulus in the top layers is significantly higher than in lower layers, which also leads to a higher load transfer from those layers in the elastic domain. However, the report of Valkonen et al. (2007) does not contain detailed information on the compression progress and hence its explanatory power is limited.
6.2. The distribution function
The distribution functions for the elastic strain modulus and the tensile hardening modulus are based on the
and the maintenance of the tensile and plate bending stiffness determined in previous simulations or experiments. As Fig. 10 indicates, an average hardening modulus of H = 20 MPa appears to be in the right order of magnitude to represent the target stiffness of the beam and hence the factor seven is chosen for Ht, avg = 3.07 MPa. Furthermore, simulations are conducted with smaller beams - to reduce the computational costs - with lower Ht, max/Ht, avg ratios. As Fig. 21 shows the stiffness for factor Ht, max/Ht, avg = 5 appears to be too low and the
force displacement curve will not meet the target load for the target displacement.
The distribution functions (Figs. 11 and 12) show a large gradient over the top layers and from there a more moderate decrease towards the bottom. The high strength in the first three grain layers refers to the layers above waterline. This seems to comply with some full scale temperature measurements, e.g. Nakawo and Sinha (1981) or Petrich and Eicken (2009), which may however change with small temperature gradients or snow covers (see e.g. Gerland et al., 1999).
Data other than those of von Bock und Polach and Ehlers (2013) on the yield strength of model scale ice are not available. A graphic illustration on the idealized rheology of ice by Karr and Das (1983) indicates that the yield strength in compression is around twice as high as in tension. The ratio of compressive yield strength to tensile yield strength in the here presented model for model-scale ice is larger (« 3), but it in principle the model complies with the common phenomenon of a significantly larger yield strength in compression than in tension.
The presented ice property distributions are valid for an ice thickness of 25 mm and ab = 59 kPa. It is uncertain, whether lower bending strength values manifest solely in different properties on top or over the entire thickness. This topic will be subjected to future research.
Fig. 18. a) The target bending force together with the experiments and the simulations, including tip flooding; b) bending force for dry top surface.
Tables
Compressive material properties for varying yield strength.
Thickness/z-coordinate Elastic strain Hardening Critical damage
[mm] modulus, E [MPa] modulus, H [MPa] value, dc
0 to 2 770 0.84 0.0065
2 to 5 197 0.025
5 to 10 118 0.09
10 to 15 78 0.1
15 to 20 56 0.15
20 to 25 43 0.15
Furthermore, it is the relative impact of both cooling history and temperature distribution on the ice properties that is to be clarified.
6.3. Numerical implementation
The currently available material model in LS DYNA does not allow assigning properties for tension and compression for one single element. Instead the elements responding with compressive and tensile properties need to be assigned prior to simulations and do not change during the simulation run, e.g. after the tensile layers fail. This is the reason for the loads in Figs. 16 and 18 not fully dropping to zero. Once the tensile layers fail the upper layers of the remaining cross-section
are loaded in tension, but have assigned compressive property values. The compressive strength and strength parameters are significantly higher than the tensile ones so that they do not fail yet and the load remains at a certain level.
6.4. Tip flooding
The deflection of the cantilever beam causes flooding of the top surface. Fig. 9 shows the failed cantilever beam together with the flooded top at the point of failure. In the calculations the area 2a is considered as flooded, since it can be defined clearly and is most likely to be flooded at the time of beam failure. The flooded area on top of the beam is of significance at the point of failure reducing the load-cell force in the calculations. The image used in Fig. 9 is taken as the crack can be identified. However, the failure of the beam is likely to occur before the crack is visible, since it requires additional downward movement of the plunger to open the crack and make it visible.
The cutting process may lead to a small accumulation of model ice chips that create a slight barrier on the beam edges. The height of this artificial embankment is estimated less than 1 mm, which is not removed to avoid possible structural damage to the specimen. Nevertheless, this embankment and the surface tension of the cold water might cause a delay in wetting the surface, even though parts of the beam may be already below water level. Additionally, the plunger
-Exponential distribution, model F Constant distribution, model C (von Bock und Polach and Ehlers, 2013)
Exponential distribution, model F Constant distribution, model C (von Bock und Polach and Ehlers, 2013)
5 10 15 20 Tensile hardening modulus [MPa]
200 400 600 800
Strain modululus of elasticity [MPa]
Fig. 19. a) Through thickness distribution of the elastic strain modulus and b) the tensile hardening modulus.
E 0.01
Yield stress [kPa]
Fig. 20. Through thickness distribution of the elastic strain modulus and the tensile hardening modulus.
Displacement [mm] Fig. 21. Force-displacement curves for different Ht, max/Ht, avg ratios.
Fig. 22. Force-displacement curve for a cantilever beam with a dry top surface and a flooded tip.
pushing down the beam is acting as additional embankment and prevents immediate flooding of the area directly behind it.
Fig. 22 shows a simulated force-displacement curve together with a run with dry top surface. Fig. 22 shows the expected reduction in the maximum force due to the tip flooding. The difference in maximum force between the dry beam and the flooded tip is about 2%.
7. Conclusion
A numerical analysis of model-scale ice is presented, in which material parameters are determined to reassemble tensile, compressive and
flexural failure of model-scale ice. The material properties are defined by state of the art calculation methods in connection with experiments.
The presented analysis shows that the distribution of the ice properties through thickness is essential to model an in-situ beam bending experiment. In bending, the strains over the cross-section ofthe cantilever beam are not equal, but are high on top and bottom and zero in the neutral axis. In von Bock und Polach and Ehlers (2013) it has been shown that the material behavior of the model-scale ice may be considered following a stress-strain relationship, with a low yield point, so that the elastic strain modulus plays a minor role in the compressive and tensile tests. However, in bending, fibers near the neutral axis are subjected to small strains and respond elastic for a longer time of the
Table 6
Summary of material properties.
Thickness/z-coordinate [mm] Elastic strain modulus, E [MPa] Yield strength, ay [kPa] Hardening modulus, H [MPa] dc
Tension Compression Tension Compression Tension Compression
0 to 2 770 7.5 23 21.5 0.84 0.04 0.0065
2 to 5 197 4.5 5.2 0.025
5 to 10 118 2.5 3.0 0.09
10 to 15 78 1.8 1.8 0.1
15 to 20 56 1.4 1.16 0.15
20 to 25 43 1.2 0.76 0.15
deformation process, which assigns a high significance to the elastic strain modulus and its distribution. In consequence, the presented analysis and model show the importance of accounting for through thickness property distributions in model-scale ice bending tests. The model presented in this paper reflects the model-scale ice properties of an ice sheet of 25 mm thickness and 59.1 kPa ITTC bending strength. Table 6 summarizes the derived material parameters for elements on grain size level (0.68 mm).
Acknowledgments
The process of developing the presented model has been long and laborious and involved the time and thoughts of many people. First of all the author expresses his gratitude to the professors Soren Ehlers (NTNU) and Pentti Kujala (Aalto) for their support as well as Dr. Daniel Hilding (Dynamore), who advised the numerical modeling. Additional thanks goes to Merenkulkusaatio and the Norwegian Shipowners Foundation for financing the travels between Helsinki and Trondheim. Furthermore, the author thanks the following people for advising and contributing thoughts to the solution process: Dr. Garry Timco, Prof. Heikki Remes, Prof. Jani Romanoff, Prof. Erland Schulson, Prof. Andreas Echtermeyer, as well as Mikael Korgassar.
The beam bending experiments are conducted with a load-cell recording the force and a displacement transducer for the displacement. These two items are factory calibrated, which is cross-checked in the laboratory at working temperature. The accuracy of these devices is high and the related measurement error is considered negligible in relation to the component errors of the measurement of the specimen dimensions, relevant for calculating the flexural strength according to ITTC (2002). The measurement error is calculated as component error following NASA (2010).
The equation for the flexural strength according to ITTC (2002) is defined as:
ct b =-
6■ Fl
The measurements of the dimensions beam length, l, beam width, b, and ice thickness, h, are taken by measurement tapes and calipers and hence the estimated accuracies, A, are as follows:
• A( = 2 mm
• Ab = 1 mm
• Ah = 1 mm
Appendix A. Measurement data and measurement uncertainty
The data of the cantilever beam experiments conducted are compiled in Table 7. The measurement uncertainty in Table 7 is determined according to ITTC (2005). The bending force is recorded with a factory calibrated load-cell, which is cross-checked in the laboratory at working temperature. The obtained measurement accuracy is of 99.9996%. The measurements of the dimensions beam length, l, beam width, b, and ice thickness, h, are taken by measurement tapes and calipers. It is acknowledged that the high pressure sensitivity of the material adds uncertainty and therefore the estimated accuracies, A, are as follows:
• Ai = 2 mm
• Ab = 1 mm
• Ah =1 mm
The equation for the flexural strength according to ITTC (2002) is defined as:
ct b = -
The total uncertainty, Atot, is defined according to ITTC (2005):
A tot =
/fdCTbAFiVCTbAiiV^CTblv+ftb
Atot, average '■
±4.73 kPa.
The inherent uncertainty of each particular measurement is found in Table 7.
Table 7
Data of the experimental cantilever beam tests including average values. The beam length, l, is the difference between the beam length after failure, lb, and the distance between plunger and beam tip, dp.
Test h [mm] b [mm] lb [mm] dp [mm] l [mm] F[N] CTb[kPa] Atot [kPa]
IB 1 26 58 173 19 154 2.64 62.13 4.96
IB 2 27 58 170 25 145 2.74 56.47 4.37
IB3 24 59 165 19.5 145.5 2.13 54.78 4.72
IB 4 26 58 175 17 158 2.55 61.55 4.91
IB 5 25 60 174 15 159 2.41 61.26 5.06
IB 6 27 56 172 19 153 2.19 49.36 3.82
IB 7 26 56.5 177 17 160 3.06 76.88 6.14
IB 8 27 59 173 15 158 2.29 50.40 3.88
Average: 59.10 4.73
The total error, Atot, is defined as
A dCTb dCTb dCTb A Atot =-sr Al + -db Ab + -gh~ Ah
•tot, average ;
6.3 kPa.
Appendix B. Mesh sensitivity of the elastic strain modulus
The significance of the through thickness properties requires a revision of the previously determined elastic strain modulus, since it was previously determined with a mesh size of 25 mm (E = 148 MPa). In the numerical specimen tests for tension, compression and bending the mesh size is aligned to the grain size of 0.68 mm. Due to the high computational costs and the low significance of the linear elastic regime for the tensile and compressive test simulations the mesh size convergence was not investigated in von Bock und Polach et al. (2013). However, in the beam bending experiments the actual value of the elastic strain modulus is of higher significance than in tension or compression, because layers near the neutral axis experience little strain and may respond linear-elastically until failure. Furthermore, access to the supercomputer Vilje1 allowed the employment of larger models with a finer discretization. The deflection test of the infinite plate on elastic foundation from von Bock und Polach et al. (2013) is repeated with different mesh discretization, but unchanged material properties. Fig. 23 displays the sensitivity of the displacement of the infinite plate for different element sizes and a constant elastic modulus. An element size correlating to the grain size as used in the beam bending simulations could not be tested due to memory limitations. However, Fig. 23 indicates convergence starting from 6 mm element size.
The elastic strain modulus is ultimately found to be 110 MPa for element size of 6.25 mm to maintain the same displacement, i.e. 0.055 mm, as in the experiments and the 25 mm element size model with 148 MPa as in Fig. 23.
Appendix C. Numerical model
The material law to represent the flexural failure is based on the damage law of Lemaitre and Desmorat (2001) and is analogous to the
1 https://www.notur.no/hardware/vilje.
E 0.052 JE
1 °.°5 iE
o 0.048 J2 ü CO
¡5 0.046
Element size [mm]
Fig. 23. Mesh size impact on the bending displacement of an infinite ice sheet.
one used in von Bock und Polach and Ehlers (2013). The analysis is performed in LS-DYNA with the material model 153-Damage_3. The cylindrical plunger pushing the beam downwards consists of shell elements and is modeled as rigid material. The contact is modeled as AUTOMATIC_NODES_ TO_SURFACE, while the plunger speed at contact is modeled with a smooth velocity profile to avoid high impact loads at the initial contact. The buoyancy effect is modeled with Motion_Nodes, where displacements cause a reaction force like springs. Since the material model allows only the parameters for either compression or tension, the upper layers are created as tensile layers and the lower ones as compressive layers. The neutral axis, na (see Fig. 6), as their interface is calculated based on the hardening moduli.
References
Cox, G., Weeks, W., 1983. Equations for determining the gas and brine volumes in sea ice
samples. J. Glaciol. 29, 306-316. Ehlers, S., Kujala, P., 2014. Optimization-based material parameter identification for the numerical simulation of sea ice in four-point bending. J. Eng. Marit. Environ. 228 (1), 70-80.
Enkvist, E., 1972. On the ice resistance encountered by ships operating in the continuous mode of ice breaking. Technical report. The Swedish Academy of Engineering Sciences in Finland.
Fish, A.M., Zaretsky, Y.K., 1997a. Ice strength as an function of hydrostatic pressure and
temperature. Technical report, CRREL. Fish, A.M., Zaretsky, Y.K., 1997b. Temperature effect on strength of ice under triaxial
compression vol. 2 pp. 415-422. Gerland, S., Winther, J.-G., Orbaek, J., Ivanov, B., 1999. Physical properties, spectral reflectance and thickness development of first year fast ice in Kongsfjorden, Svalbard. Polar Res. 18.
ITTC, 2002. Ice property measurements, 7.5-02-04-02. http://ittc.sname.org.
ITTC, 2005. Experimental uncertainty analysis for ship resistance in ice tank testing, 7.502-04-02.5. http://ittc.sname.org. Karr, D., Das, S., 1983. Ice strength in brittle and ductile failure modes. J. Struct. Eng. 109, 2802-2811.
Kerr, A., Palmer, W., 1972. The deformations and stresses in floating ice plates. Acta Mech. 15,57-72.
Lemaitre, J., Desmorat, R., 2001. Handbook of Materials Behavior Models. Academic Press. Lindqvist, G., 1989. A straightforward method for calculation of ice resistance of ships. The 10th International Conference on Port and Ocean Engineering Under Arctic Conditions 2, pp. 722-735. Nakawo, M., Sinha, N., 1981. Growth rate and salinity profile of first-year sea ice in the
high arctic. J. Glaciol. 27 (96). NASA, 2010. Measurement uncertainty analysis principles and methods. NASA-HDBK-8739.19-3. National Aeronautics and Space Administration (NASA), Washington DC, p. 20546.
Petrich, C., Eicken, H., 2009. Sea ice, chapter growth, structure and properties of sea ice.
Wiley-Blackwell, Oxford, UK. Timco, G., Weeks, W., 2010. A review of the engineering properties of sea ice. Cold Reg.
Sci. Technol. 60 (2), 107-129. Valanto, P., 2001. The resistance of ships in level ice. SNAME Transactions 119 pp. 53-83. Valkonen, J., Harris, T., Izumiyama, K., Takimoto, T., Ritch, R., 2007. WP10 field trials vol. 2. Varsta, P., 1983. On the mechanics of ice load on ships in level ice in the Baltic Sea PhD
thesis. Helsinki University of Technology. von Bock und Polach, R., Ehlers, S., 2013. Model scale ice — part b: numerical model. Cold Reg. Sci. Technol. 94.
von Bock und Polach, R., Ehlers, S., Kujala, P., 2013. Model scale ice — part a: experiments.
Cold Reg. Sci. Technol. 94, 53-60. von Bock und Polach, R., Ehlers, S., Erikstad, S., 2014. A decision-based design approach for ships operating in open water and ice. J. Ship Prod. Des. 30 (3), 1-11. http://dx.doi. org/10.5957/JSPD.30.3.140001. Weeks, W., Anderson, D., 1958. An experimental study of strength of young sea ice. Trans. Am. Geophys. Union 39 (4).
