3:00 PM - 3:15 PM

# [SCG74-18] Numerical modeling of fracture of porous Maxwell fluid by phase-field method

Keywords:Explosive eruption, Fragmentation, Phase-field method

Brittle fragmentation of vesicular magma is a key process in explosive eruption. Recent estimation on the decompression time in real explosive events indicates that the style of fragmentation is to be “brittle-like fragmentation” (Kameda et al. JVGR 2013), which is defined as the solid-like fracture of the material whose bulk rheological properties is close to fluid state. We also found the fact that the spatial inhomogeneity of bubble distribution is a major source of crack development that may lead to brittle-like fragmentation (Kameda et al. in preparation).

In order to model our findings by numerical approach, we propose a continuum description of fracture in viscoelastic fluid using a phase-field method (Spatschek et al. Phil. Mag. 2011). The phase-field modes with sharp interfaces consist in incorporating a continuous field variable so-called “order parameter,” by which the magma and the bubbles (or the cracks) are distinguished. In this study, the evolution of the order parameter is described using the Allen-Cahn equation consisting of local elastic energy, surface energy and additional numerical parameters. Evolution of crack is driven by local elastic energy. A commercial multiphysics solver (COMSOL) is used as the platform of numerical simulation. The evolution of the order parameter is solved using PDE solver of COMSOL in which weak formulation of the governing equation is described by ourselves. The evolution of strain-stress field in the system is calculated using conventional finite element analysis (FEM) prepared in the solver, in which the rheology of unbroken magma is assumed to be a linear Maxwell fluid.

We consider a sphere consisting of Maxwell fluid and gas bubbles (Fig. 1). In order to reduce the computational cost, one-eighth portion of the sphere is used as the computational domain. A large bubble is placed at the center of the sphere. Another satellite bubble is placed beneath the large bubble. We tested two arrangements about the location of the satellite bubble. In one case, it is placed at the position equal distance away from three planes of symmetry. In another case, the position is biased toward one of the symmetry plane. Physical properties of the materials are determined according to our previous laboratory experiment (Kameda et al. 2013). Isotropic (spherical) negative pressure (gradually increases with time) is loaded outer surface of the computational domain as a decompression.

Time evolution of Maxwell fluid/gas interface is displayed in Fig. 2. Figure 2 indicates that the arrangement of the bubbles remarkably affects the evolution of the crack. The process is divided into four parts: First, a tip is formed on the satellite bubble toward the large bubble. Second, the tip reaches the surface of the large bubble. Third, only in the case where the position of the satellite bubble is biased, a plane crack beneath the large bubble is rapidly developed toward the plane of symmetry. Finally, the plane crack propagates toward outer boundary, then a sharp crack is opened.

We found that stress concentration due to mutual interaction of neighboring bubbles is a dominant factor to evolution of cracks in porous material. Differential stress around a solitary bubble by decompression is inversely proportional to the cube of relative distance in radial direction normalized by its radius (Zhang Nature 1999). In the field of multiple bubbles, the local differential stress is close to the sum of stresses induced by neighboring bubbles. This means that the maximum differential stress is found at the surface of small bubbles close to the large bubble. Evolution of the crack shown in Fig. 2 is explained well with distribution local stress concentration.

In conclusion, this simulation demonstrates that the evolution of crack is very sensitive to the arrangement of bubbles. Therefore, we should incorporate the effect of bubble arrangement into the criteria for the onset of fragmentation of vesicular magma.

In order to model our findings by numerical approach, we propose a continuum description of fracture in viscoelastic fluid using a phase-field method (Spatschek et al. Phil. Mag. 2011). The phase-field modes with sharp interfaces consist in incorporating a continuous field variable so-called “order parameter,” by which the magma and the bubbles (or the cracks) are distinguished. In this study, the evolution of the order parameter is described using the Allen-Cahn equation consisting of local elastic energy, surface energy and additional numerical parameters. Evolution of crack is driven by local elastic energy. A commercial multiphysics solver (COMSOL) is used as the platform of numerical simulation. The evolution of the order parameter is solved using PDE solver of COMSOL in which weak formulation of the governing equation is described by ourselves. The evolution of strain-stress field in the system is calculated using conventional finite element analysis (FEM) prepared in the solver, in which the rheology of unbroken magma is assumed to be a linear Maxwell fluid.

We consider a sphere consisting of Maxwell fluid and gas bubbles (Fig. 1). In order to reduce the computational cost, one-eighth portion of the sphere is used as the computational domain. A large bubble is placed at the center of the sphere. Another satellite bubble is placed beneath the large bubble. We tested two arrangements about the location of the satellite bubble. In one case, it is placed at the position equal distance away from three planes of symmetry. In another case, the position is biased toward one of the symmetry plane. Physical properties of the materials are determined according to our previous laboratory experiment (Kameda et al. 2013). Isotropic (spherical) negative pressure (gradually increases with time) is loaded outer surface of the computational domain as a decompression.

Time evolution of Maxwell fluid/gas interface is displayed in Fig. 2. Figure 2 indicates that the arrangement of the bubbles remarkably affects the evolution of the crack. The process is divided into four parts: First, a tip is formed on the satellite bubble toward the large bubble. Second, the tip reaches the surface of the large bubble. Third, only in the case where the position of the satellite bubble is biased, a plane crack beneath the large bubble is rapidly developed toward the plane of symmetry. Finally, the plane crack propagates toward outer boundary, then a sharp crack is opened.

We found that stress concentration due to mutual interaction of neighboring bubbles is a dominant factor to evolution of cracks in porous material. Differential stress around a solitary bubble by decompression is inversely proportional to the cube of relative distance in radial direction normalized by its radius (Zhang Nature 1999). In the field of multiple bubbles, the local differential stress is close to the sum of stresses induced by neighboring bubbles. This means that the maximum differential stress is found at the surface of small bubbles close to the large bubble. Evolution of the crack shown in Fig. 2 is explained well with distribution local stress concentration.

In conclusion, this simulation demonstrates that the evolution of crack is very sensitive to the arrangement of bubbles. Therefore, we should incorporate the effect of bubble arrangement into the criteria for the onset of fragmentation of vesicular magma.