Le plissement est généralement observé dans les structures minces ayant un comportement de type membrane. Ces structures minces ne supportent pas d'effort de flexion et sont donc sollicitées en traction. Dans cette thèse, nous avons développé une technique de réduction de modèle pour la modélisation du plissement des membranes minces. Cette technique, basée sur les séries de Fourier à double échelle, permet de déduire d'un modèle complet de membrane, un modèle réduit capable de prendre en compte les instabilités globales et locales. Les valeurs critiques de charge et longueur d'onde sont déterminées analytiquement puis numériquement. Des exemples numériques nous ont permis de valider le modèle numérique par rapport au modèle analytique. Les modèles numériques étudiés prennent en compte le modèle complet et le modèle réduit de la membrane. Le modèle complet est simulé dans Abaqus et résolu numériquement à l'aide de la méthode de la longueur d'arc et le modèle réduit est implémenté dans Matlab et résolu numériquement à l'aide de la méthode asymptotique numérique. Nous avons étudié le comportement de la membrane sous sollicitation mécanique, thermique et thermo-mécanique. Les résultats obtenus montrent que le modèle réduit est capable de se substituer au modèle complet dans la détermination des contraintes critiques et longueurs d'onde correspondantes. Le gain en temps de calcul obtenu est important, ceci grâce à la très faible densité de maillage requis par le modèle réduit. Le modèle réduit est très sensible aux conditions aux bords de la membrane et requiert d'avoir une longueur d'onde des plis quasiment constante dans la largeur de la membrane / Wrinkling is an instability phenomenon generally observed in thin structures with membrane's behavior. Those thin structures have no rigidity to flexion and are therefore used in traction. In this thesis, we developed a reduction model's technique for the modeling of wrinkling phenomenon in thin membranes. This technique, based on the double scale Fourier series, allow us to deduce from a full membrane model, a reduced membrane model that is able to take into account the global and local instability of the structure. The critical load and critical wavelength are determined analytically on one side, then numerically on the other side. Numerical exemples are conducted to validate the numerical model towards the analytical one. Numerical models studied take into account both full and reduce membrane models. The full model is simulated in Abaqus and solved numerically using the arc length method and the reduced model is implemented in Matlab and solved numerically using the asymptotic numerical method. We studied the membrane behavior under mechanical, thermal and thermo-mechanical loading. The results obtained show that the full membrane model can be replaced by the reduced one in determining critical loads and corresponding wavelengths. The gain in computation time obtained is important, due to the coarse mesh required by the reduced model. The reduced model is very sensitive to membrane's boundaries conditions and requires to have a quasi constant wavelength along the membrane width

