22.7. Meta-material with negative \(\nu\)

Sebastian Hirnschall and Joachim Schöberl

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

t=.5
w=3.5
h=6
phi = 30

d = t/sin((90-phi)/360*2*pi)
a = w/tan((90-phi)/360*2*pi)
wp = MoveTo(w,a).LineTo(0,0).LineTo(0,h).LineTo(w,h-a).Finish().Offset(t/2)
f1 = Face(wp.Wire())

wp = MoveTo(w-t/2,a+d/2).LineTo(t/2,d/2).LineTo(t/2,h-d/2).LineTo(w-t/2,h-a-d/2).LineTo(w-t/2,h-a+d/2).LineTo(t/2,h+d/2).LineTo(-t/2,h+d/2).LineTo(-t/2,-d/2).LineTo(+t/2,-d/2).LineTo(w-t/2,a-d/2).Close()
#DrawGeo(f1);
f1=wp.Face()
DrawGeo(f1);

faces = []
nx,ny = 9,7
for j in range(ny):
    for i in range(nx):
        faces.append(f1.Move( (2*i*w,2*j*(h-a),0) ))
    for i in range(nx):
        faces.append(f1.Move( ((2*i+1)*w,(2*j+1)*(h-a),0)))
        
faces = sum(faces)

DrawGeo (faces);

wp = MoveTo(3.255,3.5).Rectangle(49.490,43).Reverse().MoveTo(-10,-10).Rectangle(100,100).Face()
geo = faces-wp
DrawGeo(geo);

sol = Prism(geo, (0,0,5))

xmin = sol.faces.Min(X).center[0]
xmax = sol.faces.Max(X).center[0]
print ("xmin/max =", xmin, xmax)

sol.faces[X<xmin+1e-4].name = "left"
sol.faces[X>xmax-1e-4].name = "right"

xmin/max = 3.2549999999999955 52.745000000000005

DrawGeo (sol);

mesh = Mesh(OCCGeometry(sol).GenerateMesh(maxh=1))
mesh.ne

40092

Draw (mesh);

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

force = CF( (1,0,0) )


def Stress(strain):
    return 2*mu*strain + lam*Trace(strain)*Id(3)    

fes = VectorH1(mesh, order=2, dirichlet="left")
u,v = fes.TnT()
gfu = GridFunction(fes)

with TaskManager():
    a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
    a.Assemble()
    inv = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
    f = LinearForm(force*v*ds("right")).Assemble()
    
gfu = GridFunction(fes)
gfu.vec.data = inv * f.vec

Draw (gfu);

observe the limits of geometrically linearized elasticity