On the mechanical behavior of compacted pack ice : A theoretical and numerical investigation

Pack ice is an aggregate of ice floes drifting on the sea surface. The forces controlling the motion and deformation of pack ice are air and water drag forces, sea surface tilt, Coriolis force and the internal force due to the interaction between ice floes. In this thesis, the mechanical behavior of compacted pack ice is investigated using theoretical and numerical methods, focusing on the three basic material properties: compressive strength, yield curve and flow rule. A high-resolution three-category sea ice model is applied to investigate the sea ice dynamics in two small basins, the whole Gulf Riga and the inside Pärnu Bay, focusing on the calibration of the compressive strength for thin ice. These two basins are on the scales of 100 km and 20 km, respectively, with typical ice thickness of 10-30 cm. The model is found capable of capturing the main characteristics of the ice dynamics. The compressive strength is calibrated to be about 30 kPa, consistent with the values from most large-scale sea ice dynamic studies. In addition, the numerical study in Pärnu Bay suggests that the shear strength drops significantly when the ice-floe size markedly decreases. A characteristic inversion method is developed to probe the yield curve of compacted pack ice. The basis of this method is the relationship between the intersection angle of linear kinematic features (LKFs) in sea ice and the slope of the yield curve. A summary of the observed LKFs shows that they can be basically divided into three groups: intersecting leads, uniaxial opening leads and uniaxial pressure ridges. Based on the available observed angles, the yield curve is determined to be a curved diamond. Comparisons of this yield curve with those from other methods show that it possesses almost all the advantages identified by the other methods. A new constitutive law is proposed, where the yield curve is a diamond and the flow rule is a combination of the normal and co-axial flow rule. The non-normal co-axial flow rule is necessary for the Coulombic yield constraint. This constitutive law not only captures the main features of forming LKFs but also takes the advantage of avoiding overestimating divergence during shear deformation. Moreover, this study provides a method for observing the flow rule for pack ice during deformation.

Asteroid identification using statistical orbital inversion methods

An efficient and statistically robust solution for the identification of asteroids among numerous sets of astrometry is presented. In particular, numerical methods have been developed for the short-term identification of asteroids at discovery, and for the long-term identification of scarcely observed asteroids over apparitions, a task which has been lacking a robust method until now. The methods are based on the solid foundation of statistical orbital inversion properly taking into account the observational uncertainties, which allows for the detection of practically all correct identifications. Through the use of dimensionality-reduction techniques and efficient data structures, the exact methods have a loglinear, that is, O(nlog(n)), computational complexity, where n is the number of included observation sets. The methods developed are thus suitable for future large-scale surveys which anticipate a substantial increase in the astrometric data rate. Due to the discontinuous nature of asteroid astrometry, separate sets of astrometry must be linked to a common asteroid from the very first discovery detections onwards. The reason for the discontinuity in the observed positions is the rotation of the observer with the Earth as well as the motion of the asteroid and the observer about the Sun. Therefore, the aim of identification is to find a set of orbital elements that reproduce the observed positions with residuals similar to the inevitable observational uncertainty. Unless the astrometric observation sets are linked, the corresponding asteroid is eventually lost as the uncertainty of the predicted positions grows too large to allow successful follow-up. Whereas the presented identification theory and the numerical comparison algorithm are generally applicable, that is, also in fields other than astronomy (e.g., in the identification of space debris), the numerical methods developed for asteroid identification can immediately be applied to all objects on heliocentric orbits with negligible effects due to non-gravitational forces in the time frame of the analysis. The methods developed have been successfully applied to various identification problems. Simulations have shown that the methods developed are able to find virtually all correct linkages despite challenges such as numerous scarce observation sets, astrometric uncertainty, numerous objects confined to a limited region on the celestial sphere, long linking intervals, and substantial parallaxes. Tens of previously unknown main-belt asteroids have been identified with the short-term method in a preliminary study to locate asteroids among numerous unidentified sets of single-night astrometry of moving objects, and scarce astrometry obtained nearly simultaneously with Earth-based and space-based telescopes has been successfully linked despite a substantial parallax. Using the long-term method, thousands of realistic 3-linkages typically spanning several apparitions have so far been found among designated observation sets each spanning less than 48 hours.

Comparative study of different formulations and solution methods for two phase flow in saturated porous media

Six models (Simulators) are formulated and developed with all possible combinations of pressure and saturation of the phases as primary variables. A comparative study between six simulators with two numerical methods, conventional simultaneous and modified sequential methods are carried out. The results of the numerical models are compared with the laboratory experimental results to study the accuracy of the model especially in heterogeneous porous media. From the study it is observed that the simulator using pressure and saturation of the wetting fluid (PW, SW formulation) is the best among the models tested. Many simulators with nonwetting phase as one of the primary variables did not converge when used along with simultaneous method. Based on simulator 1 (PW, SW formulation), a comparison of different solution methods such as simultaneous method, modified sequential and adaptive solution modified sequential method are carried out on 4 test problems including heterogeneous and randomly heterogeneous problems. It is found that the modified sequential and adaptive solution modified sequential methods could save the memory by half and as also the CPU time required by these methods is very less when compared with that using simultaneous method. It is also found that the simulator with PNW and PW as the primary variable which had problem of convergence using the simultaneous method, converged using both the modified sequential method and also using adaptive solution modified sequential method. The present study indicates that pressure and saturation formulation along with adaptive solution modified sequential method is the best among the different simulators and methods tested.