La fabrication additive (FA) est devenue une nouvelle alternative pour la fabrication des pièces dans l'industrie. Cependant, il existe encore des limites pour ce procédé, en particulier la forme finale défavorable et les propriétés macroscopiques indésirables des pièces métalliques construites dans les systèmes de FA. La distorsion ou la fissure due à la contrainte résiduelle de ces pièces pose généralement de graves problèmes pour certains types de technologie de la FA métallique. Dans un système de FA, la qualité finale d'une pièce métallique dépend de nombreux paramètres de procédé, qui sont normalement optimisés par une série d'expériences sur des machines de FA. La simulation macroscopique dédiée au procédé de FA est une alternative potentielle pour les pièces métalliques fabriquées par la fabrication additive. Dans cette thèse, nous étudions d'abord le pré-processing de la simulation de FA par la méthode des éléments finis (FEM). Le procédé de fabrication additive est un phénomène multi-physique des champs couplés (champs thermique, mécanique et métallurgique). La simulation macroscopique est réalisée à deux niveaux différents. Au niveau de la couche, la reconstruction du modèle 3D est effectuée à partir du fichier de chemin de balayage de la machine de FA, basée sur la manipulation inverse de l'algorithme d'offsetting-clipping. Au niveau de la pièce, le modèle 3D de CAO est reconstruit dans un maillage des voxels, ce qui est pratique pour une pièce avec une géométrie complexe. Avec les températures de préchauffage différentes et les paramètres du procédé différents, la contrainte résiduelle d'une pièce est analysée. Ces simulations impliquent la technique potentielle pour réduire la contrainte résiduelle par l'optimisation des paramètres du procédé, au lieu de moyens traditionnels par augmenter la température de préchauffage. Basées sur la plateforme de simulation de FEM ci-dessus, deux simulations au niveau de ligne sont également étudiées dans cette thèse, visant à la relation entre le procédé de FA et la qualité finale de la pièce. Ces exemples démontrent la possibilité d'utiliser des simulations macroscopiques pour améliorer le contrôle de la qualité pendant le procédé de FA. Dans la première tâche, l'ensemble de données des paramètres de chauffage et la contrainte résiduelle sont générés par la simulation de FA. La corrélation entre eux est étudiée en utilisant des algorithmes de régression, tel que le réseau neuronal artificiel. Dans la deuxième tâche, un contrôleur de PID pour la boucle de rétroaction puissance-température est intégré dans la simulation de procédé de FA et l'auto-réglage de PID est numériquement étudié au lieu d'utiliser la machine de FA. Les deux tâches montrent le rôle important de la simulation de procédé macroscopique de FA, qui peut remplacer ou combiner les nombreuses expériences essai-erreur dans la fabrication additive métallique. / Additive manufacturing (AM) has become a new option for the fabrication of metallic parts in industry. However, there are still some limitations for this application, especially the unfavourable final shape and undesired macroscopic properties of metallic parts built in AM systems. The distortion or crack due to the residual stress of these parts leads usually to severe problems for some kinds of metal AM technology. In an AM system, the final quality of a metallic part depends on many process parameters, which are normally optimized by a series of experiments on AM machines. In order to reduce the considerable time consumption and financial expense of AM experiments, the numerical simulation dedicated to AM process is a prospective alternative for metallic part fabricated by additive manufacturing. Because of the multi-scale character in AM process and the complex geometrical structures of parts, most of the academic researches in AM simulation concentrated on the microscopic melting pool. Consequently, the macroscopic simulation for the AM process of a metallic part becomes a current focus in this domain. In this thesis, we first study the pre-processing of AM simulation on Finite Element Method (FEM). The process of additive manufacturing is a multi-physics problem of coupled fields (thermal, mechanical, and metallurgical fields). The macroscopic simulation is conducted in two different levels with some special pre-processing work. For the layer level, the reconstruction of 3D model is conducted from the scan path file of AM machine, based on the inverse manipulation of offsetting-clipping algorithm. For the part level, the 3D model from CAD is reconstructed into a voxel-based mesh, which is convenient for a part with complex geometry. The residual stress of a part is analysed under different preheat temperatures and different process parameters. These simulations imply the potential technique of reducing residual stress by the optimisation of process parameters, instead of the traditional way by increasing preheat temperature. Based on the FEM simulation platform above, two simulations at line level are also studied in this thesis, aiming at the relation between the AM process and part's final quality. These examples demonstrate the feasibility of using macroscopic simulations to improve the quality control during the AM process. In the first task, dataset of heating parameters and residual stress are generated by AM simulation. The correlation between them is studied by using some regression algorithm, such as artificial neural network. In the second task, a PID controller for power-temperature feedback loop is integrated into AM process simulation and the PID auto-tuning is numerically investigated instead of using AM machine. Both of the two tasks show the important role of AM macroscopic process simulation, which may replace or combine with the numerous trial and error of experiments in metal additive manufacturing.

