#!/usr/bin/env python3 # -*- coding: UTF8 -*- # # Copyright (C) 2018-2020 Konstantinos Poulios, Yves Renard. # # This file is a part of GetFEM++ # # GetFEM++ is free software; you can redistribute it and/or modify it # under the terms of the GNU Lesser General Public License as published # by the Free Software Foundation; either version 2.1 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public # License for more details. # You should have received a copy of the GNU Lesser General Public License # along with this program; if not, write to the Free Software Foundation, # Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. # ############################################################################ import getfem as gf import numpy as np gf.util_trace_level(1) # Input data NX = 80 # number of elements in horizontal direction NY = 60 # number of elements in vertical direction LX = 1. # [mm] Length LY = 1. # [mm] Height E = 210e3 # [N/mm^2] nu = 0.3 rho = 8e-9 # [ton/mm^3] Gc0 = 2.7e0 # [N/mm] chi = 0.89 # [-] hydrogen damage coefficient Dgb0 = 30e6 # [N*mm/mol] Gibbs free energy jump C0 = 0.5 # [wt ppm] initial hydrogen concentration VH = 2000. # [mm^3/mol] molar volume of hydrogen in solution R = 8.314e3 # [N*mm/(K*mol)] universal gas constant T = 300. # [K] absolute temperature ll = 0.015 # [mm] length scale strain_rate = 2e-4 init_time_step = 0.4 steps = 250 press_fem_order = 1 # pressure finite element order disp_fem_order = 2 # displacements finite element order damage_fem_order = 2 # phase field finite element order C_fem_order = 1 # hydrogen concentration finite element order # mesh generation NX = 2*((NX+1)//2) NY = 2*((NY+1)//2) XX1 = np.linspace(-1., 0., NX//4+1) XX2 = np.linspace(0., 1., (NX-NX//4)+1)[1:] XX = np.concatenate((LX/2*(0.2*XX1+0.8*np.sign(XX1)*np.power(np.abs(XX1),1.5)), LX/2*(1.0*XX2+0.0*np.sign(XX2)*np.power(np.abs(XX2),1.5)))) YY = np.linspace(-1., 1., NY+1) YY = LY/2*(0.2*YY+0.8*np.sign(YY)*np.power(np.abs(YY),1.5)) mesh = gf.Mesh("cartesian", XX, YY) # mesh regions T_RG = 6 TB_RG = 7 EXT_RG = 8 bottom_faces = mesh.outer_faces_in_box([-LX/2-1e-5,-LY/2-1e-5],[LX/2+1e-5,-LY/2+1e-5]) top_faces = mesh.outer_faces_in_box([-LX/2-1e-5,LY/2-1e-5],[LX/2+1e-5,LY/2+1e-5]) mesh.set_region(T_RG, top_faces) mesh.set_region(TB_RG, bottom_faces) mesh.extend_region(TB_RG, top_faces) mesh.set_region(EXT_RG, mesh.outer_faces()) # FEM mfp = gf.MeshFem(mesh, 1) mfp.set_classical_fem(press_fem_order) mfu = gf.MeshFem(mesh, 2) mfu.set_classical_fem(disp_fem_order) mfdir = mfu mfd = gf.MeshFem(mesh, 1) mfd.set_classical_fem(damage_fem_order) mfC = gf.MeshFem(mesh, 1) mfC.set_classical_fem(C_fem_order) mfout = gf.MeshFem(mesh) mfout.set_classical_discontinuous_fem(2) # Integration methods mim4 = gf.MeshIm(mesh, 3) mim9 = gf.MeshIm(mesh, 5) mimd4 = gf.MeshImData(mim4, -1) mimd9 = gf.MeshImData(mim9, -1) # Model md = gf.Model("real") md.add_fem_variable("u", mfu) # displacements field md.add_fem_variable("d", mfd) # fracture phase field md.add_fem_variable("p", mfp) # hydrostatic pressure field md.add_fem_variable("C", mfC) # hydrogen concentration field md.set_variable("C", C0*np.ones(mfC.nbdof())) md.add_fem_data("u_prev", mfu) md.add_fem_data("v_prev", mfu) md.add_im_data("psi0_max", mimd9) md.add_initialized_data("Dt", init_time_step) md.add_initialized_data("kappa", E/(3.*(1.-2.*nu))) md.add_initialized_data("mu", E/(2*(1+nu))) md.add_initialized_data("rho", rho) md.add_initialized_data("Gc0", Gc0) md.add_initialized_data("l", ll) md.add_initialized_data("C0", C0) md.add_initialized_data("chi", chi) md.add_initialized_data("c1", np.exp(-Dgb0/(R*T))/5.5e-5) md.add_initialized_data("c2", VH/(R*T)) # initial crack md.set_variable("psi0_max", md.interpolation("1e4*Gc0*min(1,max(0,1-(abs(X(2))+pos_part(X(1)))/(l/2)))", mimd9)) md.add_linear_term(mim9, "rho/Dt*((u-u_prev)/Dt-v_prev).Test_u") md.add_macro("degradation", "sqr(1-d)+1e-7") md.add_macro("deveps", "Sym(Grad(u))-Div(u)/3*Id(2)") md.add_macro("psi0", "(0.5*kappa*sqr(Div(u))+mu*Norm_sqr(deveps))") md.add_macro("Gc", "(1-chi*C/(C+c1))*Gc0") md.add_nonlinear_term(mim9, "degradation*(kappa*Div(u)*Id(2)+2*mu*deveps):Grad(Test_u)") md.add_nonlinear_term(mim9, "(-2*(1-d)*max(psi0_max,psi0)*Test_d" "+Gc*(d/l*Test_d+l*Grad(d).Grad(Test_d)))") md.add_nonlinear_term(mim4, "(p+degradation*kappa*Div(u))*Test_p") md.add_nonlinear_term(mim4, "(Grad(C)+c2*C*Grad(p)).Grad(Test_C)" "+1e3*pos_part(2*d-1)*(C-C0)*Test_C") # loading and boundary conditions md.add_fem_data("dirichlet_data", mfu); ibdir = md.add_Dirichlet_condition_with_multipliers(mim9, "u", mfdir, TB_RG, "dirichlet_data") dirmultname = md.mult_varname_Dirichlet(ibdir) md.add_Dirichlet_condition_with_multipliers(mim4, "C", mfC, EXT_RG, "C0") print("Displacement dofs: %i\nTotal model dofs: %i" % (mfu.nbdof(), md.nbdof())) eps = 0. with open("gf_phasefield_forces.dat", "w") as f: for step in range(steps): eps_old = eps while True: # loop until a solution without temporary rate-dependance is achieved eps = eps_old + (step > 0)*strain_rate*md.variable("Dt") X = md.from_variables() print("Step %i with eps=%e" % (step, eps)) md.set_variable("dirichlet_data", md.interpolation("%.15g*[0;X(2)]" % eps, mfu)) try: nit, conv = md.solve("noisy", "lsolver", "mumps", "max_iter", 20, "max_res", 1e-9, "lsearch", "simplest", "alpha max ratio", 1e9, "alpha min", 1., "alpha mult", 0.1, "alpha threshold res", 1e9) except (KeyboardInterrupt, SystemExit): raise except: conv = False if conv: md.set_variable("u_prev", md.variable("u")) md.set_variable("v_prev", md.interpolation("(u-u_prev)/Dt", mfu)) md.set_variable("psi0_max", md.interpolation("max(psi0_max,psi0)", mimd9)) if nit <= 4: md.set_variable("Dt", 1.5*md.variable("Dt")) break else: md.to_variables(X) md.set_variable("Dt", 0.5*md.variable("Dt")) # post-processing output out = (mfu, md.variable("u"), "Displacements", mfd, md.variable("d"), "Damage", mfp, md.variable("p"), "Pressure", mfC, md.variable("C"), "C") for i,j in [[1,1],[2,2],[1,2]]: SIGMA = md.interpolation("degradation*(kappa*Div(u)*Id(2)+2*mu*deveps)(%i,%i)" % (i,j), mfout) out += (mfout, SIGMA, "Cauchy Stress %i%i" % (i,j), mfout, md.local_projection(mim9, "psi0_max", mfout), "psi0_max", mfout, md.interpolation("psi0", mfout), "psi0", mfout, md.interpolation("Gc", mfout), "Gc") mfout.export_to_vtk("gf_phasefield_%i.vtk" % step, *out) f.write(("step=%i average strain=%e reaction force=%e\n") % (step, eps, gf.asm_generic(mim9, 0, "%s(2)" % dirmultname, T_RG, md))) f.flush()