Frequency Domain Based Solution for Certain Class of Wave Equations: An exhaustive study of Numerical Solutions

The paper discusses the frequency domain based solution for a certain class of wave equations such as: a second order partial differential equation in one variable with constant and varying coefficients (Cantilever beam) and a coupled second order partial differential equation in two variables with constant and varying coefficients (Timoshenko beam). The exact solution of the Cantilever beam with uniform and varying cross-section and the Timoshenko beam with uniform cross-section is available. However, the exact solution for Timoshenko beam with varying cross-section is not available. Laplace spectral methods are used to solve these problems exactly in frequency domain. The numerical solution in frequency domain is done by discretisation in space by approximating the unknown function using spectral functions like Chebyshev polynomials, Legendre polynomials and also Normal polynomials. Different numerical methods such as Galerkin Method, Petrov- Galerkin method, Method of moments and Collocation method or the Pseudo-spectral method in frequency domain are studied and compared with the available exact solution. An approximate solution is also obtained for the Timoshenko beam with varying cross-section using Laplace Spectral Element Method (LSEM). The group speeds are computed exactly for the Cantilever beam and Timoshenko beam with uniform cross-section and is compared with the group speeds obtained numerically. The shear mode and the bending modes of the Timoshenko beam with uniform cross-section are separated numerically by applying a modulated pulse as the shear force and the corresponding group speeds for varying taper parameter in are obtained numerically by varying the frequency of the input pulse. An approximate expression for calculating group speeds corresponding to the shear mode and the bending mode, and also the cut-off frequency is obtained. Finally, we show that the cut-off frequency disappears for large in, for epsilon > 0 and increases for large in, for epsilon < 0.

Interphase anisotropy effects on lamellar eutectics: A numerical study

In directional solidification of binary eutectics, it is often observed that two-phase lamellar growth patterns grow tilted with respect to the direction z of the imposed temperature gradient. This crystallographic effect depends on the orientation of the two crystal phases alpha and beta with respect to z. Recently, an approximate theory was formulated that predicts the lamellar tilt angle as a function of the anisotropy of the free energy of the solid(alpha)-solid(beta) interphase boundary. We use two different numerical methods-phase field (PF) and dynamic boundary integral (BI)-to simulate the growth of steady periodic patterns in two dimensions as a function of the angle theta(R) between z and a reference crystallographic axis for a fixed relative orientation of alpha and beta crystals, that is, for a given anisotropy function (Wulff plot) of the interphase boundary. For Wulff plots without unstable interphase-boundary orientations, the two simulation methods are in excellent agreement with each other and confirm the general validity of the previously proposed theory. In addition, a crystallographic ``locking'' of the lamellae onto a facet plane is well reproduced in the simulations. When unstable orientations are present in the Wulff plot, it is expected that two distinct values of the tilt angle can appear for the same crystal orientation over a finite theta(R) range. This bistable behavior, which has been observed experimentally, is well reproduced by BI simulations but not by the PF model. Possible reasons for this discrepancy are discussed.

Numerical investigation of dynamic stall of an oscillating aerofoil

The flow structure around an NACA 0012 aerofoil oscillating in pitch around the quarter-chord is numerically investigated by solving the two-dimensional compressible N-S equations using a special matrix-splitting scheme. This scheme is of second-order accuracy in time and space and is computationally more efficient than the conventional flux-splitting scheme. A 'rigid' C-grid with 149 x 51 points is used for the computation of unsteady flow. The freestream Mach number varies from 0.2 to 0.6 and the Reynolds number from 5000 to 20,000. The reduced frequency equals 0.25-0.5. The basic flow structure of dynamic stall is described and the Reynolds number effect on dynamic stall is briefly discussed. The influence of the compressibility on dynamic stall is analysed in detail. Numerical results show that there is a significant influence of the compressibility on the formation and convection of the dynamic stall vortex. There is a certain influence of the Reynolds number on the flow structure. The average convection velocity of the dynamic stall vortex is approximately 0.348 times the freestream velocity.

Numerical-Value Perturbation Method with Application of Solving Convective-Diffusion Integral Equation

