#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from math import pi from scipy import sparse import getfem as gf import gmsh gf.util('trace level', 1) def GenerateFullDisk(outputFile, diskRadius, meshSize): gmsh.initialize() model_name=outputFile gmsh.model.add(model_name) cm = 1. Lc1 = meshSize #Create points factory = gmsh.model.occ factory.addPoint(0 * cm, 0 * cm, 0, Lc1, 1) factory.addPoint(diskRadius * cm, diskRadius * cm, 0, Lc1, 2) factory.addPoint(0 * cm, 2*diskRadius * cm, 0, Lc1, 3) factory.addPoint(-diskRadius * cm, diskRadius * cm, 0, Lc1, 4) factory.addPoint(0 * cm, diskRadius * cm, 0, Lc1, 5) #Create lines factory.addCircleArc(1, 5, 2, 10) factory.addCircleArc(2, 5, 3, 11) factory.addCircleArc(3, 5, 4, 12) factory.addCircleArc(4, 5, 1, 13) factory.addCurveLoop([10,11,12,13], 20) #Create surfaces factory.addPlaneSurface([20], 30) factory.synchronize() #Crete physical groups (for boundary conditions) gmsh.model.addPhysicalGroup(2, [30], 100) gmsh.model.addPhysicalGroup(1, [11,12], 102) gmsh.model.addPhysicalGroup(1, [10,13], 103) #Assign a mesh size to all the points: #gmsh.model.mesh.setSize(gmsh.model.getEntities(0), Lc1) gmsh.model.mesh.setAlgorithm(2, 30, 2) gmsh.option.setNumber("Mesh.ElementOrder", 2) gmsh.model.mesh.generate(2) gmsh.option.setNumber("Mesh.MshFileVersion", 2) gmsh.write("%s.msh" % model_name) gmsh.finalize() def SolverHertzContactProblem(mesh,materialParameters,boundaryTags,ponctualForce): mfu = gf.MeshFem(mesh, 2) mfu.set_classical_fem(3) mflambda = gf.MeshFem(mesh, 1) mflambda.set_classical_fem(2) mim = gf.MeshIm(mesh, 4) mim_c=gf.MeshIm(mesh, gf.Integ('IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(4),2)')) model = gf.Model('real') model.add_fem_variable('u', mfu) #Define behaviour law E,nu=materialParameters["Young"],materialParameters["Poisson"] clambda = E*nu/((1+nu)*(1-2*nu)) cmu = E/(2*(1+nu)) clambdastar = clambda # Plane strain assumption model.add_initialized_data('cmu', [cmu]) model.add_initialized_data('clambdastar', [clambdastar]) model.add_linear_generic_assembly_brick(mim, "clambdastar*(Div_u*Div_Test_u)+cmu*((Grad_u+(Grad_u)'):Grad_Test_u)") #Define boundary conditions dataunitv_x = np.array([1,0]) model.add_initialized_data('pts_vec', np.array([0.0,100.0])) model.add_initialized_data('dataunitv_x', dataunitv_x) model.add_pointwise_constraints_with_multipliers('u', 'pts_vec', 'dataunitv_x') model.add_initialized_data('pts_vec2', np.array([0.0,300.0])) model.add_pointwise_constraints_with_multipliers('u', 'pts_vec2', 'dataunitv_x') mfu = model.mesh_fem_of_variable('u') basicDofNodes=mfu.basic_dof_nodes() id_x,id_y=np.where( np.abs(basicDofNodes[1,:]-400.0) < 1e-6)[0] explicitRhs=np.zeros_like(basicDofNodes[0,:]) explicitRhs[id_x]=0 explicitRhs[id_y]=-ponctualForce model.add_explicit_rhs('u', explicitRhs) #Define unilateral contact condition model.add_initialized_data('r', [0.2]) model.add_initialized_data('N1', [0.,-1.0]) model.add_filtered_fem_variable("lambda", mflambda, boundaryTags["CONTACT_BOUNDARY"]) model.add_linear_generic_assembly_brick(mim_c, "lambda*(Test_u.N1)", boundaryTags["CONTACT_BOUNDARY"]) model.add_nonlinear_generic_assembly_brick(mim_c, "-(1/r)*(-lambda +neg_part(-lambda-r*(u.N1-X(2))))*Test_lambda", boundaryTags["CONTACT_BOUNDARY"]) #Solve problem model.solve('max_res', 1E-5,'very noisy','max_iter',60) print(gf.asm_generic(mim, 0, "lambda", boundaryTags["CONTACT_BOUNDARY"], model)) return model def ExportSolution(exportFile,model): displacement,mfu = model.variable('u'),model.mesh_fem_of_variable('u') contactMultipliers,mflambda= model.variable('lambda'),model.mesh_fem_of_variable('lambda') mfu.export_to_pos(exportFile, mfu, displacement,'Displacement',mflambda,-contactMultipliers,"Lagrange multiplier") def ExtractContactSolution(model,mesh,boundaryTags): contactMultipliers,mflambda= model.variable('lambda'),model.mesh_fem_of_variable('lambda') mflambda=model.mesh_fem_of_variable('lambda') # dummyModel=gf.Model('real') # mim = gf.MeshIm(mesh, 4) # dummyModel.add_filtered_fem_variable("ulambda", mflambda,boundaryTags["CONTACT_BOUNDARY"]) # dummyModel.add_linear_generic_assembly_brick(mim, 'ulambda.Test_ulambda', boundaryTags["CONTACT_BOUNDARY"]) # dummyModel.assembly() # data=dummyModel.tangent_matrix().csc_val() # rows,nbval=dummyModel.tangent_matrix().csc_ind() # nbdof=rows.shape[0]-1 # massMat=sparse.csc_matrix((data, nbval, rows), shape=(nbdof, nbdof)) # contactMultipliers=-massMat.dot(contactMultipliers) refIndices=mflambda.basic_dof_on_region(boundaryTags["CONTACT_BOUNDARY"]) nodesOnContactBoundary=mflambda.basic_dof_nodes().T[refIndices] xCoord=nodesOnContactBoundary[:,0] xAxisOrder=np.argsort(xCoord) return xCoord[xAxisOrder],contactMultipliers[xAxisOrder] def ComputeAnalyticalSolution(xCoord,diskRadius,verticalForce,materialParameters): E,nu=materialParameters["Young"],materialParameters["Poisson"] maxRadiusContact=np.power(3*verticalForce*diskRadius*(1-nu**2)/(4*E), 1.0/3.0) maxPressure=3*verticalForce/(2*pi*maxRadiusContact**2) pressureRepartition=1- (xCoord/maxRadiusContact)**2 pressureRepartition[pressureRepartition < 0] = 0 analyticalContactSolution=maxPressure*np.sqrt(pressureRepartition) return analyticalContactSolution def ComputeAnalyticalDeltaDisp(diskRadius,verticalForce,materialParameters): E,nu=materialParameters["Young"],materialParameters["Poisson"] deltaDisp=np.power( (9*verticalForce**2) / (16*diskRadius* (E/(1-nu**2))**2) ,1.0/3.0) return deltaDisp def ComputeNumericalDeltaDisp(model,diskRadius): displacement,mfu = model.variable('u'),model.mesh_fem_of_variable('u') basicDofNodes=mfu.basic_dof_nodes() _,id_y=np.where( np.abs(basicDofNodes[1,:]-2*diskRadius) < 1e-6)[0] return displacement[id_y] if __name__ =="__main__": diskRadius=200 GenerateFullDisk(outputFile="Full2DDisk", diskRadius=diskRadius, meshSize=10) mesh=gf.Mesh('import', 'gmsh', "Full2DDisk.msh") materialParameters={"Young":25e6,"Poisson":0.25 } boundaryTags={"CONTACT_BOUNDARY":103,"TOP":102} verticalForce=20 model=SolverHertzContactProblem(mesh=mesh,materialParameters=materialParameters,boundaryTags=boundaryTags,ponctualForce=verticalForce) ExportSolution("Full2DDiskSolution.pos",model) xCoord,numericalSolution=ExtractContactSolution(model=model,mesh=mesh,boundaryTags=boundaryTags) analyticalContactSolution=ComputeAnalyticalSolution(xCoord,diskRadius,verticalForce,materialParameters) deltaDispAnalytical=ComputeAnalyticalDeltaDisp(diskRadius,verticalForce,materialParameters) deltaDispNumerical=ComputeNumericalDeltaDisp(model,diskRadius) #print(deltaDispAnalytical," ",deltaDispNumerical) print(numericalSolution) print(analyticalContactSolution)