Gradient-Based Dovetail Joint Shape Optimization for Stiffness

It is common to manufacture an object by decomposing it into parts that can be assembled. This decomposition is often required by size limits of the machine, the complex structure of the shape, etc. To make it possible to easily assemble the final object, it is often desirable to design geometry that enables robust connections between the subcomponents. In this project, we study the task of dovetail-joint shape optimization for stiffness using gradient-based optimization. This optimization requires a differentiable simulator that is capable of modeling the contact between the two parts of a joint, making it possible to reason about the gradient of the stiffness with respect to shape parameters. Our simulation approach uses a penalty method that alternates between optimizing each side of the joint, using the adjoint method to compute gradients. We test our method by optimizing the joint shapes in three different joint shape spaces, and evaluate optimized joint shapes in both simulation and real-world tests. The experiments show that optimized joint shapes achieve higher stiffness, both synthetically and in real-world tests.


INTRODUCTION
The dovetail joint has been used in manufacturing for millennia; Figure 1, for example, is a dovetailed box from an Egyptian tomb Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page.Copyrights for third-party components of this work must be honored.For all other uses, contact the owner/author(s).SCF '23, October 8-10, 2023, New York City, NY, USA  [Ruiz et al. 1984], as shown in Figure 2, are commonly used to connect components in woodworking, turbine blades, 3D printing, etc.In woodworking, there are widely-used heuristics regarding the appropriate angles for hardwood versus softwood [Schwarz and Moskowitz 2009], but how can we determine the optimal dovetail shape under more general conditions and with different materials?
In this work, we study the problem of efficiently finding the optimal dovetail joint design that maximizes the stiffness, given a design space and specified external loads.We choose to study the dovetail joint for its simplicity and ubiquity, but the method we develop should work for any joint shape that is piecewise linear.To find the optimal design, we formalize the task as an optimization problem, with the stiffness as a function of design parameters that characterize the shape of the dovetail joint in the design space.We calculate the gradient of stiffness with respect to the design parameters, and directly optimize them using gradient descent.However, optimizing the shape of a dovetail joint is challenging.First, we need a simulator that can simulate the deformation of the joint given external loads, while considering the contact between two parts of the joint.Second, we need to compute the gradient of stiffness with respect to the design parameters.
To address the first challenge, we perform contact simulations using the penalty approach [Huněk 1993], alternating between the two sides of the joint.In every iteration, we compute the deformation of one side of the joint while considering the other side as rigid, and we penalize the elastic side's penetrating the rigid side.By iteratively applying this alternating approach on the two sides, we have a sufficiently accurate simulator for joint deformation simulation.
As the simulation can be written as partial differential equations (PDEs), the stiffness maximization task can now be understood as a PDE-constrained optimization [De los Reyes 2015].To perform this optimization, we use the adjoint method [Cao et al. 2003;Errico 1997] to calculate the gradients of stiffness with respect to design parameters.
To demonstrate the effectiveness of our method, we test it both in simulation and in real tests.We first verify the accuracy of the gradients by comparing the adjoint method with the finite difference method.We then run optimization on three different dovetail joint design spaces, each using two different initial designs.Experiments show that optimized designs provide greater stiffness compared to the initial ones, both in simulation and in real tests on 3D-printed structures.We finally study the sensitivity of the optimization result with respect to material parameters by optimizing using different Poisson's ratios.Experiments show that the optimization result is not sensitive to different Poisson's ratios.

RELATED WORK 2.1 Dovetail joint shape optimization
Researchers have explored the optimization of dovetail joints using different approaches.The most straightforward approach is simply testing different design parameters.Kogo et al. [2002] and Kogo et al. [2019] conducted tensile and shear tests on carbon-carbon composite dovetail joints with different dovetail angles.Miyauchi et al. [2006] tested wooden dovetail joints with different inclinations and base widths.Jeong et al. [2012] tested different wooden dovetail joints for maximum tension load.Estenlund et al. [2022] studied the dovetail design for mounting coils on rotors by building a simulator and enumerating different dovetail angles.Another approach for dovetail optimization is applying a gradient-free optimizer on a simulator.Hu et al. [2022] tested dovetails with different combinations of tenon length, width, thickness, and angle, and studied the effect of each design parameter using the linear model.Yang et al.
[2018] used commercial FEM software for simulation and optimized dovetail shapes for aero-engines using several different gradientfree optimizers.Some researchers build (differentiable) surrogate models for the simulators and optimize the surrogate models.For example, Hahn and Cofer IV [2012] first optimized design parameters using a surrogate model and further used other gradient-free optimizers to optimize design parameters that are sensitive.In this work, instead of enumerating, using gradient-free optimizers or surrogate models, we directly optimize dovetail design parameters in the FEM simulation, using the adjoint method to compute the gradients.

