4. Meta-material with negative Poisson’s ratio \(\nu\)

Sebastian Hirnschall and Joachim Schöberl

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

4.1. Geometry with OCC

We will first construct a single cell before building the complete geometry with a pattern.

4.1.1. Single Cell

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()
#Draw(f1);
f1=wp.Face()
Draw(f1);

4.1.2. Whole Geometry with Pattern

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)
Draw (faces);

4.1.3. Cleanup

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

4.1.4. Boundaries

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<xmin+1e-4].col = (112/255, 163/255, 204/255)
sol.faces[X>xmax-1e-4].name = "right"
sol.faces[X>xmax-1e-4].col = (204/255, 154/255, 121/255)

xmin/max = 3.2549999999999955 52.745000000000005

Draw (sol);

4.2. Mesh

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

40543

#mesh.ngmesh.Export('metamaterial.stl','STL Format')

Draw (mesh);

4.3. Solve Linear Elasticity

Similar to https://docu.ngsolve.org/ngs24/SaS/linearelasticity.html

Note that there are limitations of linear elasticity visible for large deformations.

E, nu = 210, 0.2 #material parameters for base material
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);