I am calculating the contact pressure between an elastic body (2D semicircle) and a rigid body (arc) due to thermal expansion. It is known that the temperature distribution of the elastic body is Delta_T. For 1D contact (uniform temperature, same expansion, same contact pressure on the circumference), it can be equivalent to point-to-point contact: in this case, binary search can be used. When the radius of the elastic body after expansion is greater than the radius of the rigid body, a uniformly distributed traction force P directed towards the center of the circle is applied on the circumference of the elastic body iteratively increased until there is no overlap of the radii. At this point, P is the contact pressure. This is my code, and it works well：
# 1D
from dolfin import *
# Create mesh
mesh = Mesh("fuel.xml")
boundary = MeshFunction("size_t", mesh, "fuel_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "fuel_physical_region.xml")
# Parameters
E = Constant(2e5) # MPa
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = Constant(1.2e-5)
def eps(v):
return 0.5*(grad(v) + grad(v).T)
def sigma(v, dT):
return (lmbda*tr(eps(v))- alpha*(3*lmbda+2*mu)*dT)*Identity(2) + 2.0*mu*eps(v)
ds = Measure('ds', domain=mesh, subdomain_data=boundary)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
Vu = VectorFunctionSpace(mesh, 'CG', 2)
du = TrialFunction(Vu)
u_ = TestFunction(Vu)
u = Function(Vu, name="Displacement")
bcu = DirichletBC(Vu.sub(0), Constant(0) , boundary, 3)
n = FacetNormal(mesh)
Delta_T =1000 #Uniform temperature
P = 100 #Uniform contact pressure
Wint = inner(sigma(du, Delta_T), eps(u_))*dx
aM = lhs(Wint)
LM = rhs(Wint) + dot(Constant(P) * n, u_)*ds(4)
solve(aM == LM, u, bcu)
For 2D contact (non-uniform temperature, different expansion, different contact pressure on the circumference), I want to mimic the above method. Although the expansion on the circumference is not uniform, I assume that all nodes expand and overlap with the rigid body. At this point, a discrete traction force vector P = (P1, P2… PN) is applied on the circumference, where N is the number of nodes. Using binary search, each component of P is simultaneously increased. When a component corresponding to a node no longer overlaps, that component remains unchanged, and the process is repeated for the remaining components. The resulting vector P is a rough estimate of the non-uniform contact pressure in 2D. It is no need to worry about the correctness of my idea, just that I encountered a problem in the implementation process. As shown in the above code, I can only implement 1D contact. For 2D contact, I cannot apply a non-uniform vector directed towards the center of the circle. I have learned from you the use of non-uniform heat flux density boundary conditions, which can assign different traction forces to nodes, but I cannot guarantee that they point towards the center of the circle, or the effect is not very good (because when I assign non-uniform discrete pressure to the nodes, no matter how I change the pressure component sizes, the displacement of the nodes after expansion remains almost unchanged). I hope everyone can provide assistance.（I believe the potential areas where errors may have occurred are:vertex_to_dof_map = vertex_to_dof_map(VT)）
# 2D：Additional code to modify P into a non-uniform vector based on the 1D code
from dolfin import *
import numpy as np
import pandas as pd
import joblib
# Create function space of Temperature
VT = FunctionSpace(mesh, "CG", 1)
# Define Uniform temperature boundary condition
Tf = Function(VT)
Tf_value = [984, 983, 982, 981, 980, 979, 978, 977, 976, 975, 974, 973, 972, 971, 970, 969, 968, 967, 966, 965, 964, 963, 962, 961, 960,
960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984]
f_coordinates = joblib.load('fuel_coordinates.joblib') # Coordinates of elastic body circumference nodes
vertex_to_dof_map = vertex_to_dof_map(VT) # Perhaps the problem lies here, where the assignment of traction force boundary conditions should be changed to Vu, but I tried it and the program reported an error.
class Coupleboundary(SubDomain):
def inside(self, x, on_boundary):
return boundary.where_equal(4) and on_boundary
f_s = Coupleboundary()
vertex_mark = MeshFunction("size_t", mesh, 0, 0)
f_s.mark(vertex_mark, 1)
vertices = vertex_mark.where_equal(1)
dof_coords = VT.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array([dof_coords[boundary_dof] for boundary_dof in boundary_dofs])
points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(f_coordinates, 0)
distances = np.sum(np.square(points_A - points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
Tf.vector()[boundary_dofs[row]] = Tf_value[col]
print(boundary_coords[row][1], Tf_value[col])
bcT = DirichletBC(VT, Tf, boundary, 4)
# Define variational problem
T_ = TestFunction(VT)
dT = TrialFunction(VT)
aT = dot(grad(dT), grad(T_)) * dx
LT = Constant(0) * T_ * dx
# Compute Delta_T
Delta_T = Function(VT, name="Temperature-fuel")
solve(aT == LT, Delta_T, bcT)
# Mechanical calculation
# Define Uniform traction boundary condition
P = Function(Vu)
P_0 = [100] * 50 # As an example of uniformity
for row, col in zip(*positions):
P.vector()[boundary_dofs[row]] = P_0[col]
Wint = inner(sigma(du, Delta_T), eps(u_))*dx
aM = lhs(Wint)
LM = rhs(Wint) + dot(P, u_)*ds(4)
solve(aM == LM, u, bcu)
This is the geo mesh file：
// Gmsh project created on Sun Apr 21 15:44:11 2024
SetFactory("OpenCASCADE");
N = 50;
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 7.06, 0, 1.0};
//+
Point(3) = {0, -7.06, 0, 1.0};
//+
Circle(1) = {2, 1, 3};
//+
Line(2) = {2, 3};
//+
Transfinite Curve {1} = N Using Progression 1;
//+
Curve Loop(1) = {2, -1};
//+
Plane Surface(1) = {1};
//+
Physical Curve("sym", 3) = {2};
//+
Physical Curve("surface", 4) = {1};
//+
Physical Surface("fuel", 5) = {1};
The geometry is as shown in the figure, the dashed part of the rigid body is not modeled, and the above process only needs to determine the relationship between the displacement after expansion and the radius of the rigid body:
