Computation of reaction force

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