### Mixed Linear Elasticity and Helmholtz equation

279

views

0

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$

$\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$∫

where:

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

I sum both equations to have:

$F\left(\left(u,\rho\right),\left(v,w\right)\right)=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

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$

`r`^{2}∫_{Ω}(∇`ρ`·∇`w`)`d``x`+∫_{Ω}(`ρ`−`ρ`_{N})`d``x`=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$∫

_{Ω}`S``I``M``P`(`ρ`)`μ`2λ (∇`u`+∇`u`^{T})·(∇`v`+∇`v`^{T})`d``x`+∫_{Ω}`S``I``M``P`(`ρ`)(∇·`u`)(∇·`v`)`d``x`−∫_{Ω}(1λ )(`b`·`v`)`d``x`=0where:

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

`S``I``M``P`(`ρ`)=`ϵ`+(1−`ϵ`)`ρ`^{p}I sum both equations to have:

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

`F`((`u`,`ρ`),(`v`,`w`))=0And 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

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
6 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
6 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
6 months ago by
EduardoMS

Please login to add an answer/comment or follow this question.