Cette thèse s'intéresse à la problématique du contact entre matériaux hétérogènes. L'industrie (automobile, aéronautique, spatiale, ...) s'intéresse de plus en plus à ces types de matériaux. Il s'agira par exemple des alliages métalliques, des matériaux poreux, matériaux composites (composites tissés, interlocks 3D, interlocks 2D), des billes céramiques contenant des impuretés (porosités/précipités),... Dans ce manuscrit, un modèle de contact basé sur les méthodes semi-analytiques a été développé. Un algorithme de gradient conjugué est utilisé afin de résoudre rapidement le problème de contact. Le modèle permet de prendre en compte la présence d'une ou de plusieurs hétérogénéités isotropes/anisotropes dans le problème de contact. Une approche inspirée de la méthode de l'inclusion équivalente proposée par Eshelby est utilisée dans le solveur de contact pour prendre en compte l'effet de ces hétérogénéités. Les méthodes de transformées de Fourier rapides (FFT) permettent d'accélérer les calculs. Une méthode numérique a été mise en œuvre afin de prendre en compte l'interaction entre plusieurs hétérogénéités. Le massif peut être élastique ou viscoélastique. L'approche développée dans la thèse peut résoudre à la fois les problèmes d'indentation, de roulement/glissement ou de fretting en présence de matériaux élastiques hétérogènes, viscoélastiques homogènes ou hétérogènes. Les solutions sont données en termes de champs de pressions, de cisaillements et de contraintes. Dans le cas des matériaux viscoélastiques le code de calcul est capable de fournir le coefficient de frottement apparent ainsi que toutes les variables de contact aussi bien en régime permanent que transitoire. Le modèle a été validé par comparaison avec la méthode des éléments finis classiques en utilisant le logiciel commercial Abaqus v6.11. Le temps de calcul ainsi que l'espace mémoire nécessaire sont considérablement réduits par rapport à la méthode éléments finis. La parallélisation a été introduite dans le code de contact afin de réduire toujours plus le temps de calcul. Il s'agit d'un code robuste, rapide et facilement utilisable en Bureau d'Etudes. Une approche expérimentale originale a été mise en place afin de mesurer les champs de déplacements à l'interface des corps en contact. De bonnes corrélations essais/calculs ont été obtenues. Enfin quelques applications industrielles ont été présentées. Un couplage entre un code éléments finis structurel et le code semi-analytique de résolution de contact a été également réalisé. / The present PhD thesis deals with contact problems between heterogeneous materials. Nowadays heterogeneous materials are extensively used in several industrial domains (automotive, aeronautics, aerospace, ...). Heterogeneous materials involve porous materials, aluminum alloys, composites materials (woven composites, interlocks 3D, interlocks 2D), metallic or ceramics materials containing impurities (porosities/precipitates). In this work, a contact model based on semi-analytical method is proposed. A conjugate gradient algorithm is used for a fast resolution of contact equations. The model can account for one or more isotropic/anisotropic inhomogeneities. An approach taking inspiration from the Eshelby equivalent inclusion method is used in the contact solver to account for the effect of inhomogeneities. 2D and 3D Fast Fourier Transforms (FFT) are used to speed up the computation. A numerical method is implemented in order to take into account interactions between many heterogeneities. The semi-infinite space/ matrix can be either elastic or visco-elastic. The model developed in the present PhD thesis can solve indentation, rolling/sliding or fretting contact problems between heterogeneous elastic materials, homogeneous or heterogeneous visco-elactic materials. In the case of visco-elastic materials, the model permits to get the solution in terms of contact pressure distribution, subsurface stresses, apparent friction coefficient, both in the transient and then steady-state regimes. The model has been validated by performing a comparison with the results of a finite element model. The CPU time and memory necessary are greatly reduced in comparison with the classical finite element method. The model developed is fast, robust and extremely easy to use. An original experimental approach was proposed in order to measure the displacement fields at interface of two contacting bodies. A good agreement between experimental results and numerical simulations is obtained. Finally, the model is applied on some industrial applications. A coupling between a finite element model and the semi-analytical code allow to take into account the effects of structure on contact problem.

Les engrenages roues et vis sans fin sont une solution avantageuse pour transmettre le couple entre des axes perpendiculaires non concourants. Ces engrenages offrent une solution simple et efficace en terme de coût dans les applications de transmission de puissance, où un grand rapport de réduction est nécessaire, en comparaison avec les engrenages classiques à axes parallèles qui nécessitent normalement deux ou trois étapes pour obtenir les mêmes réductions avec une augmentation conséquente de complexité et du nombre de pièces. L’usure de surface est un des modes de défaillance observés dans la vie des engrenages roues et vis sans fin qui influe sur la portée de contact, les caractéristiques de transmission et le bruit résultant. La première étape de ces travaux est la mise au point d’un modèle numérique pour étudier le comportement quasi statique des engrenages roues et vis sans fin avec une roue en bronze et une vis en acier. Le modèle est basé sur la résolution des équations de compatibilité des déplacements ainsi que sur la méthode des coefficients d’influence. Les effets globaux de flexion et les effets locaux de contact ont été séparés. Les effets de contact ont été obtenus par la théorie de Boussinesq. Les coefficients de flexion sont estimés par la combinaison d’un calcul Éléments Finis et des fonctions d’interpolation, permettant d’une part de prendre en compte l’environnement de l'engrenage (la géométrie des arbres, des jantes et des voiles, l’emplacement des roulements,...) et d’autre part de réduire significativement les temps de calculs. Dans une seconde étape, une méthodologie est proposée pour modéliser l’usure de la surface de dent de la roue. Le modèle de contact quasi-statique de la répartition des charges est combiné avec un modèle d’usure d’Archard. Ce modèle suppose que la profondeur d’usure est directement proportionnelle à la pression de contact et à la distance de glissement et inversement proportionnelle à la dureté du matériau. Cette loi d’usure est modifiée pour prendre en compte l’influence des conditions de lubrification en utilisant un coefficient d’usure local, dépendant de l’épaisseur du film lubrifiant, rapportée à l’amplitude des rugosités des surfaces. L’enlèvement de matière par l’usure du flanc de la roue influe sur la répartition des pressions et donc les modifications de la géométrie des dents doivent être incluses dans la prédiction de l’usure. Le calcul des pressions de contact est ainsi mis à jour pour tenir compte des changements de géométrie. Enfin, pour valider le modèle développé des comparaisons du modèle avec des résultats expérimentaux issus de la bibliographie ont été effectuées. / Worm gears are one of the technical devices for transmitting torque between spatial crossed axes. They provide a simple and cost effective solution in power transmission applications, where a high reduction ratio is required. Comparable conventional parallel axis gearing would normally require two or three stages to achieve the same reduction, with a consequent increase in complexity and number of parts. Surface wear is one of the failure modes observed in life worm gear sets which affects the contact patterns, the other transmission characteristics and the resultant noise. The first step of this work is the development of a numerical model to study provide the quasi-static behavior of worm gears with bronze wheel and steel worm. The model is based on solving of the equation of displacement compatibility and the influence coefficient method. The global effects of bending and local effects of contact are separated. The contact effects are obtained with the theory of Boussinesq. Bending effects are estimated by the combination of one standard FEM computation and interpolation functions. These methods allow, on the one hand, to take into account the environment of the gear (shaft shape, rim, web, bearing location ...) and on the other to reduce significantly the computation time. In a second step, a methodology is proposed for predicting the wear of the wheel tooth surface. In this process, a quasi-static contact model of the load distribution is combined with Archard's wear model. This model assumes that the wear depth is directly proportional to the contact pressure and sliding distance and inversely proportional to the hardness of the material. The wear law is modified to take into account the influence of the lubrication conditions using a local wear coefficient, depending on the lubricant film thickness, relative to the amplitude of surface roughness. Removal of material by wear on the wheel flank affects the pressure distribution, therefore the changes in teeth must be included in the prediction of wear. The calculation of contact pressures must also be updated to take into account the modification of the gear flank geometry. The last step concerns the validation of the numerical. Comparisons have been carried out between the model results and experimental ones issued from the bibliography.

Actuellement, les méthodes classiques (éléments finis, endommagement, critère de rupture) pour analyser la tenue des structures sous des chargements extrêmes sont très dépendantes de la taille de maille du modèle et nécessitent à la fois un savoir-faire spé- cifique dans le domaine et des études de sensibilité au maillage. De nouvelles approches basées sur la méthode des éléments finis étendus permettent de traiter des propagations de fissure sur des structures de petites tailles et volumiques. Cependant, la propagation sur de grandes longueurs avec des modèles volumiques demande une puissance de calcul importante, souvent inaccessible dans le cadre industriel. Cette thèse a pour but de cou- pler des éléments finis de coque avec la méthode des éléments finis étendue (X-FEM). On peut ainsi diminuer la taille des modèles et gagner en temps de calcul. La fissure peut éga- lement évoluer librement dans le maillage. Après avoir fait le choix d’un élément fini de coque simple et de bonne qualité, l’objectif est de modifier cet élément afin de permettre la description d’une fissure au sein même de celui-ci. Ensuite, l’enjeu est d’adapter les critères de propagation qui existent déjà pour des modèles plans ou volumiques pour les matériaux dits "ductiles" afin de les utiliser dans le cadre d’une modélisation coque. Ces critères sont basés sur l’analyse des champs de contrainte et déformation sur un demi- disque aval à la pointe de fissure. Le calcul de la contrainte équivalente extraite de ces champs servant de seuil pour déclencher ou non la propagation est un point clef de ce travail. Cette étude se place dans le cadre de la plasticité généralisée et fait l’hypothèse d’une fissure initialement traversante dans l’épaisseur de la coque. La phase d’amorçage de la fissure n’est pas prise en compte et le défaut initial est supposé préexistant au sein de la structure. En vue de valider le couplage coque/X-FEM et le critère de propagation, des essais de fissuration sur des structures minces sont réalisés et présentés dans ce document. / In shipbuilding industry, classical methods to analyze the behavior of structures under extreme loadings are very dependent on the size of the mesh. Moreover, propagation over long lengths with volumetric models requires huge processing power, often inaccessible within this framework. In order to manage these issues and due to the geometry to be considered, a coupling between shell finite element and the extended finite element method (X-FEM) using an adapted propagation criterion is proposed. The developments are made in the fast explicit dynamic finite element code EUROPLEXUS, CEA Saclay. For shell structures involving significant thickness such as submarines, Mindlin-Reissner theory is needed to enable shear strain. Therefore, locking-free element are used to avoid the numerical issue of shear-locking that appears when the shell becomes too thin. The fracture of Mindlin-Reissner plates based on the X-FEM discrete approximation framework is studied by Dolbow and Belytschko with the MITC4. A four node shell element using the same formulation is here only enriched with a step function along the crack line to take into consideration the discontinuity of the displacement field across the crack. The calculation remains accurate without the asymptotic enrichment functions near the crack-tip, as long as the mesh is refined near the crack tip. The numerical integration issue for elements cut by the crack is solved by a partitioning strategy developed by Elguedj. Since the crack is contained in the shell for which the mid plane's position is entirely known, only one information left is needed to locate it. Therefore, a crack is represented by several line segments on the three-dimensional mesh. Only through thickness cracks are considered so far. As regards to the crack propagation, a local criteria proposed by Haboussa is used based on the calculation of mechanical equivalent quantities in the vicinity of the crack tip. The maximum of the equivalent stress tensor near the crack tip is used to decide if the crack propagates as well as its propagation direction, and the Kaninen equation gives the crack velocity.

