tmain.tex - seaice-experiments - sea ice experiments using Granular.jl
(HTM) git clone git://src.adamsgaard.dk/seaice-experiments
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tmain.tex (20080B)
---
1 %\documentclass[11pt,a4paper]{article}
2 \documentclass[11pt]{article}
3
4 %\usepackage{a4wide}
5 \usepackage[margin=1.4in]{geometry}
6
7 \usepackage{graphicx}
8 %\usepackage[german, english]{babel}
9 %\usepackage{tabularx}
10 %\usepackage{cancel}
11 %\usepackage{multirow}
12 %\usepackage{supertabular}
13 %\usepackage{algorithmic}
14 %\usepackage{algorithm}
15 %\usepackage{amsthm}
16 %\usepackage{float}
17 %\usepackage{subfig}
18 \usepackage{rotating}
19 \usepackage{amsmath}
20
21 \usepackage[T1]{fontenc} % Font encoding
22 \usepackage{charter} % Serif body font
23 \usepackage[charter]{mathdesign} % Math font
24 \usepackage[scale=0.9]{sourcecodepro} % Monospaced fontenc
25 \usepackage[lf]{FiraSans} % Sans-serif font
26 \usepackage{textgreek}
27
28 \usepackage{listings}
29
30 \usepackage{hyperref}
31
32 \usepackage{soul} % for st strikethrough command
33
34 %\usepackage[round]{natbib}
35 \usepackage[natbib=true, style=authoryear, bibstyle=authoryear-comp,
36 maxbibnames=10,
37 maxcitenames=2, backend=bibtex8]{biblatex}
38 \bibliography{/Users/ad/articles/own/BIBnew.bib}
39
40
41 \begin{document}
42
43 \title{Discrete-element simulation of sea-ice mechanics:\\
44 Parameterization of contact mechanics}
45
46 \author{\small Anders Damsgaard}
47 \date{\small \today}
48
49 \maketitle
50
51 \section{Rationale}
52 Lagrangian parameterizations of sea ice offer several advantages to Eulerian
53 continuum methods. The continuum assumption of sea-ice mechanics becomes
54 increasingly questionable as spatial resolution and cell size approaches the
55 size of ice floes. Ice-floe scale discretizations are natural for Lagrangian
56 models, which additionally offer the convenience of being able to handle
57 arbitrary sea-ice concentrations. This is likely to improve model performance
58 in ice-marginal zones with strong advection. Furthermore, phase transitions in
59 granular rheology around the jamming limit, such as observed when sea ice moves
60 through geometric confinements, includes sharp thresholds in effective
61 viscosity, which are typically ignored in Eulerian models.
62 Granular jamming is a stochastic process dependent on having the right grains in
63 the right place at the right time, and the jamming likelihood over time can be
64 described by a probabilistic model
65 \citep[e.g.][]{Tang2009}. Difficult to parameterize in continuum formulations
66 \citep[e.g.][]{Rallabandi2017}, jamming comes natural to dense granular systems
67 simulated in a Lagrangian framework, and is a very relevant process complicating
68 sea-ice transport through narrow straits \citep[e.g.][]{Kwok2010}.
69 However, computational performance is a central concern for coupled
70 ocean-atmosphere models, with traditional discrete-element method (DEM)
71 formulations being too costly for global-scale models. Here we investigate the
72 behavior of a Lagrangian sea-ice model and assess how mechanical
73 parameterizations of differing complexity affect collective behavior at scales
74 much larger than the individual ice floes.
75
76
77 \section{Contact models}
78 The presented experiments compare jamming behavior between two differing
79 ice-floe contact models. Common to both models, the resistive force
80 $\boldsymbol{f}_\text{n}$ to axial compressive strain between to cylindrical ice
81 floes $i$ and $j$ is provided by (Hookean) linear-elasticity based on the
82 overlap distance $\boldsymbol{\delta}_\text{n}$:
83 \begin{equation}
84 \boldsymbol{f}_\text{n}^{ij} =
85 A^{ij} E^{ij} \boldsymbol{\delta}_\text{n}^{ij}
86 \quad \text{when} \quad
87 0 > |\boldsymbol{\delta}_\text{n}^{ij}| \equiv |\boldsymbol{x}^i -
88 \boldsymbol{x}^j| - (r^i + r^j)
89 \label{eq:f_n}
90 \end{equation}
91 where the contact cross-sectional area $A^{ij} = R^{ij} \min(T^i, T^j)$ is
92 determined by the harmonic mean $R^{ij} = 2 r^i r^j/(r^i + r^j)$ of the ice-floe
93 radii $r^i$ and $r^j$, as well as the smallest of the involved ice-floe
94 thickness $T^i$ and $T^j$. The linear elastic force resulting from axial strain
95 of a distance $|\boldsymbol{\delta}_\text{n}^{ij}|$ of the contact is scaled by
96 the harmonic mean of Young's modulus $E^{ij}$. The stiffness scaling based on
97 Young's modulus is scale invariant \citep[e.g.][]{Obermayr2013}, and is based on
98 the principle that the elastic stiffness of the ice itself is constant,
99 regardless of ice-floe size.
100
101 Nonlinear elasticity models based on Hertzian contact mechanics may
102 alternatively be applied to determine the stresses resulting from contact
103 compression \citep[e.g.][]{Herman2013}. However, with nonlinear models the
104 numerical stability of the explicit temporal integration scheme will depend on
105 the stress and packing state of the granular assemblage, and will under
106 localized compressive extremes require very small discretization in time.
107
108 As demonstrated later, models based compressive strength alone are not
109 sufficient to cause granular jamming. We explore two different additions to the
110 contact model presented in Eq.~\ref{eq:f_n}. The first approach is typical to
111 DEM models and is based on resolving including shear resistance through
112 tangential (contact parallel) elasticity and associated Coulomb frictional
113 limits. An alternative approach, fundamentally complementary to compressive
114 elasticity and shear friction, is tensile strength of ice-floe contacts which
115 leads to a cohesive bulk granular rheology.
116
117 \subsection{Contact-parallel elasticity with Coulomb friction}
118 The contact-parallel (tangential) elasticity is resolved by determining the
119 contact travel distance $\boldsymbol{\delta}_\text{t}$ on the contact plane for
120 the duration of the contact $t_\text{c}$:
121 \begin{equation}
122 \boldsymbol{\delta}_\text{t} = \int_0^{t_\text{c}}
123 \left[
124 (\boldsymbol{v}^i - \boldsymbol{v}^j) \cdot \hat{\boldsymbol{t}}^{ij} -
125 R^{ij} \left(\omega^i + \omega^j\right)
126 \right]
127 \label{eq:d_t}
128 \end{equation}
129 where $\boldsymbol{v}$ and $\omega$ denotes linear and rotational velocity,
130 respectively.
131 The deformation distance $\boldsymbol{\delta}_\text{t}$ is corrected for contact
132 rotation over the lifetime of the interaction, and is used to determine the
133 contact-tangential elastic force:
134 \begin{equation}
135 \boldsymbol{f}_\text{t} = \frac{E^{ij} A^{ij}}{R^{ij}}
136 \frac{2 (1 - \nu^2)}{(2 - \nu)(1 + \nu)} \boldsymbol{\delta}_\text{t}
137 \end{equation}
138 with $\nu$ being the dimensionless Poisson's ratio characteristic for the
139 material. However, the frictional force above is restricted to the
140 Coulomb-frictional limit, relating the magnitude of the normal force to the
141 magnitude of the tangential force:
142 \begin{equation}
143 |\boldsymbol{f}_\text{t}| \leq \mu |\boldsymbol{f}_\text{n}|
144 \label{eq:coulomb_friction}
145 \end{equation}
146 The Coulomb-frictional coefficient $\mu$ may be parameterized as state
147 dependent, for example with different values for sliding and stationary
148 contacts. However, for the following a single and constant value is used for
149 the frictional coefficient $\mu$.
150
151 In the case of slip ($|\boldsymbol{f}_\text{t}| > \mu
152 |\boldsymbol{f}_\text{n}|$) the length of the contact-parallel travel path
153 $\boldsymbol{\delta}_\text{t}$ is reduced to be consistent to be consistent with
154 the Coulomb limit. This accounts for the tangential contact plasticity and
155 irreversible work associated with sliding.
156
157 Since the above model of tangential shear resistance is based on deformation
158 distance on the inter-floe contact plane, it requires solving for ice-floe
159 rotational kinematics of each ice floe and a bookkeeping algorithm for storing
160 contact histories.
161
162 \subsection{Tensile contact strength}
163 Cohesion is introduced by parameterizing resistance to extension beyond the
164 overlap distance between a pair of ice floes (i.e.\ $\delta_\text{n}^{ij} > 0$).
165 In real settings tensile strength may be provided by refreezing processes at the
166 ice-floe interface or mechanical ridging. An ideal formulation of bond
167 deformation includes resistance to bond tension, shear, twist, and rolling
168 \citep[e.g.][]{Potyondy2004, Obermayr2013, Herman2016}. However, for this study
169 we explore the possibility of using bond \emph{tension} alone as a mechanical
170 component contributing to bulk granular shear strength.
171
172 Tensile strength is parameterized by applying Eq.~\ref{eq:f_n} for the extensive
173 regime ($\delta_\text{n} > 0$). Eq.~\ref{eq:f_n} is applied until the tensile
174 stress exceeds the tensile strength $\sigma_\text{c}$ defined for the bonds.
175 The implementation allows for time-dependent tensile strengthening based on
176 contact age, but for the following we parameterize the bonds to obtain full
177 tensile strength as soon as a pair of ice floes first undergo compression
178 ($\delta_\text{n} < 0$).
179
180
181 \subsection{Implementation}
182 The intent of this work is to develop and apply an offline sea-ice dynamical
183 code in order to explore strengths and limitations of different methods related
184 to sea-ice mechanics. Once the appropriate parameterization has been identified
185 it will be built in to the GFDL model of coupled ocean/atmosphere dynamics.
186
187 The presented experiments are performed with the
188 \texttt{SeaIce.jl}\footnote{\url{https://github.com/anders-dc/SeaIce.jl}} Julia
189 package for offline simulations of sea-ice mechanics, which includes the above
190 parameterizations. Furthermore, the model includes the option to include linear
191 (Newtonian) contact viscosity in the normal and tangential components, but an
192 analysis of these is not included here. Simulation-specific scripts are
193 provided in the repository
194 \texttt{SeaIce-experiments}\footnote{\url{https://github.com/anders-dc/SeaIce-experiments}}.
195
196 The offline model implementation identifies ice-floe contacts and placement in
197 the ocean/atmosphere grid through a cell-based spatial discretization, which
198 reduces the computational overhead significantly (Fig.~\ref{fig:performance},
199 and additionally mirrors the current Lagrangian ice-berg implementation in the
200 GFDL ocean model.
201
202 \begin{figure}
203 \begin{center}
204 \includegraphics[width=0.7\textwidth]{graphics/profiling-cpu.pdf}
205 \end{center}
206 \caption{\label{fig:performance} Single-threaded performance of the DEM
207 algorithm with all-to-all contact search, or a contact search using spatial
208 decomposition based on the ocean/atmosphere grid (Profiling script:
209 \url{https://github.com/anders-dc/SeaIce.jl/blob/master/test/profiling.jl}).}
210 \end{figure}
211
212 The experiments rely on pseudo-random number generation for generating ice-floe
213 particle size distributions (PSDs). In order to assess the statistical
214 probability of granular jamming we seed the system with different values and
215 repeat each experiment ten times.
216
217 For the presented runs the ice floes are forced with a uniform and constant wind
218 field oriented from north to south, with a simplified channel-like constriction
219 in the middle. The ocean also flows from north to south but is constrained by
220 the basin geometry. The spatial velocity pattern is determined by a stream
221 function under mass conservation, meaning that currents increase as
222 the channel narrows. Ice floes are initially placed in a pseudo-random
223 arrangement north of the channel. All parameters are kept constant between the
224 experiment sets unless explicitly noted.
225
226 \section{Results}
227
228 \subsection{Experiment set 1}
229
230 Increasing frictional coefficients, increasing tensile strength. Constant
231 uniform size distribution.
232 Generated by running \texttt{`make`} from
233 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-comparison}
234 Animation with \texttt{seed=1} can be found at:
235 \url{https://youtu.be/JmWvoi9Z6Zk}
236
237 \begin{sidewaysfigure}
238 \begin{center}
239 \includegraphics[width=0.99\textheight]%
240 {graphics/cohesion-friction-comparison.png}
241 \includegraphics[width=0.99\textheight]%
242 {{graphics/cohesion-friction-comparison-jam_fraction.png}}
243 \end{center}
244 \caption{\label{fig:set1}%
245 \textbf{Experiment set 1}:
246 Comparison between the frictional and cohesive model with different
247 Coulomb-frictional coefficient values and tensile strengths.
248 First row: Ice fluxes over time with increasing friction ($\mu = [0.0, 0.1,
249 0.2, 0.3, 0.4]$).
250 Second row: Ice fluxes over time with increasing tensile strength
251 ($\sigma_\text{c} = [10, 100, 200, 400, 800]$ kPa).
252 Third row: Fraction of runs jammed over time with increasing friction.
253 Fourth row: Fraction of runs jammed over time with increasing tensile
254 strength.
255 }
256 \end{sidewaysfigure}
257
258 While the frictionless runs rarely jam, we observe that increasing either
259 friction or cohesion greatly increases the likelihood of jamming
260 (Fig.~\ref{fig:set1}). Under cohesive runs with large tensile strengths the ice
261 floes north of the assemblage behave like a single rigid block, while the
262 corresponding frictional runs show more spread in the jam-fraction probability.
263 At intermediate to low frictional coefficients ($\mu \leq 0.3$ or
264 $\sigma_\text{c} \leq$ 200 kPa) there is good agreement between frictional and
265 cohesive models, where the runs with $\mu = 0.3$ most closely resemble the
266 jam-fraction distribution of $\sigma_\text{c} = 200$ kPa. We base further runs
267 on comparing the differences between these two parameter sets under varying
268 circumstances.
269
270 \subsection{Experiment set 2}
271 Constant friction coefficient ($\mu = 0.3$) for frictional runs, constant
272 tensile strength ($\sigma_\text{c} = 200$ kPa) for cohesive runs. Variable
273 channel width. Monodisperse ice-floe size distribution.
274 Generated by running \texttt{`make`} from
275 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-width-monodisperse-comparison}.
276 Animation with \texttt{seed=1} can be found at:
277 \url{https://youtu.be/tPMSxdw1UZ8}.
278
279 \begin{sidewaysfigure}
280 \begin{center}
281 \includegraphics[width=0.99\textheight]%
282 {graphics/cohesion-friction-width-monodisperse-comparison.png}
283 \includegraphics[width=0.99\textheight]%
284 {{graphics/cohesion-friction-width-monodisperse-comparison-jam_fraction.png}}
285 \end{center}
286 \caption{\label{fig:set2}%
287 \textbf{Experiment set 2}:
288 Comparison between the frictional ($\mu = 0.3$ and $\sigma_\text{c}$ = 0
289 kPa) and cohesive ($\mu = 0$ and $\sigma_\text{c}$ = 200 kPa) models with a
290 monodisperse (single) ice-floe size and increasing channel width.
291 First row: Ice fluxes over time with friction ($\mu = 0.3$) and increasing
292 channel width.
293 Second row: Ice fluxes over time with cohesion ($\sigma_\text{c} = 200$ kPa)
294 and increasing channel width.
295 Third row: Fraction of frictional runs jammed over time.
296 Fourth row: Fraction of cohesive runs jammed over time.
297 }
298 \end{sidewaysfigure}
299
300 The jamming probability of two-dimensional granular assemblages under gravity is
301 well described in the literature, especially for monodisperse (same-sized)
302 granular materials \citep[e.g.][]{Beverloo1961, To2001, Tang2009}.
303 While the ocean/atmosphere forcing is different from the gravity drag in the
304 silo flows in these publications, we observe the same relationship between
305 conduit width and grain size in monodisperse granular assemblages as previously
306 reported. The consistency between the frictional and cohesive model
307 (Fig.~\ref{fig:set2}) indicates that the simpler cohesive model is performing
308 well in this setting.
309
310 \subsection{Experiment set 3}
311 Constant friction coefficient ($\mu = 0.3$) for frictional runs, constant
312 tensile strength ($\sigma_\text{c} = 200$ kPa) for cohesive runs. Variable
313 channel width. Constant and uniformly-distributed ice-floe size distribution.
314 Generated by running \texttt{`make`} from
315 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-width-comparison}.
316
317 \begin{sidewaysfigure}
318 \begin{center}
319 \includegraphics[width=0.99\textheight]%
320 {graphics/cohesion-friction-width-comparison.png}
321 \includegraphics[width=0.99\textheight]%
322 {{graphics/cohesion-friction-width-comparison-jam_fraction.png}}
323 \end{center}
324 \caption{\label{fig:set3}%
325 \textbf{Experiment set 3}:
326 Comparison between the frictional ($\mu = 0.3$ and $\sigma_\text{c}$ = 0
327 kPa) and cohesive ($\mu = 0$ and $\sigma_\text{c}$ = 200 kPa) models with a
328 uniform ice-floe size distribution and increasing channel width.
329 First row: Ice fluxes over time with friction ($\mu = 0.3$) and increasing
330 channel width.
331 Second row: Ice fluxes over time with cohesion ($\sigma_\text{c} = 200$ kPa)
332 and increasing channel width.
333 Third row: Fraction of frictional runs jammed over time.
334 Fourth row: Fraction of cohesive runs jammed over time.
335 }
336 \end{sidewaysfigure}
337
338 While the monodisperse ice-floe size used in Set 2 provides a comparative basis
339 to the literature on granular jamming, we have for Set 3 used a wide grain size
340 distribution and varied the channel width as before. We observe that the width
341 of ice-floe size distribution decreases the likelihood of jamming, and that
342 these runs are less likely to jam than the corresponding monodisperse ensemble
343 runs. Again we see a good correspondence between the frictional and cohesive
344 model runs (Fig.~\ref{fig:set3}).
345
346 \subsection{Experiment set 4}
347 Constant friction coefficient ($\mu = 0.3$) for frictional runs, constant
348 tensile strength ($\sigma_\text{c} = 200$ kPa) for cohesive runs. Constant
349 channel width. Different widths in power-law distributed ice-floe size
350 distribution.
351 Generated by running \texttt{`make`} from
352 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-psd-comparison}.
353 Animation with \texttt{seed=1} can be found at:
354 \url{https://youtu.be/I9P8XsA1S0I}
355
356 \begin{sidewaysfigure}
357 \begin{center}
358 \includegraphics[width=0.99\textheight]%
359 {graphics/cohesion-friction-psd-comparison.png}
360 \includegraphics[width=0.99\textheight]%
361 {{graphics/cohesion-friction-psd-comparison-jam_fraction.png}}
362 \end{center}
363 \caption{\label{fig:set4}%
364 \textbf{Experiment set 4}:
365 Comparison between the frictional ($\mu = 0.3$ and $\sigma_\text{c}$ = 0
366 kPa) and cohesive ($\mu = 0$ and $\sigma_\text{c}$ = 200 kPa) models with a
367 increasingly wide power-law ice-floe size distribution and constant channel
368 width.
369 First row: Ice fluxes over time with friction ($\mu = 0.3$) and increasing
370 ice-floe size span.
371 Second row: Ice fluxes over time with cohesion ($\sigma_\text{c} = 200$ kPa)
372 and increasing ice-floe size span.
373 Third row: Fraction of frictional runs jammed over time.
374 Fourth row: Fraction of cohesive runs jammed over time.
375 }
376 \end{sidewaysfigure}
377
378 Ice-floe distributions in natural sea-ice settings have been observed to follow
379 a power-law probability distribution with an exponent of $\alpha \approx -1.8$
380 \citep[e.g.][]{Steer2008, Herman2013}. With constant channel size and an
381 increasing width of the ice-floe size distribution, we observe that the
382 probability of jamming decreases, and observe good agreement between the
383 frictional and cohesive model runs.
384
385 \section{Summary}
386 We construct a flexible discrete-element framework for simulating Lagrangian
387 sea-ice dynamics at the ice-floe scale, forced by ocean and atmosphere velocity
388 fields. While frictionless contact models based on tensile stiffness alone are
389 very unlikely to jam, we describe two different approaches based on friction and
390 tensile strength, both resulting in increased bulk shear strength of the
391 granular assemblage. We demonstrate that the discrete-element approach is able
392 to undergo dynamic jamming when forced through an idealized confinement, where
393 the probability of jamming is determined by the channel width, ice-floe size,
394 and ice-floe size variability. The frictionless but cohesive contact model can
395 with certain tensile strength values display jamming behavior which on the large
396 scale is highly similar to a model with contact friction and ice-floe rotation.
397 The frictionless and cohesive model is likely an ideal candidate for coupled
398 computations due to its reduced algorithmic complexity.
399
400 \printbibliography{}
401
402 \end{document}