Note

Go to the end to download the full example code.

Transient thermal temperature minimization

Minimal example for transient thermal topology optimization

Three heaters are placed in the 2D domain, which are turned on at different times. The optimization problem then minimizes the average temperatures of the heated areas throughout the simulated time.

This example contains some specific modules used in transient problems

References (2D case): - M.J.B. Theulings, R. Maas, L. Noël, F. van Keulen, M. Langelaar

Reducing time and memory requirements in topology optimization of transient problems International Journal for Numerical Methods in Engineering 125(14), 2024. DOI: https://doi.org/10.1002/nme.7461

 23 import numpy as np
 24 import pymoto as pym
 25
 26 nx, ny, nz = 80, 120, 0  # Set nz to zero for the 2D problem, nz > 0 runs a 3D problem
 27 Lx, Ly, Lz = 2/30, 0.1, 0
 28 xmin = 1e-6
 29 filter_radius = 3.0
 30 volfrac = 0.5
 31 # End time [s], as this value decreases the optimal design (of the 2D example) will disconnect from the boundary, as the
 32 # heat has no time to diffuse to the boundary. For 10s, a large connection is formed at the boundary, at 1s a small
 33 # connection, and at 0.1s no connection.
 34 end_time = 0.1
 35 n_steps = 200
 36 dt = end_time / n_steps
 37 theta = 0.5  # time-stepping algorithm, 0.0 for forward Euler, 0.5 for Crank-Nicolson, 1.0 for backward Euler
 38 ndof = 1
 39
 40 if __name__ == "__main__":
 41     print(__doc__)
 42
 43     if nz == 0:  # 2D analysis
 44         # Generate a grid
 45         domain = pym.VoxelDomain(nx, ny, unitx=Lx/nx, unity=Ly/ny)
 46
 47         # Get dof numbers at the boundary
 48         boundary_dofs = domain.nodes[:, 0].flatten()
 49
 50         # Make a force vector of three heaters in the middle of the domain
 51         h_xs = slice(int(nx * (1/2 - 1/15)), int(nx * (1/2 + 1/15)) + 1)
 52         h1_dofs = domain.nodes[h_xs, int(ny * (1/4 - 1/30)) : int(ny * (1/4 + 1/30)) + 1].flatten()
 53         h2_dofs = domain.nodes[h_xs, int(ny * (1/2 - 1/30)) : int(ny * (1/2 + 1/30)) + 1].flatten()
 54         h3_dofs = domain.nodes[h_xs, int(ny * (3/4 - 1/30)) : int(ny * (3/4 + 1/30)) + 1].flatten()
 55         q = np.zeros((domain.nnodes * ndof, int(end_time / dt + 1)))
 56         q[h1_dofs, int(2*n_steps/3):] = 6e3 * dt / h1_dofs.size
 57         q[h2_dofs, int(n_steps/3):] = 7.5e2 * dt / h2_dofs.size
 58         q[h3_dofs, :] = 1e3 * dt / h3_dofs.size
 59         heat_dofs = np.concatenate((h1_dofs, h2_dofs, h3_dofs))
 60
 61     else:
 62         domain = pym.VoxelDomain(nx, ny, nz)
 63         boundary_nodes = domain.nodes[0, :, :].flatten()
 64
 65         boundary_dofs = boundary_nodes
 66         heat_dofs = domain.nodes[nx//2 - nx//8 : nx//2 + nx//8 + 1,
 67                                  ny//2 - ny//8 : ny//2 + ny//8 + 1,
 68                                  nz//2 - nz//8 : nz//2 + nz//8 + 1]
 69
 70         q = np.zeros((domain.nnodes * ndof, int(end_time / dt + 1)))
 71         q[heat_dofs, :] = 1.0  # Uniform heat input of 1.0W at all selected dofs
 72
 73     if domain.nnodes > 1e+6:
 74         print("Too many nodes :(")  # Safety, to prevent overloading the memory in your machine
 75         exit()
 76
 77     if theta == 0.0:
 78         a = 1/1000 # thermal diffusivity = k / (rho*c)
 79         if dt > np.average(domain.element_size)**2 / (2*a):
 80             raise RuntimeError("time step too large for forward Euler, numerical solution unstable")
 81
 82
 83     T_0 = np.zeros(domain.nnodes)
 84
 85     # Make heat and design vector, and fill with initial values
 86     sq = pym.Signal('q', state=q)
 87     sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
 88
 89     # Start building the modular network
 90
 91     # Filter
 92     sxfilt = pym.DensityFilter(domain, radius=filter_radius)(sx)
 93     sxfilt.tag = "Filtered design"
 94     sx_analysis = sxfilt
 95
 96     # Show the design on the screen as it optimizes
 97     if domain.dim == 2:
 98         pym.PlotDomain(domain, saveto="out/design", clim=[0, 1])(sx_analysis)
 99
100     # SIMP material interpolation
101     sSIMP = pym.MathExpression(f"{xmin} + {1.0 - xmin}*inp0^3")(sx_analysis)
102
103     # System matrix assembly module
104     sK = pym.AssemblePoisson(domain, bc=boundary_dofs)(sSIMP)
105     sC = pym.AssembleMass(domain, bc=boundary_dofs, ndof=ndof, material_property=1000)(sSIMP)
106
107     # Transient linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
108     solver = None  # Default (solver is automatically chosen based on matrix properties)
109     # solver = pym.solvers.SolverSparsePardiso()  # Requires Intel MKL installed
110     # solver = pym.solvers.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
111     # solver = pym.solvers.SolverSparseCholeskyScikit()  # Requires scikit installed
112     sT = pym.TransientSolve(dt=dt, end=end_time, x0=T_0, solver=solver)(sq, sK, sC)
113     sT.tag = "temperatures"
114
115     # Output the design, temperature, and heat field to a Paraview file
116     pym.WriteToVTI(domain, saveto='out/dat.vti')(sx_analysis, sT[:, -1], sq[:, -1])
117
118     # Output the design, temperature over time, and heat field over time to Paraview series
119     pym.SeriesToVTI(domain, saveto='out/dat.vti', delta_t=dt, interval=2)(sx_analysis, sT, sq)
120
121     # Total temperature at heat input Tt = sum(T[heat dofs]) for all time steps
122     stemps = pym.EinSum('ij->')(sT[heat_dofs, :])
123     # stemps = pym.EinSum('ij->')(sT)  # Variation: minimize average temperature of the entire domain
124     stemps.tag = 'average temperatures'
125
126     # MMA needs correct scaling of the objective
127     sg0 = pym.Scaling(scaling=100.0)(stemps)
128     sg0.tag = "objective"
129
130     # Calculate the volume of the domain by adding all design densities together
131     svol = pym.EinSum('i->')(sx_analysis)
132     svol.tag = 'volume'
133
134     # Volume constraint
135     sg1 = pym.MathExpression(f'10*(inp0/{domain.nel} - {volfrac})')(svol)
136     sg1.tag = "volume constraint"
137
138     # Maybe you want to check the design-sensitivities?
139     do_finite_difference = False
140     if do_finite_difference:
141         pym.finite_difference(sx, [sg0, sg1], dx=1e-4)
142         exit()
143
144     pym.PlotIter()(sg0, sg1)  # Plot iteration history
145
146     # Do the optimization with MMA
147     pym.minimize_mma(sx, [sg0, sg1])
148
149     # Here you can do some post processing
150     print("The optimization has finished!")
151     print(f"The average temperature value obtained is {np.average(sT.state[heat_dofs, :])}")
152     print(f"The maximum temperature is {np.max(sT.state)}")

Gallery generated by Sphinx-Gallery