Note

Go to the end to download the full example code.

Coupled thermo-mechanical loads

Example of the design of a thermoelastic structure with combined heat and mechanical load

First the heat equations are solved to determine the temperature distribution. After that, the mechanical load due to thermal expansion is calculated based on the temperatures. The compliance of the heat-expansion load combined with the mechanical load is minimized.

This example contains the following specific modules

References:

Gao, T., & Zhang, W. (2010). Topology optimization involving thermo-elastic stress loads. Structural and multidisciplinary optimization, 42, 725-738. DOI: https://doi.org/10.1007/s00158-010-0527-5

 24 import numpy as np
 25 import pymoto as pym
 26
 27 # Problem settings
 28 nx, ny = 60, 80  # Domain size
 29 xmin, filter_radius = 1e-9, 2
 30
 31 load = -100.0  # Mechanical force
 32 heatload = 1.0  # Thermal heat load
 33
 34 volfrac = 0.25
 35
 36
 37 if __name__ == "__main__":
 38     # Set up the domain
 39     domain = pym.VoxelDomain(nx, ny)
 40
 41      # Node and dof groups
 42     nodes_right = domain.nodes[-1, :]
 43     nodes_left = domain.nodes[0, :]
 44
 45     dofs_right = domain.get_dofnumber(nodes_right, [0, 1], ndof=2).flatten()
 46     dofs_left_x = domain.get_dofnumber(nodes_left, 0, ndof=2).flatten()
 47
 48     fixed_dofs = np.unique(np.hstack([dofs_left_x, dofs_right]))
 49
 50     # Setup rhs for loadcase
 51     f = np.zeros(domain.nnodes * 2)  # Generate a force vector
 52     f[2 * domain.nodes[0, 0] + 1] = load
 53     q = np.zeros(domain.nnodes)  # Generate a heat vector
 54     q[domain.nodes[0, 0]] = heatload
 55
 56     # Initial design
 57     s_x = pym.Signal('x', state=volfrac * np.ones(domain.nel))
 58
 59     # Setup optimization problem
 60     with pym.Network() as fn:
 61         # Density filtering
 62         s_xfilt = pym.DensityFilter(domain, radius=filter_radius)(s_x)
 63         s_xfilt.tag = 'Filtered density'
 64
 65         # Plot the design
 66         pym.PlotDomain(domain, saveto="out/design")(s_xfilt)
 67
 68         # RAMP with q = 2
 69         s_xRAMP = pym.MathExpression(f"{xmin} + {1 - xmin}*(inp0 / (3 - 2*inp0))")(s_xfilt)
 70
 71         # Assemble stiffness and conductivity matrix
 72         s_K = pym.AssembleStiffness(domain, bc=fixed_dofs)(s_xRAMP)
 73         s_KT = pym.AssemblePoisson(domain, bc=nodes_right)(s_xRAMP)
 74
 75         # Solve for temperature
 76         s_T = pym.LinSolve()(s_KT, q)
 77         s_T.tag = "temperature"
 78
 79         # Determine thermo-mechanical load
 80         s_Telem = pym.ElementAverage(domain)(s_T)
 81         s_xT = pym.MathExpression("inp0 * inp1")(s_Telem, s_xfilt)
 82         s_thermal_load = pym.ThermoMechanical(domain, alpha=1.0)(s_xT)
 83
 84         # Combine thermo-mechanical and purely mechanical loads
 85         s_load = pym.MathExpression("inp0 + inp1")(f, s_thermal_load)
 86
 87         # Solve mechanical system of equations
 88         s_disp = pym.LinSolve()(s_K, s_load)
 89         s_disp.tag = "displacement"
 90
 91         # Compliance
 92         s_compliance = pym.EinSum('i,i->')(s_disp, s_load)
 93
 94         # Objective
 95         s_objective = pym.Scaling(scaling=10)(s_compliance)
 96         s_objective.tag = "Objective"
 97
 98         # Output to Paraview VTI format
 99         pym.WriteToVTI(domain, saveto="out/dat.vti")(s_xfilt, s_T, s_disp)
100
101         # Volume
102         s_volume = pym.EinSum(expression='i->')(s_xfilt)
103
104         # Volume constraint
105         s_volume_constraint = pym.Scaling(scaling=1, maxval=volfrac * domain.nel)(s_volume)
106         s_volume_constraint.tag = "Volume constraint"
107
108         # Show iteration history
109         pym.PlotIter()(s_objective, s_volume_constraint)
110
111     # Optimization
112     pym.minimize_mma(s_x, [s_objective, s_volume_constraint], verbosity=3, maxit=100)

Gallery generated by Sphinx-Gallery