PDE-constrained optimization
PDE-constrained optimization is a type of constrained optimization problems whose constraints can be written as partial differential equations reflecting the physics that determine the behavior of a system.See, e.g., De los Reyes [2015] for a comprehensive introduction.In a PDE-constrained optimization problem, the objective function depends on both the design variable and the state variable, and the constraints between them can be written as PDEs.There are two major types of approaches to solving PDE constrained optimization problems [Herzog and Kunisch 2010]: black-box and all-at-once.all-at-once methods treat the design variable and the state variable as being independent during optimization, and researchers may use algorithms such as sequential quadratic programming (SQP) [Boggs and Tolle 1995] to solve the problem.black-box methods treat only the design as the independent variable during optimization, and researchers may use the adjoint method to calculate the total derivative of the objective with respect to the design variable.See Givoli [2021] for a tutorial on the adjoint method.In this work, we formalize the dovetail joint shape optimization problem as a PDE-constrained optimization task and use the adjoint method with gradient descent to solve it.

METHOD
The pipeline of our approach is shown in Figure 3.For an initial set of shape parameters, we first calculate the displacement of the corresponding dovetail joint given a fixed external load on the two sides.We simulate this by employing a contact solver with the penalty approach, alternating between sides while assuming the other side is rigidly fixed in the deformed configuration of the previous iteration.After this alternating solver converges, we use the adjoint method to calculate the derivative of the displacement with respect to design parameters and use line search to find an appropriate step size in the direction of steepest descent.

Alternating penalty contact simulator
In this subsection, we describe how we simulate the deformation of a specific joint given external loads.Assuming the material is isotropic, we use the linear elastic model from Langtangen and Logg [2017].Denoting the body as Ω, we have the equations governing the deformation on Ω as where  is the stress tensor,  is the body force (0 in our case),  is the strain tensor, u is the displacement vector,  and  are Lamé parameters, and  is the identity matrix.Under the assumption of planar stress, we have where  is the Young's modulus of the material (we set it to 1 GPa), and  is the Poisson's ratio of the material (we set it to 0.4).Equivalently, we are minimizing the total potential energy Π which is defined as [Alnaes et al. 2015]: where the colon is the dot product between tensors and  is the traction force.For reasonable deformations, we set  to 0.001 GPa for single dovetail joints and 0.003 GPa for complex and double dovetail joints (see Figure 4).We apply equal traction forces on two sides of the joint as Neumann boundary conditions.
One difficulty is to model the contact between the two sides of the joint.As in Bleyer [2018], we use a penalty approach-solving one side while assuming the other side is rigid, and applying the penalty directly on the displacement field u.Denote the two sides of the joint as Ω  and Ω  , and consider the case that we want to solve the deformation on Ω  while considering Ω  rigid.For simplicity and due to the piecewise-linear boundary of the joints, we fit lines to all the contacting edges (see § 4.1 for more details) of Ω  and penalize u on Ω  if it collides with the fitted lines.We have our penalized total potential energy Π  for Ω  as where  pen is the weight of the penalization term, which we set to 1, softplus() = (ln(1 + exp())/) 2 and  is a scale factor that we set to 50, sdf is the signed distance function and sdf(u; Ω  , Ω  ) measures the signed distance from the deformed left-hand side to the fitted lines of the deformed right-hand side (positive if outside and negative if inside).We create meshes with a mesh step size of 0.5 mm using pygmsh [Schlömer 2021] and implement the simulator using FEniCS [Alnaes et al. 2015].We find that four iterations on both sides are usually enough for u to converge.Note that the simulator works for any piecewise linear joints, not only for dovetail joints.

Gradient calculation and optimization
In this subsection, we describe how we optimize the shape design parameters  .With the simulator, we can compute the displacement field u( ).We define the length change of the joint  (u( )) as the difference of average displacement in the horizontal direction between the left and right edges of the joint.Note that maximizing the stiffness is equivalent to minimizing the displacement given fixed traction.We define the optimization objective function L ( ) as where L min_l (•) and L min_w (•) are regularizers, and  min_l and  min_w are weights for the regularizers, both of which we set to 1.The minimum contact length regularizer penalizes contact between edges that are too short: where contact( ) is the set of all contacting edges (see § 4.1 for more details), |•| measures the length of an edge, and we set min_len to 1.5 mm.The minimum width regularization penalty prevents the joint from being too small: where width(•) measures the width of a joint (see § 4.1 for more details), and we set min_width to 3.5 mm.
To calculate the gradient of the objective function with respect to the design parameters, i.e.

dL (𝜃 )
d , we use the adjoint method, which provides the following result: where  (u,  ) = 0 is the PDE corresponding to the simulator.We implement the automatic gradient calculation using dolfin-adjoint [Dokken et al. 2020;Mitusch et al. 2019] and PyTorch [Paszke et al. 2019].
To optimize  , we use gradient descent.For every step, we perform a line search along the gradient direction using SciPy [Virtanen et al. 2020], finding a step size that satisfies strong Wolfe conditions [Wolfe 1969[Wolfe , 1971].If the line search fails, we randomly sample a step size from N (0, 0.5 2 ), where a negative step size indicates moving in the opposite direction.This approach introduces controlled noise to prevent entrapment in local minima.We perform 15 optimization steps and keep the step with minimum  (u( )).
Lastly, to prevent landing in design parameters that are sensitive to manufacturing errors, every time we evaluate the objective function or its gradient, we apply independent random noise sampled from N (0, 0.01 2 ) to every dimension of  three times and take the average.

EXPERIMENTS
In this section, we study the effectiveness of our method.We present different dovetail joint design spaces for all the experiments in § 4.1.We then provide simulated results ( § 4.2) and check the correctness of gradients on mesh vertex coordinates ( § 4.3).We show the main result-the optimization results in § 4.4, evaluating them both synthetically and in real experiments.Finally, to better understand the effect of material parameters, we show the sensitivity test of optimization results with respect to the Poisson's ratio in § 4.5.

Shapes for experiments
As shown in Figure 4, we use three different joint design spaces for all the experiments.We only visualize the left part of the joint, as the right part is complementary to it.All joints are top-to-bottom symmetric, so all simulation is only performed on the lower half of the mesh for lower computational cost.The first design space is named single dovetail joint (Figure 4a), which contains the simplest form of a dovetail joint and three degrees of freedom.The second design space is called complex dovetail joint (Figure 4b), which still contains only one dovetail shape but with a structure that is more complex and has six degrees of freedom.The third design space is double dovetail joint (Figure 4c), which contains two dovetail structures, a non-vertical boundary in the middle, and six degrees of freedom.We annotate all the contacting edges for each design space, which are used in the contact penalizer and the minimum contact length regularizer.We also label the width of every design space, which is used in the minimum width regularize.

Alternating penalty contact simulator results
In this subsection, we visualize simulation results from the alternating penalty contact simulator to demonstrate it is producing reasonable results, and the results are shown in Figure 5.In the first iteration, we undertake the two-step process: first, solving the left-hand side while keeping the right-hand side fixed, leading to a slight leftward movement of the left side.Subsequently, in the same iteration, we address the right-hand side while maintaining the left side in a fixed state, causing the dovetail section on the right side to bend in opposite directions.In the second iteration, we observe a progressive leftward pull of the left-hand side when the right-hand side remains fixed, followed by increased bending of the right-hand side, consequently creating more space for the left-hand side to experience a further leftward pull in the third iteration.With each iteration, the results evolve and eventually converge as the number of iterations increases.

Gradient correctness check
In this subsection, we check the correctness of our gradient computation.We calculate the derivative of the displacement with respect to the coordinates of vertexes from the mesh in two different approaches: the adjoint method and the finite difference method, as shown in Figure 6.For the finite difference method, we use a step size of 10 −4 .We only visualize gradients on three edges, as interior gradients should be all zero.The two sets of gradients are indistinguishable, which infers our gradient calculation is consistent with the finite difference method.The average relative difference from the adjoint method results to the finite difference results is 1.82 × 10 −4 , which is negligibly small.

Shape optimization for stiffness
Fabrication and experimental setup.We use Fusion 3601 to draw 3D shapes and CHITUBOX2 and UVtools3 for slicing.We use the ELEGOO Saturn S resin 3D printer to print laminates of ABS-like resin with a height of 5 mm.For the printed parts to be assemblable, we introduce a gap of 0.1 mm between the joints.To measure the stiffness of a joint, we use an Instron 600DX universal testing machine, which is set to move at a speed of 20 mm/min until the tested object breaks and produces a position-load curve.We also include tabs on two sides of the joint for the machine to clamp, as shown in Figure 7: the parts on the left-hand side of the solid line are the tabs, and the solid lines indicate the position of the clamp.
Optimization and test results.For each design space, we randomly select two different initial sets of shape parameters and optimize them for fifteen gradient descent steps.We select the iteration with the smallest simulated displacement, and the results are shown in Figure 8.The optimized results are similar though the initializations are quite different.The simulated results are shown in Figure 9.We print three copies of each design and test them on the universal testing machine with position-load curves recorded.We measure the stiffness using the position change between a load of 30 N and a load of 60 N, and the results are listed in Table 1.The significant differences between the initial and optimized designs indicate the effectiveness of our algorithm.

Poisson's ratio sensitivity test
One interesting question raised by our method is the effect of material parameters on the optimization results.As Young's modulus is obviously not affecting our optimization process (if simultaneously increasing the traction), we study whether changing Poisson's ratio would change the results, and the results are shown in Figure 10.
As different Poisson's ratios produce almost the same optimization result, we observe that the results are not sensitive to Poisson's ratio.

DISCUSSION
In this project, we studied the task of dovetail joint shape optimization to maximize its stiffness, and, as existing works, we formalized the task as an optimization problem, viewing the stiffness of the joint as a function of shape design parameters.Existing works use search algorithms, gradient-free optimizers, or surrogate models for optimization, which have limited efficiency.To use gradient-based optimizers, we first built our own contact simulator by alternatively simulating the deformation of one side of the joint while considering the other side as rigid, using the penalty approach.We then  use the adjoint method and gradient descent for optimization.For experiments, we first verified the gradients were correct by comparing them to the gradients calculated from the finite difference method.We then tested the optimized joint shapes on different initializations in different dovetail shape design spaces, both synthetically and in real experiments, showing that the optimized joints are much stiffer than the initial ones.Note that our method is not restricted to dovetail joint shape optimization but works for joints with piecewise linear joint boundaries.
We would also like to discuss some limitations and future directions of our approach.First, the simulator has limited accuracy as several assumptions and simplifications are made, e.g., plane stress, fitting boundaries using lines for the penalty term, etc.Real experiments showed that the simulator is still reasonably accurate, but future work can be done on more accurate simulators.Besides, the optimization is not in real-time, and most optimization (15 gradient descent steps) in this project takes 10 to 20 minutes to finish on a laptop computer.However, compared to the manufacturing time, this is acceptable, and there are amortized approaches that

Figure 1 :
Figure 1: Dovetailed box from Hatnefer's tomb, ca.1492-1473 B.C. (Public domain image via The Metropolitan Museum of Art.) Figure 2: Initial and optimized designs of single and complex dovetail joints with average stiffness measured over three batches and external forces applied on the two sides.(a) (c): initial (randomly chosen) single and complex dovetail joints.(b) (d): optimized single and complex dovetail joints, which provide greater stiffness.

Figure 3 :
Figure 3: Given a set of shape parameters, we first use the alternating penalty contact simulator ( § 3.1) to calculate the deformation of the joint, and then use the adjoint method to calculate the gradient and use line search to find the step size for gradient descent ( § 3.2).

Figure 4 :
Figure 4: Three dovetail joint shape design spaces that are used in the experiments.All designs are top-to-bottom symmetric.(a) single dovetail joint is the simplest design; three degrees of freedom.(b) complex dovetail joint with a design that is more complex; six degrees of freedom.(c) double dovetail joint with two dovetails and a non-vertical boundary in the middle; six degrees of freedom.All contact edges and the width are annotated for every design space.
Figure 5: Simulation results on a specific dovetail joint design.The alternating simulator produces reasonable results as the number of iterations increases, and the results converge in the end.

Figure 6 :Figure 7 :
Figure6: Gradient directions on three contacting edges calculated using the adjoint method and the finite difference method.The two results are indistinguishable, which indicates the gradient calculation is correct.

Figure 8 :
Figure 8: Initial and optimized designs of dovetail joints.The optimized designs are similar for the same design space though the initializations are very different.

Figure 9 :
Figure 9: Simulated results and displacements of initial and optimized designs.The displacements of the optimized designs are much smaller than those of the initial designs, which indicates the optimization algorithm is working effectively.

Table 1 :
Real measured stiffness of initial and optimized designs over three batches.The significant increase in stiffness indicates that our algorithm successfully optimizes the joint shape design.