Forces and Moments

22.6. Forces and Moments#

We want to fix a beam at the left side $\Gamma_L$ by imposing mean value constraints:

\[ \int_{\Gamma_L} u_x = 0 \qquad \int_{\Gamma_L} u_y = 0 \]

Fix also the rotation around the origin via $$ \int_{\Gamma_L} u \wedge x = \int_{\Gamma_L} u_x y - u_y x = 0 $$

We minimize energy under these constraints. Therefore, we introduce the Lagrange function

\[ L(u, \lambda_x, \lambda_y, \lambda_M) = \frac{1}{2} \int \sigma(\varepsilon(u)) : \varepsilon(v) dx - \int f v dx + \frac{\lambda_x}{t} \int_{\Gamma_L} u_x ds + \frac{\lambda_y}{t} \int_{\Gamma_L} u_y ds + \frac{12 \lambda_M}{t^3} \int_{\Gamma_L} u \wedge x \, ds \]

The factors are chosen such that a linearized rigid body displacement

\[\begin{split} u(x,y) = \left( \begin{array}{c} a \\ b \end{array} \right) + c \left( \begin{array}{c} y \\ -x \end{array} \right) \end{split}\]

leads to $\lambda_x a + \lambda_y b + \lambda_M c$, which leads to the physical meaning of force and angular momentum of the Lagrange parameters.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw

thickness = 0.1
bar = MoveTo(0,-thickness/2).Rectangle(1,thickness).Face()
bar.edges.Min(X).name="left"
bar.edges.Max(X).name="right"

mesh = Mesh(OCCGeometry(bar, dim=2).GenerateMesh(maxh=thickness/3))
Draw (mesh);

E = 100    # Young modulus 
nu = 0.2   # Poisson ration (Querkontraktionszahl)
# Lame parameters via plane strain assumptions
mu = E/2/(1+nu)
lam = 2*nu*mu/(1-nu)

def stress(strain):
    return 2*mu*strain + lam*Trace(strain)*Id(mesh.dim)

Lagrange system:

Find $u, \lambda_x, \lambda_y, M$ such that

\[\begin{split} \begin{array}{ccccccccll} A(u,v) & + & \tfrac{\lambda_x}{t} \int_{\Gamma_L} v_x & + & \tfrac{\lambda_y}{t} \int_{\Gamma_L} v_y & + & \tfrac{12 \lambda_M}{t^3} \int_{\Gamma_L} v \wedge x & = & f(v) \quad & \forall \, v \\ \tfrac{\mu_x}{t} \int_{\Gamma_L} u_x &&&&&&& = & 0 & \forall \, \mu_x \\ \tfrac{\mu_y}{t} \int_{\Gamma_L} u_y &&&&&&& = & 0 & \forall \, \mu_y \\ \tfrac{12 \mu_M}{t^3} \int_{\Gamma_L} u \wedge x &&&&&&& = & 0 & \forall \, \mu_M \end{array} \end{split}\]

fes = VectorH1(mesh, order=3)
fesR = NumberSpace(mesh)

X = fes * fesR * fesR * fesR
u,lamx,lamy,lamM = X.TrialFunction()
v,mux,muy, muM = X.TestFunction()

a = BilinearForm(X)
a += InnerProduct(stress(Sym(Grad(u))), Sym(Grad(v)))*dx
a += 1/thickness*(u*CF((mux,muy)) + v*CF((lamx,lamy))) * ds("left")
a += 12/thickness**3*(u*muM+v*lamM) * CF((y,-x)) * ds("left")
a.Assemble();
inv = a.mat.Inverse(X.FreeDofs())

Volume force in horizontal direction, total force is $1N$:

f = LinearForm(1/thickness*v[0]*dx).Assemble()
gfu = GridFunction(X)
gfu.vec[:] = inv*f.vec
print ("Fx, Fy, M:", list(gfu.vec[-3:]))

Fx, Fy, M: [1.000000000000007, -2.0863152884223066e-15, 1.1758452081569025e-15]

Volume force in vertical direction, total force is $1 N$. Center of mass is $(0.5,0)$. Angular momentum is $0.5 Nm$.

f = LinearForm(1/thickness *v[1]*dx).Assemble()
gfu = GridFunction(X)
gfu.vec[:] = inv*f.vec
print ("Fx, Fy, M:", list(gfu.vec[-3:]))

Fx, Fy, M: [-5.865286511959579e-13, 1.0000000000036993, -0.5000000000029882]