The convective--diffusion equation is of primary importance in such fields as fluid dynamics and heat transfer hi the numerical methods solving the convective-diffusion equation, the finite volume method can use conveniently diversified grids (structured and unstructured grids) and is suitable for very complex geometry The disadvantage of FV methods compared to the finite difference method is that FV-methods of order higher than second are more difficult to develop in three-dimensional cases. The second-order central scheme (2cs) offers a good compromise among accuracy, simplicity and efficiency, however, it will produce oscillatory solutions when the grid Reynolds numbers are large and then very fine grids are required to obtain accurate solution. The simplest first-order upwind (IUW) scheme satisfies the convective boundedness criteria, however. Its numerical diffusion is large. The power-law scheme, QMCK and second-order upwind (2UW) schemes are also often used in some commercial codes. Their numerical accurate are roughly consistent with that of ZCS. Therefore, it is meaningful to offer higher-accurate three point FV scheme. In this paper, the numerical-value perturbational method suggested by Zhi Gao is used to develop an upwind and mixed FV scheme using any higher-order interpolation and second-order integration approximations, which is called perturbational finite volume (PFV) scheme. The PFV scheme uses the least nodes similar to the standard three-point schemes, namely, the number of the nodes needed equals to unity plus the face-number of the control volume. For instanc6, in the two-dimensional (2-D) case, only four nodes for the triangle grids and five nodes for the Cartesian grids are utilized, respectively. The PFV scheme is applied on a number of 1-D problems, 2~Dand 3-D flow model equations. Comparing with other standard three-point schemes, The PFV scheme has much smaller numerical diffusion than the first-order upwind (IUW) scheme, its numerical accuracy are also higher than the second-order central scheme (2CS), the power-law scheme (PLS), the QUICK scheme and the second-order upwind(ZUW) scheme.

Numerical study on heavy rigid particle motion in a plane wake flow by spectral element method

Experimental particle dispersion patterns in a plane wake flow at a high Reynolds number have been predicted numerically by discrete vortex method (Phys. Fluids A 1992; 4:2244-2251; Int. J. Multiphase Flow 2000; 26:1583-1607). To address the particle motion at a moderate Reynolds number, spectral element method is employed to provide an instantaneous wake flow field for particle dynamics equations, which are solved to make a detail classification of the patterns in relation to the Stokes and Froude numbers. It is found that particle motion features only depend on the Stokes number at a high Froude number and depend on both numbers at a low Froude number. A ratio of the Stokes number to squared Froude number is introduced and threshold values of this parameter are evaluated that delineate the different regions of particle behavior. The parameter describes approximately the gravitational settling velocity divided by the characteristic velocity of wake flow. In order to present effects of particle density but preserve rigid sphere, hollow sphere particle dynamics in the plane wake flow is investigated. The evolution of hollow particle motion patterns for the increase of equivalent particle density corresponds to that of solid particle motion patterns for the decrease of particle size. Although the thresholds change a little, the parameter can still make a good qualitative classification of particle motion patterns as the inner diameter changes.

A Hybrid Scheme Based on Finite Element/Volume Methods for Two Immiscible Fluid Flows

We have successfully extended our implicit hybrid finite element/volume (FE/FV) solver to flows involving two immiscible fluids. The solver is based on the segregated pressure correction or projection method on staggered unstructured hybrid meshes. An intermediate velocity field is first obtained by solving the momentum equations with the matrix-free implicit cell-centered FV method. The pressure Poisson equation is solved by the node-based Galerkin FE method for an auxiliary variable. The auxiliary variable is used to update the velocity field and the pressure field. The pressure field is carefully updated by taking into account the velocity divergence field. This updating strategy can be rigorously proven to be able to eliminate the unphysical pressure boundary layer and is crucial for the correct temporal convergence rate. Our current staggered-mesh scheme is distinct from other conventional ones in that we store the velocity components at cell centers and the auxiliary variable at vertices. The fluid interface is captured by solving an advection equation for the volume fraction of one of the fluids. The same matrix-free FV method, as the one used for momentum equations, is used to solve the advection equation. We will focus on the interface sharpening strategy to minimize the smearing of the interface over time. We have developed and implemented a global mass conservation algorithm that enforces the conservation of the mass for each fluid.

Multiscale model reduction methods for deterministic and stochastic partial differential equations

Partial differential equations (PDEs) with multiscale coefficients are very difficult to solve due to the wide range of scales in the solutions. In the thesis, we propose some efficient numerical methods for both deterministic and stochastic PDEs based on the model reduction technique.

For the deterministic PDEs, the main purpose of our method is to derive an effective equation for the multiscale problem. An essential ingredient is to decompose the harmonic coordinate into a smooth part and a highly oscillatory part of which the magnitude is small. Such a decomposition plays a key role in our construction of the effective equation. We show that the solution to the effective equation is smooth, and could be resolved on a regular coarse mesh grid. Furthermore, we provide error analysis and show that the solution to the effective equation plus a correction term is close to the original multiscale solution.

For the stochastic PDEs, we propose the model reduction based data-driven stochastic method and multilevel Monte Carlo method. In the multiquery, setting and on the assumption that the ratio of the smallest scale and largest scale is not too small, we propose the multiscale data-driven stochastic method. We construct a data-driven stochastic basis and solve the coupled deterministic PDEs to obtain the solutions. For the tougher problems, we propose the multiscale multilevel Monte Carlo method. We apply the multilevel scheme to the effective equations and assemble the stiffness matrices efficiently on each coarse mesh grid. In both methods, the $\KL$ expansion plays an important role in extracting the main parts of some stochastic quantities.

For both the deterministic and stochastic PDEs, numerical results are presented to demonstrate the accuracy and robustness of the methods. We also show the computational time cost reduction in the numerical examples.

hp-Adaptive simulation and inversion of magnetotelluric measurements

xlix, 121 p.