L'étape de déliantage est une étape importante et parfois critique pour le procédé Moulage par Injection des Poudres. A cet effet, des analyses thermogravimétriques (TGA) ont été réalisées pour bien comprendre les mécanismes du déliantage thermique sous atmosphère imposée (Argon). Les méthodes de Kissinger et Ozawa ont été utilisées, en se basant sur les résultats des analyses thermogravimétriques, afin d’estimer les paramètres cinétiques nécessaires pour la simulation numérique, notre modèle se propose de décrire les phénomènes physiques liés à la dégradation du polymère, le transfert de chaleur de la déformation du composant pendant le déliantage thermique.La deuxième partie de la thèse est dédiée à la compréhension des mécanismes et du comportement du fritage des composants en tungstène sous une atmosphère d’hydrogène jusqu’à une température de 1 700°C. Des appareils expérimentaux, ont été mis en place afin de constituer une base de données physiques nécessaire pour l’identification des différents paramètres. L’identification de l’ensemble des lois de comportement a été réalisée en prenant en compte les spécificités physiques des poudres utilisées. Un modèle de comportement de type thermo élasto-viscoplastique est formulé pour représenter la loi de densification par diffusion solide, puis appliqué pour les différentes tailles de poudres de tungstène. La dernière étape consiste à valider des simulations numériques avec ABAQUS pour une meilleure détermination des densités et des retraits finaux des composants injectés. / Thermal debinding is one of the most important steps In Powder Injection Moulding process. Thermogravimetric analyses (TGA) are employed to analyze the physics and kinetics of thermal debinding behaviour under argon atmosphere. The Kissinger and Ozawa method have been used to estimate the kinetic parameters from thermogravimetric experiments. To set up the numerical simulations of thermal debinding stage using finite element method, a coupled mathematical has been developed. The basic steps of the proposed model consist to solve the following sequences of coupled problems: themal degradation of binder coupled with heat transfer and deformation phenomena by finite element method using Comsol Multiphysics software.In the second part of this thesis, sintering behaviour of tungsten powders injection moulded component, under pure hydrogen atmosphere at temperature up to 1700°C. The experimental tests are used to determine the material parameters in the parameters in the viscoplastic constitutive law, which is incorporated with the identified parameters in order to simulate the final shrinkages and densities of tungsten injection moulded components during the sintering process. Comparison between the numerical simulations results and experimental ones, in term of shrinkages and sintered densities, shows a good agreement.

