The model is defined on the square $\Omega$ with length $L=1\,\text{m}$. An homogenous heat source of $q=100\,\text{W}/\text{m}^2$ is applied to the solution domain. Simple isotropic Fourier's law $$ \boldsymbol{h} = \kappa \mathrm{grad}\left( T \right) $$ connects diffusive heat flux $\boldsymbol{h}$ and temperature gradient via the scalar conductivity $\kappa=400\,\mathrm{W}/\mathrm{m\,K}$.
Homogenous Neumann boundary conditions are prescribed at the south and east edges.
A constant Dirichlet temperature boundary condition $T_D= 0\,\mathrm{K}$ is prescribed on the west and north edges.
Then, the Poisson-like PDE is defined as:
$$\begin{align} \mathrm{div}\left( \boldsymbol{h} \right) + q = 0 \quad &\in \Omega,\\ \boldsymbol{h}\cdot\boldsymbol{n} = 0 \quad &\in \partial\Omega_\text{N},\\ T = T_\text{D} \quad &\in \partial\Omega_\text{D}. \end{align}$$
import torch
from scipy.optimize import bisect
from tqdm import tqdm
from torchfem import PlanarHeat
from torchfem.materials import IsotropicConductivity2D
from torchfem.mesh import rect_quad
torch.set_default_dtype(torch.double)
Base problem without optimization¶
Material model
# Thermal conductivity
kappa = 400.0
# Material model
material = IsotropicConductivity2D(kappa=kappa)
Geometry and model
# Edge length of domain
L = 1.0
# Number of elements per edge
n_elements = 200
# Create mesh
nodes, elements = rect_quad(n_elements + 1, n_elements + 1, L, L)
# Create model
model = PlanarHeat(nodes, elements, material)
Boundary conditions
# Heat source
q = 1000
# Dirichlet temperature
T_D = 0.0
# Apply Dirichlet boundary conditions at west and north edges
west = torch.isclose(nodes[:, 0], nodes[:, 0].min())
north = torch.isclose(nodes[:, 1], nodes[:, 1].max())
dirichlet_edges = west | north
model.constraints[dirichlet_edges] = True
model.temperatures[dirichlet_edges] = T_D
# Distribute the total heat source uniformly over the domain
areas = model.integrate_field()
model.heat_flux = model.integrate_body_load(q / areas.sum())
Solve the basic model
T, *_ = model.solve()
Visualize basic temperature distribution
model.plot(
node_property={"Temperature (T)": T},
cmap="magma",
linewidth=0.0,
bcs=False,
colorbar=True,
)
Optimization¶
Instead of sensitivity filtering, we penalize the conventional objective $\mathcal{F}$ using Tikhonov regularization.
$$\mathcal{F} = \underline{T}\, \underline{f}^\top + \alpha \int_\Omega \left\lVert \frac{\partial \rho} {\partial \boldsymbol{x}} \right\rVert^2 \, \mathrm{d}V $$
The first term is the "thermal compliance" Yet, some authors argue that it bears little physical interpretability (in contrast to mechanical compliance) and adjusting the objective to minimizing the temperature variance being a more natural choice [1].
Under a minimal volume constraint the optimization problem reads:
$$ \min_{\rho \left( \boldsymbol{x} \right)} \left( \mathcal{F} \right)\quad\text{s.t.}\, \int_\Omega \rho \, \mathrm{d}V \leq V_0. $$
To compute meaningful spatial gradients of the design varibale, we have to initialize $\rho$ at nodes and interpolate it to element centroids in the forward pass.
# Volume fraction
vol_frac = 0.4
# SIMP exponent
p = 3.0
# Move limit
move = 0.1
# Regularization parameter
alpha = 0.5
# Minimum thermal conductivity
kappa_min = 1.0
# Compute sensitivities with autograd instead of analytically
TORCH_SENS = True
# evaluate N, N,x and detJ at element centroids
element_centroid = model.etype.ipoints.sum(dim=0) # exploit regularity of the mesh
N, B, detJ = model.eval_shape_functions(element_centroid)
# Initial, minimum, and maximum values of design variables
rho_nodes0 = vol_frac * torch.ones(len(nodes))
rho_min = kappa_min / kappa * torch.ones_like(rho_nodes0)
rho_max = torch.ones_like(rho_nodes0)
# Volume constraint
V_0 = vol_frac * areas.sum()
# Unit-thickness element conductivity matrices for the analytical sensitivity
k0 = torch.einsum("i,ijk->ijk", 1.0 / model.thickness, model.k0())
# Initialize design variable and temperature history
rho_nodes = [rho_nodes0]
temperature = []
# Iterate solutions
for k in tqdm(range(50)):
# Make rho_nodes differentiable for autograd
rho_k_nodes = rho_nodes[k].requires_grad_(True)
# Interpolate rho at element centroids
rho_k = torch.einsum("EN, N -> E", rho_k_nodes[elements], N)
# only for isotropic material -> same effect as updating kappa
model.thickness = rho_k**p
# Compute solution
u_k, f_k, _, _, _ = model.solve(differentiable_parameters=rho_k_nodes)
# Gradient of rho for the Tikhonov regularization
grad_rho = torch.einsum("EiN, EN -> Ei", B, rho_k_nodes[elements])
if TORCH_SENS:
# Sensitivity via automatic differentiation
regularization = torch.einsum("Ei, Ei, E ->", grad_rho, grad_rho, detJ)
objective = torch.inner(f_k.ravel(), u_k.ravel()) + alpha * regularization
sensitivity = torch.autograd.grad(objective, rho_k_nodes)[0]
else:
# Analytical sensitivity dC/drho_E = -p rho_E^(p-1) u.k0.u of the compliance
# and 2 alpha grad(rho).B detJ of the regularization, both chained from the
# element to the nodal design variables
u_e = u_k[elements].reshape(model.n_elem, -1)
w_k = torch.einsum("...i,...ij,...j", u_e, k0, u_e)
dC = torch.einsum("E, N -> EN", -p * rho_k ** (p - 1.0) * w_k, N)
dR = 2.0 * alpha * torch.einsum("Ei, EiN, E -> EN", grad_rho, B, detJ)
sensitivity = torch.zeros_like(rho_k_nodes)
sensitivity.index_add_(0, elements.ravel(), (dC + dR).ravel())
with torch.no_grad():
# For a certain value of mu, apply the iteration scheme
def make_step(mu):
G_k = -sensitivity / mu
# reasoning of OC may be violated here, clamping works to some extent :o
G_k = G_k.clamp(min=0.0)
upper = torch.min(rho_max, (1 + move) * rho_k_nodes)
lower = torch.max(rho_min, (1 - move) * rho_k_nodes)
rho_trial = G_k**0.5 * rho_k_nodes
return torch.maximum(torch.minimum(rho_trial, upper), lower)
# Constraint function
def g(mu):
rho_nodes_it = make_step(mu)
rho_k = torch.einsum("EN, N -> E", rho_nodes_it[elements], N)
return torch.inner(areas, rho_k) - V_0
# Find the root of g(mu)
mu = bisect(g, 0.0, 100000.0)
# Track the design variable and temperature history
rho_nodes.append(make_step(mu))
temperature.append(u_k.detach())
100%|██████████| 50/50 [00:17<00:00, 2.86it/s]
from matplotlib import pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2)
model.plot(
node_property=rho_nodes[-1],
cmap="gray_r",
linewidth=0.0,
bcs=False,
ax=ax1,
title="Topology",
)
model.plot(
node_property=temperature[-1],
cmap="magma",
linewidth=0.0,
bcs=False,
ax=ax2,
title="Temperature",
)
plt.tight_layout()
