diff --git a/modules/xfem/include/postprocessors/OxideLayerThickness.h b/modules/xfem/include/postprocessors/OxideLayerThickness.h new file mode 100644 index 000000000000..d22ce5510134 --- /dev/null +++ b/modules/xfem/include/postprocessors/OxideLayerThickness.h @@ -0,0 +1,38 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#ifndef OXIDELAYERTHICKNESS_H +#define OXIDELAYERTHICKNESS_H + +#include "ElementPostprocessor.h" +#include "MovingLineSegmentCutSetUserObject.h" +// Forward Declarations +class OxideLayerThickness; + +template <> +InputParameters validParams(); + +class OxideLayerThickness : public ElementPostprocessor +{ +public: + OxideLayerThickness(const InputParameters & parameters); + + void initialize() override {} + void execute() override{}; + void finalize() override{}; + void threadJoin(const UserObject & user_object) override{}; + + virtual Real getValue() override; + +protected: + const MovingLineSegmentCutSetUserObject * const _moving_line_segments; + const unsigned int _cut_data_index; +}; + +#endif // OXIDELAYERTHICKNESS_H diff --git a/modules/xfem/include/userobjects/MovingLineSegmentCutSetUserObject.h b/modules/xfem/include/userobjects/MovingLineSegmentCutSetUserObject.h index 508c87bdd1bb..d00f1131af15 100644 --- a/modules/xfem/include/userobjects/MovingLineSegmentCutSetUserObject.h +++ b/modules/xfem/include/userobjects/MovingLineSegmentCutSetUserObject.h @@ -37,6 +37,8 @@ class MovingLineSegmentCutSetUserObject : public LineSegmentCutSetUserObject virtual Real cutFraction(unsigned int cut_num, Real time) const override; + virtual Real getCutDataXCoord(unsigned int i) const { return _cut_data[i * 6]; }; + /// Pointer to XFEMMovingInterfaceVelocityBase object const XFEMMovingInterfaceVelocityBase * _interface_velocity; }; diff --git a/modules/xfem/src/postprocessors/OxideLayerThickness.C b/modules/xfem/src/postprocessors/OxideLayerThickness.C new file mode 100644 index 000000000000..ee5a1b6611ae --- /dev/null +++ b/modules/xfem/src/postprocessors/OxideLayerThickness.C @@ -0,0 +1,37 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "OxideLayerThickness.h" + +registerMooseObject("XFEMApp", OxideLayerThickness); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addRequiredParam("moving_line_segments", + "The MovingLineSegmentCutSetUserObject user object name"); + params.addParam("cut_data_index", 0, "The index of the cut data"); + return params; +} + +OxideLayerThickness::OxideLayerThickness(const InputParameters & parameters) + : ElementPostprocessor(parameters), + _moving_line_segments( + &getUserObject("moving_line_segments")), + _cut_data_index(getParam("cut_data_index")) +{ +} + +Real +OxideLayerThickness::getValue() +{ + return 6e-4 - _moving_line_segments->getCutDataXCoord(_cut_data_index); +} diff --git a/modules/xfem/src/userobjects/XFEMC4OxideVelocity.C b/modules/xfem/src/userobjects/XFEMC4OxideVelocity.C index abf0dc9a36fd..2285e0e9ac97 100644 --- a/modules/xfem/src/userobjects/XFEMC4OxideVelocity.C +++ b/modules/xfem/src/userobjects/XFEMC4OxideVelocity.C @@ -49,7 +49,7 @@ XFEMC4OxideVelocity::computeMovingInterfaceVelocity(unsigned int point_id) const Real delta = std::abs(xt - _x0); -// std::cout << "delta: " << delta << std::endl; + // std::cout << "delta: " << delta << std::endl; // Current implementation only supports the case that the interface is moving horizontally // return std::abs((_diffusivity_at_positive_level_set * grad_positive(0) - @@ -72,8 +72,8 @@ XFEMC4OxideVelocity::computeMovingInterfaceVelocity(unsigned int point_id) const const Real con_e_ox_m = 2 * con_v_ox_m; const Real con_o_ox_m(5.24e28); const Real con_o_m_ox(1.7078e28); - -// use fixed temperature at 1200C until we add a temperature diffusion kernel + + // use fixed temperature at 1200C until we add a temperature diffusion kernel const Real temperature(1473); const Real mobil_v = 4 * pow(migr_jp_l, 2) * migr_jp_f * 2 / (Kb * temperature) * @@ -87,10 +87,13 @@ XFEMC4OxideVelocity::computeMovingInterfaceVelocity(unsigned int point_id) const const Real eta = (-B - sqrt(pow(B, 2) - 4 * A * C)) / (2 * A); const Real potential = Kb * temperature * log(eta); - const Real J_v = mobil_v * potential * (con_v_ox_w - con_v_ox_m * pow(eta, 2)) / - (1 - pow(eta, 2)) / delta; + const Real J_v = + mobil_v * potential * (con_v_ox_w - con_v_ox_m * pow(eta, 2)) / (1 - pow(eta, 2)) / delta; + // const Real J_o = zircaloy_density * Na / (zirconium_atmass * 1e-3) * + // _diffusivity_at_positive_level_set * (grad_positive(0) + grad_positive(1) + + // grad_positive(2)) / 3 * 1e-6; const Real J_o = zircaloy_density * Na / (zirconium_atmass * 1e-3) * - _diffusivity_at_positive_level_set * (grad_positive(0) + grad_positive(1) + grad_positive(2)) / 3 * 1e-6; + _diffusivity_at_positive_level_set * (grad_positive(0)) * 1e-6; if (delta == 0) return sqrt(0.01126 * exp(-35890 / (1.987 * temperature)) / (2 * _t)) * (-1e-2); diff --git a/modules/xfem/test/tests/moving_interface/moving_oxide_C4_weak.i b/modules/xfem/test/tests/moving_interface/moving_oxide_C4_weak.i new file mode 100644 index 000000000000..74ce341d3047 --- /dev/null +++ b/modules/xfem/test/tests/moving_interface/moving_oxide_C4_weak.i @@ -0,0 +1,177 @@ +# test for an oxide growing on top of a ziroconium nuclear fuel cladding +# using the C4 model to compute the growth rate + +[GlobalParams] + order = FIRST + family = LAGRANGE +[] + +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 21 + ny = 20 + xmin = 0 + xmax = 6e-4 + ymin = 0 + ymax = 1e-3 + elem_type = QUAD4 +[] + +[XFEM] + qrule = volfrac + output_cut_plane = true +[] + +[UserObjects] + [./velocity] + type = XFEMC4OxideVelocity + diffusivity_at_positive_level_set = 0.2978e-5 #1e-5*0.2978 + diffusivity_at_negative_level_set = 0.6667e-11 #1e-11*0.6667 + equilibrium_concentration_jump = 0.3689 + value_at_interface_uo = value_uo + x0 = 5.9e-4 + [../] + [./value_uo] + type = PointValueAtXFEMInterface + variable = 'u' + geometric_cut_userobject = 'moving_line_segments' + execute_on = 'nonlinear' + level_set_var = ls + [../] + [./moving_line_segments] + type = MovingLineSegmentCutSetUserObject + cut_data = '5.9e-4 0 5.9e-4 1e-3 0 0' + heal_always = true + interface_velocity = velocity + [../] +[] + +[Variables] + [./u] + [../] +[] + +[ICs] + [./ic_u] + type = FunctionIC + variable = u + function = 'if(x<5.9e-4, 0.03, 0.6667)' + [../] +[] + +[AuxVariables] + [./ls] + order = FIRST + family = LAGRANGE + [../] +[] + +# Need to use XFEMTwoSideDirichlet that has been removed +[Constraints] + [./u_constraint] + type = XFEMTwoSideDirichlet + geometric_cut_userobject = 'moving_line_segments' + use_displaced_mesh = false + variable = u + value_at_positive_level_set_interface = 1 #0.2978 + value_at_negative_level_set_interface = 1 #0.6667 +# value = 0.6667 + alpha = 1e5 + [../] +[] + +[Kernels] + [./diff] + type = ConcentrationDiffusion + variable = u + diffusion_coefficient_name = 'diffusion_coefficient' + [../] + [./time] + type = TimeDerivative + variable = u + [../] +[] + +[AuxKernels] + [./ls] + type = LineSegmentLevelSetAux + line_segment_cut_set_user_object = 'moving_line_segments' + variable = ls + [../] +[] + +[Materials] + [./diffusivity_A] + type = GenericConstantMaterial + prop_names = A_diffusion_coefficient + prop_values = 1e-5 + [../] + [./diffusivity_B] + type = GenericConstantMaterial + prop_names = B_diffusion_coefficient + prop_values = 1e-11 + [../] + [./diff_combined] + type = LevelSetBiMaterialReal + levelset_positive_base = 'A' + levelset_negative_base = 'B' + level_set_var = ls + prop_name = diffusion_coefficient + [../] +[] + +[BCs] +# Define boundary conditions + [./left_u] + type = DirichletBC + variable = u + value = 0.0450 #0.03/0.6667 + boundary = left + [../] + + [./right_u] + type = DirichletBC + variable = u + value = 2.2388 #0.6667/0.2978 + boundary = right + [../] +[] + +[Postprocessors] + [./cut_data_x] + type = OxideLayerThickness + moving_line_segments = moving_line_segments + execute_on = TIMESTEP_END + [../] +[] + +[Executioner] + type = Transient + solve_type = 'PJFNK' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + line_search = 'none' + + l_tol = 1e-3 + nl_max_its = 15 + nl_rel_tol = 1e-7 + nl_abs_tol = 1e-7 + + start_time = 0.0 + dt = 10 + num_steps = 150 + max_xfem_update = 1 +[] + + +[Outputs] + execute_on = timestep_end + exodus = true + [./console] + type = Console + output_linear = true + [../] + csv = true + perf_graph = true +[] diff --git a/modules/xfem/test/tests/moving_interface/moving_oxide_CP_weak.i b/modules/xfem/test/tests/moving_interface/moving_oxide_CP_weak.i new file mode 100644 index 000000000000..0bdbbd889aa8 --- /dev/null +++ b/modules/xfem/test/tests/moving_interface/moving_oxide_CP_weak.i @@ -0,0 +1,176 @@ +# test for an oxide growing on top of a ziroconium nuclear fuel cladding +# using the Cathcart-Pawel model to compute the growth rate + +[GlobalParams] + order = FIRST + family = LAGRANGE +[] + +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 21 + ny = 20 + xmin = 0 + xmax = 6e-4 + ymin = 0 + ymax = 1e-3 + elem_type = QUAD4 +[] + +[XFEM] + qrule = volfrac + output_cut_plane = true +[] + +[UserObjects] + [./velocity] + type = XFEMCathcartPawelOxideVelocity + diffusivity_at_positive_level_set = 1e-5 + diffusivity_at_negative_level_set = 1e-11 + equilibrium_concentration_jump = 0.3689 + value_at_interface_uo = value_uo + [../] + [./value_uo] + type = PointValueAtXFEMInterface + variable = 'u' + geometric_cut_userobject = 'moving_line_segments' + execute_on = 'nonlinear' + level_set_var = ls + [../] + [./moving_line_segments] + type = MovingLineSegmentCutSetUserObject + cut_data = '5.9e-4 0 5.9e-4 1e-3 0 0' + heal_always = true + interface_velocity = velocity + [../] +[] + +[Variables] + [./u] + [../] +[] + +[ICs] + [./ic_u] + type = FunctionIC + variable = u + function = 'if(x<5.9e-4, 0.03, 0.6667)' + [../] +[] + +[AuxVariables] + [./ls] + order = FIRST + family = LAGRANGE + [../] +[] + +# Need to use XFEMTwoSideDirichlet that has been removed +[Constraints] + [./u_constraint] + type = XFEMTwoSideDirichlet + geometric_cut_userobject = 'moving_line_segments' + use_displaced_mesh = false + variable = u + value_at_positive_level_set_interface = 1.0 #0.2978 + value_at_negative_level_set_interface = 1.0 #0.6667 +# value = 0.6667 + alpha = 1e5 + [../] +[] + +[Kernels] + [./diff] + type = ConcentrationDiffusion + variable = u + diffusion_coefficient_name = 'diffusion_coefficient' + [../] + [./time] + type = TimeDerivative + variable = u + [../] +[] + +[AuxKernels] + [./ls] + type = LineSegmentLevelSetAux + line_segment_cut_set_user_object = 'moving_line_segments' + variable = ls + [../] +[] + +[Materials] + [./diffusivity_A] + type = GenericConstantMaterial + prop_names = A_diffusion_coefficient + prop_values = 1e-5 + [../] + [./diffusivity_B] + type = GenericConstantMaterial + prop_names = B_diffusion_coefficient + prop_values = 1e-11 + [../] + [./diff_combined] + type = LevelSetBiMaterialReal + levelset_positive_base = 'A' + levelset_negative_base = 'B' + level_set_var = ls + prop_name = diffusion_coefficient + [../] +[] + +[BCs] +# Define boundary conditions + [./left_u] + type = DirichletBC + variable = u + value = 0.0450 #0.03/0.6667 + boundary = left + [../] + + [./right_u] + type = DirichletBC + variable = u + value = 2.2388 #0.6667/0.2978 + boundary = right + [../] +[] + +[Postprocessors] + [./cut_data_x] + type = OxideLayerThickness + moving_line_segments = moving_line_segments + execute_on = TIMESTEP_END + [../] +[] + +[Executioner] + type = Transient + solve_type = 'PJFNK' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + line_search = 'none' + + l_tol = 1e-3 + nl_max_its = 15 + nl_rel_tol = 1e-7 + nl_abs_tol = 1e-7 + + start_time = 0.0 + dt = 10 + num_steps = 150 + max_xfem_update = 1 +[] + + +[Outputs] + execute_on = timestep_end + exodus = true + [./console] + type = Console + output_linear = true + [../] + csv = true + perf_graph = true +[] diff --git a/modules/xfem/test/tests/moving_interface/moving_oxide_CP_weak_mechanics.i b/modules/xfem/test/tests/moving_interface/moving_oxide_CP_weak_mechanics.i new file mode 100644 index 000000000000..3b0f409b9933 --- /dev/null +++ b/modules/xfem/test/tests/moving_interface/moving_oxide_CP_weak_mechanics.i @@ -0,0 +1,383 @@ +# test for an oxide growing on top of a ziroconium nuclear fuel cladding +# using the Cathcart-Pawel model to compute the growth rate + +[GlobalParams] + order = FIRST + family = LAGRANGE +[] + +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 21 + ny = 20 + xmin = 0 + xmax = 6e-4 + ymin = 0 + ymax = 1e-3 + elem_type = QUAD4 +[] + +[MeshModifiers] + [./left_bottom] + type = AddExtraNodeset + new_boundary = 'left_bottom' + coord = '0.0 0.0' + [../] +[] + +[XFEM] + qrule = volfrac + output_cut_plane = true +[] + +[UserObjects] + [./velocity] + type = XFEMCathcartPawelOxideVelocity + diffusivity_at_positive_level_set = 1e-5 + diffusivity_at_negative_level_set = 1e-11 + equilibrium_concentration_jump = 0.3689 + value_at_interface_uo = value_uo + [../] + [./value_uo] + type = PointValueAtXFEMInterface + variable = 'u' + geometric_cut_userobject = 'moving_line_segments' + execute_on = 'nonlinear' + level_set_var = ls + [../] + [./moving_line_segments] + type = MovingLineSegmentCutSetUserObject + cut_data = '5.9e-4 0 5.9e-4 1e-3 0 0' + heal_always = true + interface_velocity = velocity + [../] +[] + +[Variables] + [./u] + [../] + [./disp_x] + [../] + [./disp_y] + [../] +[] + +[ICs] + [./ic_u] + type = FunctionIC + variable = u + function = 'if(x<5.9e-4, 0.03, 0.6667)' + [../] +[] + +[AuxVariables] + [./temp] + initial_condition = 1200 + [../] + [./ls] + order = FIRST + family = LAGRANGE + [../] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./a_strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./a_strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./a_strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./b_strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./b_strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./b_strain_xy] + order = CONSTANT + family = MONOMIAL + [../] +[] + +# Need to use XFEMTwoSideDirichlet that has been removed +[Constraints] + [./u_constraint] + type = XFEMTwoSideDirichlet + geometric_cut_userobject = 'moving_line_segments' + use_displaced_mesh = false + variable = u + value_at_positive_level_set_interface = 1.0 #0.2978 + value_at_negative_level_set_interface = 1.0 #0.6667 +# value = 0.6667 + alpha = 1e5 + [../] + [./dispx_constraint] + type = XFEMSingleVariableConstraint + use_displaced_mesh = false + variable = disp_x + alpha = 1e8 + geometric_cut_userobject = 'moving_line_segments' + [../] + [./dispy_constraint] + type = XFEMSingleVariableConstraint + use_displaced_mesh = false + variable = disp_y + alpha = 1e8 + geometric_cut_userobject = 'moving_line_segments' + [../] +[] + +[Kernels] + [./diff] + type = ConcentrationDiffusion + variable = u + diffusion_coefficient_name = 'diffusion_coefficient' + [../] + [./time] + type = TimeDerivative + variable = u + [../] + [./TensorMechanics] + displacements = 'disp_x disp_y' + use_displaced_mesh = false + [../] +[] + +[AuxKernels] + [./ls] + type = LineSegmentLevelSetAux + line_segment_cut_set_user_object = 'moving_line_segments' + variable = ls + [../] + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + index_i = 0 + index_j = 0 + variable = stress_xx + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + index_i = 1 + index_j = 1 + variable = stress_yy + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + index_i = 0 + index_j = 1 + variable = stress_xy + [../] + [./a_strain_xx] + type = RankTwoAux + rank_two_tensor = A_total_strain + index_i = 0 + index_j = 0 + variable = a_strain_xx + [../] + [./a_strain_yy] + type = RankTwoAux + rank_two_tensor = A_total_strain + index_i = 1 + index_j = 1 + variable = a_strain_yy + [../] + [./a_strain_xy] + type = RankTwoAux + rank_two_tensor = A_total_strain + index_i = 0 + index_j = 1 + variable = a_strain_xy + [../] + [./b_strain_xx] + type = RankTwoAux + rank_two_tensor = B_total_strain + index_i = 0 + index_j = 0 + variable = b_strain_xx + [../] + [./b_strain_yy] + type = RankTwoAux + rank_two_tensor = B_total_strain + index_i = 1 + index_j = 1 + variable = b_strain_yy + [../] + [./b_strain_xy] + type = RankTwoAux + rank_two_tensor = B_total_strain + index_i = 0 + index_j = 1 + variable = b_strain_xy + [../] +[] + +[Materials] + [./diffusivity_A] + type = GenericConstantMaterial + prop_names = A_diffusion_coefficient + prop_values = 1e-5 + [../] + [./diffusivity_B] + type = GenericConstantMaterial + prop_names = B_diffusion_coefficient + prop_values = 1e-11 + [../] + [./diff_combined] + type = LevelSetBiMaterialReal + levelset_positive_base = 'A' + levelset_negative_base = 'B' + level_set_var = ls + prop_name = diffusion_coefficient + [../] + [./elasticity_tensor_A] + type = ComputeIsotropicElasticityTensor + base_name = A + youngs_modulus = 1e3 + poissons_ratio = 0.3 + [../] + [./strain_A] + type = ComputeSmallStrain + base_name = A + displacements = 'disp_x disp_y' + eigenstrain_names = 'eigenstrain_A' + [../] + [./thermal_expansion_strain_A] + type = ComputeThermalExpansionEigenstrain + block = 0 + stress_free_temperature = 298 + thermal_expansion_coeff = 1.0e-5 + temperature = temp + eigenstrain_name = A_eigenstrain_A + [../] + [./stress_A] + type = ComputeLinearElasticStress + base_name = A + [../] + [./elasticity_tensor_B] + type = ComputeIsotropicElasticityTensor + base_name = B + youngs_modulus = 1e2 + poissons_ratio = 0.3 + [../] + [./strain_B] + type = ComputeSmallStrain + base_name = B + displacements = 'disp_x disp_y' + eigenstrain_names = 'eigenstrain_B' + [../] + [./thermal_expansion_strain_B] + type = ComputeThermalExpansionEigenstrain + block = 0 + stress_free_temperature = 298 + thermal_expansion_coeff = 1.0e-6 + temperature = temp + eigenstrain_name = eigenstrain_B + [../] + [./stress_B] + type = ComputeLinearElasticStress + base_name = B + [../] + [./combined_stress] + type = LevelSetBiMaterialRankTwo + levelset_positive_base = 'A' + levelset_negative_base = 'B' + level_set_var = ls + prop_name = stress + [../] + [./combined_dstressdstrain] + type = LevelSetBiMaterialRankFour + levelset_positive_base = 'A' + levelset_negative_base = 'B' + level_set_var = ls + prop_name = Jacobian_mult + [../] +[] + +[BCs] +# Define boundary conditions + [./left_u] + type = DirichletBC + variable = u + value = 0.0450 #0.03/0.6667 + boundary = left + [../] + + [./right_u] + type = DirichletBC + variable = u + value = 2.2388 #0.6667/0.2978 + boundary = right + [../] + + [./bottomy] + type = PresetBC + boundary = left_bottom + variable = disp_y + value = 0.0 + [../] + [./left_x] + type = FunctionPresetBC + boundary = left + variable = disp_x + function = 0.0 + [../] +[] + +[Postprocessors] + [./cut_data_x] + type = OxideLayerThickness + moving_line_segments = moving_line_segments + execute_on = TIMESTEP_END + [../] +[] + +[Executioner] + type = Transient + solve_type = 'PJFNK' + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + line_search = 'none' + + l_tol = 1e-3 + nl_max_its = 15 + nl_rel_tol = 1e-7 + nl_abs_tol = 1e-7 + + start_time = 0.0 + dt = 10 + num_steps = 150 + max_xfem_update = 1 +[] + + +[Outputs] + execute_on = timestep_end + exodus = true + [./console] + type = Console + output_linear = true + [../] + csv = true + perf_graph = true +[] diff --git a/modules/xfem/test/tests/moving_interface/plot.py b/modules/xfem/test/tests/moving_interface/plot.py new file mode 100644 index 000000000000..b900cd7f8673 --- /dev/null +++ b/modules/xfem/test/tests/moving_interface/plot.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python + +import pandas +import matplotlib.pyplot as plt +import numpy as np +import math + +ax = plt.gca() + +CP = pandas.read_csv('moving_oxide_CP_weak_mechanics_out.csv') +time_cp = CP['time'] +thickness_cp = CP['cut_data_x'] + +C4 = pandas.read_csv('moving_oxide_C4_weak_out.csv') +time_c4 = C4['time'] +thickness_c4 = C4['cut_data_x'] + +plt.plot(time_cp,thickness_cp,label='CP (BISON XFEM)') +plt.plot(time_c4,thickness_c4,label='C4 (BISON XFEM)') + +#plt.title("Plastic strain") +plt.xlabel("time[s]") +plt.ylabel("oxide thickness [$\mu$ m]") + +ax.legend() +#ax.set_xlim(left=0) +#ax.set_ylim(bottom=0,top = 0.5) +plt.savefig('C4_and_CP_model.pdf')