De nombreux liquides subissent un changement de volume lorsqu'ils gèlent. Pour l'eau et certaines solutions aqueuses, l'expansion volumétrique au cours de la solidification peut entraîner une série de problèmes mécaniques. Dans l'industrie automobile, l'expansion de changement de phase (ECP) met en cause la sécurité des réservoirs des véhicules aux saisons froides. Une des questions les plus problématiques est l'expansion de la solution aqueuse d'urée (SAU) dans le réservoir du système SCR des véhicules diesels. Lorsque le liquide gèle, les composants intérieurs ainsi que le réservoir lui-même peuvent être endommagés dû à la pression apportée par la dilatation du liquide solidifié. Dans le centre , une méthode numérique est fortement attendue afin de prévoir la répartition de la température, des contraintes ainsi que de la déformation des composants lors d'un processus de solidification. Du fait que les informations sur la solution d'urée restent limitées, la structure de la glace cristalline ainsi que ses comportements mécaniques sont principalement passés en revue. La préférence d'orientation de croissance des grains de glace à l'interface de cristallisation met en évidence, l'hypothèse de l'ECP non-isotropique pour des problèmes de solidification. Une série de tests mécaniques a été réalisée afin d'obtenir les propriétés basiques de SAU à l'état solide à différentes températures. Une méthode « différence-volume » a été appliquée pour mesurer la variation de la densité de la SAU lors du processus de solidification. Pour la suite, des études analytiques thermiques et mécaniques sont effectuées. Pour l'aspect thermique, le problème classique de Stefan est passé en revue. Un schéma de différence-finie est proposé et il permet de calculer la position de l'interface et les profils de température pour un modèle sphérique. Pour l'aspect mécanique, un modèle sphérique similaire est établi à la base de l'ECP non-isotropique. Les solutions analytiques des contraintes et de la pression liquide sont présentées en fonction de la position de l'interface. Une méthode éléments-finis thermo-mécaniquement couplée est développée afin de simuler efficacement les contraintes thermiques, les déformations et la pression liquide dans un problème de solidification avec des relations constitutives de comportement non-linéaires. Les contraintes thermiques sont calculées en chaque point d'intégration en résolvant les équations elasto-viscoplastiques avec l'ECP non-isotropique. Le problème aux limites est résolu par la méthode de Newton-Raphson. Cette procédure est implémentée dans le package Abaqus via un UMAT. La méthode est validée d'abord pour les aspects algorithmiques par les solutions analytiques, puis pour les paramètres de comportement retenus par une série de tests expérimentaux. De plus, une étude de cas réaliste sur un réservoir de la SAU est introduite. Les avantages et les limitations de la méthode numérique lors d'une application sont évalués. / Many liquids involve a change in volume when they freeze. For water and some aqueous solutions, the volumetric expansion during solidification may invoke a series of mechanical issues. In automobile industries, the security of tanks installed in vehicles is challenged by the Phase-Change Expansion (PCE) of the freezing liquid in cold conditions. One of the most problematic issues is the expansion of Aqueous Urea Solution (AUS) in the SCR tank of diesel vehicles. As the liquid freezes, interior components may be deformed under the stress or pressure of the expanding AUS, potentially leading to failures of the storage tank. In the product center, a numerical method is of high demand to perform thermo-mechanical analysis to predict the temperature and stress distribution during a liquid solidification process in their tanks. In this work, a bibliographic study is carried out first on the basic knowledge of the ice and AUS. Due to the very limited information on urea solution in the literature, the structure and behaviors of freshwater ice are mainly reviewed. The grain orientation preference at the growth interface of polycrystalline ice provides the evidence of non-isotropic PCE for the solidification problem. A series of mechanical tests have been performed to characterize the basic properties of the solidified AUS at different temperatures. The density evolution is measured using a volume-difference method. Then, both thermal and mechanical analytical studies are performed. The classical thermal Stefan problem is reviewed and a finite-difference scheme is proposed to calculate the interface position and temperature profiles of a spherical solidification model. Mechanically, a similar spherical model is established based on the non-isotropic PCE phenomenon of ice growth. The solutions of stress distribution and liquid pressure evolution are given as a function of the solidification interface position. Finally, an efficient thermo-mechanical FEM is proposed to evaluate the thermal stress, strain, displacement and pressure in solidification problems with highly nonlinear relations. Three particular methods for treating the liquid phase with fixed-grid approaches are introduced. The thermal stress is computed at each integration point by integrating the elasto-viscoplastic constitutive equations with non-isotropic PCE. Then, the boundary value problem is solved using the full Newton-Raphson method. This procedure is implemented into the FE package Abaqus via a UMAT subroutine. The numerical model is validated first for the algorithmic aspect by the analytical solutions, and then for the parametric calibration by a series of benchmark tests. In the end, a realistic study case on a real-size AUS storage tank is introduced. Advantages and limitations of the numerical method in the application are evaluated.