#!/usr/bin/env python # -*- coding: utf-8 -*- # Python GetFEM interface # # Copyright (C) 2015-2020 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 3 of the License, or # (at your option) any later version along with the GCC Runtime Library # Exception either version 3.1 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 and GCC Runtime Library Exception 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. # ############################################################################ # # Contact of a deformable 'wheel' onto a plane deformable obstacle. # ############################################################################ import numpy as np import getfem as gf gf.util_trace_level(1) # No trace for mesh generation nor for assembly # # Physical parameters # E = 21E6 # Young Modulus (N/cm^2) nu = 0.3 # Poisson ratio clambda = E*nu/((1+nu)*(1-2*nu)) # First Lame coefficient (N/cm^2) cmu = E/(2*(1+nu)) # Second Lame coefficient (N/cm^2) applied_force = 1E7 # Force at the hole boundary (N) # # Numerical parameters # h = 1 # Approximate mesh size elements_degree = 2 # Degree of the finite element methods r_aug = E; # Augmentation parameter for the augmented Lagrangian # # Mesh generation. Meshes can also been imported from several formats. R_o = 15. R_i = 8. mo1 = gf.MesherObject('ball', [0., R_o], R_o) mo2 = gf.MesherObject('ball', [0., R_o], R_i) mo3 = gf.MesherObject('set minus', mo1, mo2) print('Meshes generation') mesh1 = gf.Mesh('generate', mo3, h, 2) #mesh1.translate([0.,-0.1]) mesh2 = gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[30,10];NOISED=0;NSUBDIV=[%d,%d];' % (int(30/h)+1, int(10/h)+1)); mesh2.translate([-15.,-10.]) # Boundary selection # fb1 = mesh1.outer_faces_in_box([-8.1, 6.9], [8.1, 23.1]) # Boundary of the hole fb2 = mesh1.outer_faces_with_direction([0., -1.], np.pi/4.5) # Contact boundary of the wheel fb3 = mesh2.outer_faces_with_direction([0., -1.], 0.01) # Bottom boundary of the foundation fb4 = mesh2.outer_faces_with_direction([0., 1.], 0.01) # Top boundary of the foundation HOLE_BOUND=1; CONTACT_BOUND=2; BOTTOM_BOUND=3; TOP_BOUND=4; mesh1.set_region(HOLE_BOUND, fb1) mesh1.set_region(CONTACT_BOUND, fb2) mesh1.region_subtract(CONTACT_BOUND, HOLE_BOUND) mesh2.set_region(BOTTOM_BOUND, fb3) mesh2.set_region(TOP_BOUND, fb4) print(mesh2.region(TOP_BOUND)) # # Definition of finite elements methods and integration method # mfu1 = gf.MeshFem(mesh1, 2) mfu1.set_classical_fem(elements_degree) mflambda = gf.MeshFem(mesh1, 2) mflambda.set_classical_fem(elements_degree-1) mflambda_C = gf.MeshFem(mesh1, 2) mflambda_C.set_classical_fem(elements_degree-1) mfu2 = gf.MeshFem(mesh2, 2) mfu2.set_classical_fem(elements_degree) mfvm1 = gf.MeshFem(mesh1, 1) mfvm1.set_classical_discontinuous_fem(elements_degree) mfvm2 = gf.MeshFem(mesh2, 1) mfvm2.set_classical_discontinuous_fem(elements_degree) mim1 = gf.MeshIm(mesh1, 4) mim1c = gf.MeshIm(mesh1, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)')) mim2 = gf.MeshIm(mesh2, 4) # # Model definition # md=gf.Model('real'); md.add_fem_variable('u1', mfu1) # Displacement of the structure 1 md.add_fem_variable('u2', mfu2) # Displacement of the structure 2 #md.add_initialized_data('cmu', [cmu]) #md.add_initialized_data('clambda', [clambda]) #md.add_isotropic_linearized_elasticity_brick(mim1, 'u1', 'clambda', 'cmu') #md.add_isotropic_linearized_elasticity_brick(mim2, 'u2', 'clambda', 'cmu') md.add_initialized_data('matparams', [clambda,cmu]) md.add_finite_strain_elasticity_brick(mim1, 'Compressible_Neo_Hookean', 'u1', 'matparams') md.add_finite_strain_elasticity_brick(mim2, 'Compressible_Neo_Hookean', 'u2', 'matparams') md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', elements_degree-1, BOTTOM_BOUND) # Contact condition (augmented Lagrangian) md.add_initialized_data('r', [1e0*r_aug]) md.add_initialized_data('f_coeff', [0.5]) md.add_filtered_fem_variable('lambda1', mflambda_C, CONTACT_BOUND) release_dist = 1 contact_implementation = 0 if contact_implementation == 0: md.add_raytracing_transformation("contact_trans", release_dist) md.define_variable_group("u", "u1", "u2") md.add_slave_contact_boundary_to_raytracing_transformation("contact_trans", mesh1, "u", CONTACT_BOUND) md.add_master_contact_boundary_to_raytracing_transformation("contact_trans", mesh2, "u", TOP_BOUND) md.add_nonlinear_term(mim1c, "-lambda1.Test_u1", CONTACT_BOUND) md.add_nonlinear_term(mim1c, "Interpolate_filter(contact_trans, lambda1.Interpolate(Test_u,contact_trans), 1)", CONTACT_BOUND) md.add_nonlinear_term(mim1c, "-(1/r)*lambda1.Test_lambda1", CONTACT_BOUND) md.add_nonlinear_term(mim1c, "Interpolate_filter(contact_trans," "(1/r)*Coulomb_friction_coupled_projection(lambda1," "Transformed_unit_vector(Grad_u1, Normal)," "u1-Interpolate(u,contact_trans)," "(Interpolate(X,contact_trans)" "+Interpolate(u,contact_trans)-X-u1)" ".Transformed_unit_vector(Grad_u1, Normal)," "f_coeff," "r).Test_lambda1, 1)", CONTACT_BOUND) elif contact_implementation == 1: ind = md.add_integral_large_sliding_contact_brick_raytracing("r", release_dist, "f_coeff", "1", 0) md.add_slave_contact_boundary_to_large_sliding_contact_brick(ind, mim1c, CONTACT_BOUND, "u1", "lambda1", "") md.add_master_contact_boundary_to_large_sliding_contact_brick(ind, mim2, TOP_BOUND, "u2", "") else: md.add_integral_contact_between_nonmatching_meshes_brick\ (mim1c, 'u1', 'u2', 'lambda1', 'r', 'f_coeff', CONTACT_BOUND, TOP_BOUND) md.add_initialized_data('xc', 0.) md.add_initialized_data('yc', R_o) md.add_variable('dy', 1) md.set_variable('dy', -0.3) md.add_variable('phi',1) md.set_variable('phi', 0.01) md.disable_variable('dy') md.disable_variable('phi') md.add_initialized_data('area', gf.asm_generic(mim1, 0, '1', HOLE_BOUND)) # Calculate the area of the region md.add_initialized_data('Fy', -applied_force) md.add_macro("rot(x)", "[cos(x)-1,sin(x);-sin(x),cos(x)-1]") md.add_macro("rot_der(x)", "[-sin(x),cos(x);-cos(x),-sin(x)]") md.add_linear_term(mim1, 'Fy/area*Test_dy', HOLE_BOUND) md.add_initialized_data('Mz', 200) md.add_filtered_fem_variable('mult', mflambda, HOLE_BOUND) md.add_nonlinear_term(mim1, '(u1-[0;dy]-rot(phi)*(X-[xc;yc])).Test_mult', HOLE_BOUND) md.add_nonlinear_term(mim1, 'Test_u1.mult', HOLE_BOUND) # # Model solve # print('Solve problem with ', md.nbdof(), ' dofs') md.solve('max_res', 1E-11, 'max_iter', 40, 'noisy' , 'lsearch', 'simplest', 'alpha min', 0.8) # print('Contact multiplier ', md.variable('lambda1')) # # Solution export # U1 = md.variable('u1') U2 = md.variable('u2') #VM1 = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u1', 'clambdastar', 'cmu', mfvm1) #VM2 = md.compute_isotropic_linearized_Von_Mises_or_Tresca('u2', 'clambdastar', 'cmu', mfvm2) VM1 = md.compute_finite_strain_elasticity_Von_Mises('Compressible_Neo_Hookean', 'u1', 'matparams', mfvm1) VM2 = md.compute_finite_strain_elasticity_Von_Mises('Compressible_Neo_Hookean', 'u2', 'matparams', mfvm2) mfvm1.export_to_vtk('displacement_with_von_mises1.vtk', mfvm1, VM1, 'Von Mises Stresses', mfu1, U1, 'Displacements', md.mesh_fem_of_variable("lambda1"), md.variable("lambda1"), 'lambda1', md.mesh_fem_of_variable("mult"), md.variable("mult"), 'mult') mfvm2.export_to_vtk('displacement_with_von_mises2.vtk', mfvm2, VM2, 'Von Mises Stresses', mfu2, U2, 'Displacements') print('You can view solutions with for instance:\nmayavi2 -d displacement_with_von_mises1.vtk -f WarpVector -m Surface -d displacement_with_von_mises2.vtk -f WarpVector -m Surface')