Mixed Linear Elasticity and Helmholtz equation


436
views
0
11 months ago by
Hello,

I'm trying to solve a topology optimization problem using the density filtering technique proposed by Lazarov. To this extend I've arrived at the following set of equations:

 $r^2\int_{\Omega}\left(\nabla\rho\cdot\nabla w\right)dx+\int_{\Omega}\left(\rho-\rho_N\right)dx=0$r2Ω(ρ·w)dx+Ω(ρρN)dx=0 
 $\int_{\Omega}SIMP\left(\rho\right)\frac{\mu}{2\lambda}\left(\nabla u+\nabla u^T\right)\cdot\left(\nabla v+\nabla v^T\right)dx+\int_{\Omega}SIMP\left(\rho\right)\left(\nabla\cdot u\right)\left(\nabla\cdot v\right)dx-\int_{\Omega}\left(\frac{1}{\lambda}\right)\left(b\cdot v\right)dx=0$ΩSIMP(ρ)μ2λ (u+uT)·(v+vT)dx+ΩSIMP(ρ)(·u)(·v)dxΩ(1λ )(b·v)dx=0 

where:

 $SIMP\left(\rho\right)=\epsilon+\left(1-\epsilon\right)\rho^p$SIMP(ρ)=ϵ+(1ϵ)ρp 

I sum both equations to have:

 $F\left(\left(u,\rho\right),\left(v,w\right)\right)=0$F((u,ρ),(v,w))=0 

And I try to solve this equation using Newton's method. However, I can't get the method to converge. Could someone help me with this issue? I'm attaching the pertinent part of the code below.

Thanks in advance!
Eduardo

from fenics import *

solver_parameters = {"newton_solver": {"maximum_iterations": 50}}

# Structure material properties
E = 1
nu = 0.2

# Geometry of the design domain
L = 4
H = 1
b = 0.01

# Load
F = 0.2

# SIMP constants
p = Constant(3.0) # penalisation
eps = Constant(1.0e-3)

# Topology Optimization constants
rho_0 = Constant(1.0) # Initial guess
r_min = Constant(0.02) # Minimum length

# FEM constants
el_size = 0.05
nx = int(L/el_size)
ny = int(H/el_size)

mesh = RectangleMesh(Point(0, 0), Point(L, H), \
                     nx, ny, 'right')

# Define function space and base functions
V = VectorElement("CG", mesh.ufl_cell(), 2)
W = FiniteElement("CG", mesh.ufl_cell(), 1)
M = FunctionSpace(mesh, V*W)
V = M.sub(0)
W = M.sub(1)

# Boundary Condition
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

left = Left()
bc = DirichletBC(V, Constant((0.0, 0.0)), left)

# Load
b = Expression(('q_x', 'q_y'), degree = 2, \
               q_x = 0.0, q_y = -F/(b*H))

dm = TrialFunction(M)
m = Function(M)
(u, rho) = split(m)
(v, w) = TestFunctions(M)

rho_n = interpolate(rho_0, W.collapse())

simp = eps + (1-eps)*(rho**p)
lmbda = nu*E/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))
K_1 = Constant(2*mu/lmbda)
K_2 = Constant(1/lmbda)

F = K_1*simp*inner(sym(nabla_grad(u)), \
                   sym(nabla_grad(v)))*dx + \
    simp*div(u)*div(v)*dx - K_2*inner(b, v)*dx + \
    (r_min**2)*inner(grad(rho), grad(w))*dx + \
    (rho-rho_n)*w*dx
J = derivative(F, m, dm)

s = Function(M)
solve(F == 0, s, bcs=bc, J=J, solver_parameters=solver_parameters)
(u, rho) = s.split()
​
Community: FEniCS Project

1 Answer


2
11 months ago by
I solve the Helmholtz eq. first and then the elasticity eq. instead of solving it all at once. You probably also like to do a "projection" after filtering to avoid intermediate densities.
I intended to use the problem solution in dolfin-adjoint. When I tried to solve the Helmholtz equation first, it complained that the design variable was not part of the forward problem. Could you show me how you did it with an example? Thank you!
written 11 months ago by EduardoMS  
My code is quite long and messy, but here are some bits and pieces from where you hopefully can get the idea. I use a continuous space (Vc) for the filtered density and a discontinuous (Vd) for the design variables and physical density.
----
def Heavy(beta,eta,x):
    return (tanh(beta*eta)+tanh(beta*(x-eta)))/(tanh(beta*eta)+tanh(beta*(1-eta)))

def filtering(rho,rmin):
   trial_rhoFil=TrialFunction(Vc)
   test_rhoFil=TestFunction(Vc)
   F=(inner(pow(rmin,2)*grad(trial_rhoFil),grad(test_rhoFil))+inner(trial_rhoFil,test_rhoFil))*dx
   Af=assemble(F)
   b=assemble(inner(rho,test_rhoFil)*dx)
   solve(Af,rhoFil.vector(),b,KrylovSolver,KrylovPrecon)
   return rhoFil

def projection(rhoFil,beta,eta):
   tmp=project(rhoFil,Vd,solver_type=KrylovSolver,preconditioner_type=KrylovPrecon)
   tmp2=project(Heavy(beta,eta,tmp),Vd,solver_type=KrylovSolver,preconditioner_type=KrylovPrecon)
   return tmp2

def stressStrain(u,rhoPhys):
    epsilon=sym(grad(u))
    E=Emin+pow(rhoPhys,pen)*(E0-Emin)
    mu = E/(2*(1+nu))
    lam = E*nu/((1+nu)*(1-2*nu))
    return (2*mu*epsilon+lam*tr(epsilon)*Identity(u.geometric_dimension()),epsilon)

def Forward(rho,beta,eta):
   rhoFil=filtering(rho,rmin)
   rhoPhys=projection(rhoFil,beta,eta)

   # Mechanics Problem
   u = Function(V,name='Displacement')
   trial_u=TrialFunction(V)
   t = TestFunction(V)
   (sigma,epsilon)=stressStrain(trial_u,rhoPhys)
   F2=inner(sigma, grad(t))*dx
   L=dot(Load,t)*ds(1)
   A, b = assemble_system(F2,L,Dbc)
   solve(A,u.vector(),b,KrylovSolver,KrylovPreconElastic)
   return (u,rhoPhys,rhoFil)
----
      parameters["adjoint"]["stop_annotating"] = False
      adj_reset()
      (u,rhoPhys,rhoFil)=Forward(rho,beta,eta)
      parameters["adjoint"]["stop_annotating"] = True
      
      (eps,sig)=stressStrain(u,rhoPhys)
      J=Functional(0.5*inner(sig,eps)*dx)
      Vol=Functional((rhoPhys-Volfrac)*dx)
      m=Control(rho)
      Jhat = ReducedFunctional(J,m)
      Vhat = ReducedFunctional(Vol,m)
----
written 11 months ago by SMadsen  
Thank you! Now it worked for me. When I attempted this approach, I was using FEniCS solver for density filtering equation instead of dolfin_adjoint solver. Your code really helped me!
written 11 months ago by EduardoMS  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »