Computation of reaction force

1

Dear all,

I wish to compute the reaction force for the multiphysics (thermoelasticity: thermal is transient and elasticity is static) problem at every time instant. I am fixing the left edge and pulling the right edge with displacement boundary condition. Also the temperature is applied to left edge. I am solving these coupled temperature-displacement using NR solver of FEniCS.

I wish to calculate the reaction force on the left edge. I have tried to calculate it from here. But it shows the reaction forces as 0.0 which is not true. I tried other way also as shown in the code. Both of them show the reaction force as 0.0. Is there something I am missing. Any help would be valuable.

from dolfin import *
from fenics import *
E = 30;
nu = 0.3000000000
alpha = 0.1;
mu = E/(2.0*(1.0+nu)) 
lamb = E*nu/((1.0+nu)*(1.0-2.0*nu))
coeff = E*alpha/((1.0-2.0*nu)*3.0);
initial=Constant(0.0)
T = 0.1            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size

def clamp(x,on_boundary):
            return on_boundary and near(x[0],0.0000)
def right(x,on_boundary):
    	return on_boundary and near(x[0],1.00)
mesh=RectangleMesh(Point(0.0,0.0),Point(1.0,1.0),40,40)
V = VectorElement("CG", mesh.ufl_cell(), 2) #displacement u
Q = FiniteElement("CG", mesh.ufl_cell(), 2) #Temperature
a = MixedElement([V,Q])
W = FunctionSpace(mesh, a)
w = Function(W) 
(u, p) = split(w)
(v, q) = TestFunctions(W)
bc1=DirichletBC(W.sub(0),Constant((0.0,0.0)),clamp)  
bc2=DirichletBC(W.sub(0).sub(0),Constant((0.01)),right)
bc3=DirichletBC(W.sub(1),Constant(1.0),clamp)    
bc=[bc1,bc2,bc3]      
d=2;    
def epsilon(u):
     return 0.5*(nabla_grad(u)+nabla_grad(u).T) 
 
def sigma(u,p):
      return lamb*nabla_div(u)*Identity(d)+2*mu*epsilon(u)-coeff*(p-initial)*Identity(d)

def sigmah(u,p):
      return tr(sigma(u,p))/3.0

p_n=project(Constant(0.0),W.sub(1).collapse())
F=inner(sigma(u,p),epsilon(v))*dx + p*q*dx + dt*inner(grad(p), grad(q))*dx - (p_n)*q*dx


x_dofs = W.sub(0).sub(0).dofmap().dofs()
y_dofs = W.sub(0).sub(1).dofmap().dofs()
t = 0.0
for n in range(num_steps):
        t += dt
        solve(F == 0, w, bc)
        (u,p)=w.split(True)
        p_n.assign(p)
        f_int = assemble(inner(sigma(u,p),epsilon(v))*dx) 
        Fx = 0
        for i in x_dofs:
              Fx += f_int[i]

        Fy = 0
        for i in y_dofs:
              Fy += f_int[i]

        print("Horizontal reaction force: %f"%(Fx))
        print("Vertical reaction force: %f"%(Fy))

        R = derivative(F, u, v)
        b=assemble(-R)
        fint=b.copy()
        Fy = 0
        for i in x_dofs:
            Fy += fint[i]
        print("######################## %f"%(Fy))

Thanks in advance.

2

Because all external forces come from BCs, and the solution is in static equilibrium, you expect the net reaction force to all BCs to be zero. If you want to get the reaction force from one BC, you need to apply the others to the residual. For example, to get the reaction to the force on the right end, you could do something like the following:

for n in range(num_steps):
        t += dt
        solve(F == 0, w, bc)
        (u,p)=w.split(True)
        p_n.assign(p)
        f_int = assemble(inner(sigma(u,p),epsilon(v))*dx)

        # Need to apply one of the BCs to isolate the force from the other.
        # (We expect the net reaction to all BCs to be zero, because the
        # problem is in equilibrium.)

        # Applying BC on left end and expecting horizontal reaction force
        # pointing to the left, to balance force exerted by BC on right end,
        # which has not been applied to the residual:
        bc1.apply(f_int)
        
        Fx = 0
        for i in x_dofs:
              Fx += f_int[i]

        Fy = 0
        for i in y_dofs:
              Fy += f_int[i]

        print("Horizontal reaction force: %f"%(Fx))
        print("Vertical reaction force: %f"%(Fy))

        # Equivalent:
        
        #R = derivative(F, u, v)
        #b=assemble(-R)
        b = assemble(F)
        bc1.apply(b)
        
        fint=b.copy()
        Fx = 0
        for i in x_dofs:
            Fx += fint[i]
        print("######################## %f"%(Fx))
3 Likes
3

Thanks a lot! Now it makes much more sense.