Discrete Adjoint Optimization Method for Low-Boom Aircraft Design Using Equivalent Area Distribution

: This paper introduces a low-boom aircraft optimization design method guided by equivalent area distribution, which effectively improves the intuitiveness and refinement of inverse design. A gradient optimization method based on discrete adjoint equations is proposed to achieve the fast solution of the gradient information of target equivalent area distribution relative to design variables and to drive the aerodynamic shape update to the optimal solution. An optimization experiment is carried out based on a self-developed supersonic civil aircraft configuration with engines. The results show that the equivalent area distribution adjoint equation can accurately solve the gradient information. After optimization, the sonic boom level of the aircraft was reduced by 13.2 PLdB, and the drag coefficient was reduced by 60.75 counts. Moreover, the equivalent area distribution adjoint optimization method has outstanding advantages, such as high sensitivity and fast convergence speed, and can take both the low sonic boom and the low drag force of the aircraft into account, providing a powerful tool for the comprehensive optimization design of supersonic civil aircraft by considering sonic boom and aerodynamic force.


Introduction
Modern advanced subsonic airliners have achieved navigation almost on a global scale.However, for long-distance intercontinental flights, commuting is still inefficient.Supersonic civil aircraft has the ability to cruise at high Mach numbers; therefore, it has become the eternal pursuit of aircraft designers.The feasibility demonstration report of the supersonic civil aircraft shows that the passenger traffic of the supersonic civil aircraft can reach more than 25% of the total market operation in the future [1][2][3][4].However, some common problems limit the application scenarios of supersonic civil aircraft.When an aircraft is flying at supersonic speed, the shock wave/expansion wave system generated in the near field propagates and evolves through the atmosphere and finally forms a strong pressure fluctuation with "N" type distribution on the ground, which is the phenomenon of sonic boom [5] characterized by loud and impulsive noise.The earliest research on supersonic civil aircraft began in the 1950s, with the birth of early models such as the Concorde, Tu-144, and B2707.The Concorde also had a glorious operational history from 1976 to 2003.Unfortunately, the initial overpressure of the Concorde's sonic boom reached 2 psf, with a loudness of approximately 107 PLdB [6].It was prohibited by many countries from flying over land and could only demonstrate its advantage of speed over the ocean, restricting economic benefits and ultimately forcing it to be stopped.Therefore, to develop a new generation of supersonic civil aircraft, the first priority is to address the issue of sonic boom reduction.
The configuration and aerodynamic shape of the aircraft to promote low sonic boom has high practical value in engineering applications, which was confirmed by NASA in the SSBD (Shaped Sonic Boom Demonstrator) project [7].For a known aerodynamic shape, parameterizing it and continuously adjusting its details to reduce the level of sonic boom is a positive design method [8][9][10][11].However, this approach largely relies on the initial configuration, and the design requirements are hardly reached due to an unreasonable initial configuration [12].A more scientific and intuitive method is the inverse design method, which starts from a predetermined design goal and reverses the low sonic boom shape [13][14][15][16].The most widely used inverse design method is currently based on the JSGD theory [17][18][19][20][21][22], proposed by Jones, Seebass, George, and Dargen in the 1970s and developed by Rallabhandi [23] and Plotkin [24,25].Lockheed Martin also applied this technology in the "N+2" project.
The principle of the JSGD method is to provide a piecewise function, called the F function, and convert it into an equivalent area distribution through Equation (1).
where Ae is the vector of equivalent area distribution, l is the length of the aircraft, x is the length along the direction of the fuselage, and F(x) is the F function.
Designers continuously adjust the F function by hand in an attempt to obtain the ideal equivalent area distribution and then continue to manually adjust the aerodynamic shape to approach this equivalent area distribution as closely as possible.
The JSGD method has certain limitations, mainly manifested as follows: a This method requires manual adjustment of the aircraft shape to match the target equivalent area distribution, which poses a "Cut and Try" problem.b With this method, it is hard to balance the design requirements of other disciplines, such as aeroelasticity, control and stability, etc. c This method can only be used in the conceptual design stage, making it difficult to meet the accuracy requirements of subsequent detailed design stages.
In recent years, adjoint optimization methods have received increasing attention in the field of low sonic boom research [26][27][28][29].It can not only efficiently optimize the appearance details of large-scale design variables but can also be combined with other disciplinary adjoint equations to carry out a multidisciplinary comprehensive optimization design, which compensates for the shortcomings of the JSGD method.
In order to overcome the limitations of traditional JSGD methods in the detailed design stage, this paper inherits the inverse design concept based on equivalent area distribution and develops a more intuitive and efficient adjoint optimization method.This method takes a reasonable target distribution of the equivalent area as the objective and quickly solves the gradient information on the aerodynamic shape design variables through the adjoint equation.Combined with the optimization algorithm to drive the aerodynamic shape update, it achieves the inverse design from the target equivalent area distribution to the low sonic boom shape.The method of combining CFD calculations with the augmented Burgers equation is used to predict far-field sonic boom signals and evaluate the optimization effect.

Equivalent Area Distribution
The equivalent area distribution is the most critical parameter for determining the sonic boom level of aircraft.In long-term exploration, low sonic boom design guided by the equivalent area distribution has received increasing attention [30][31][32][33][34], which is also the reason why this parameter is chosen as the starting point of this paper.
The disturbances caused by the supersonic aircraft to the surrounding atmosphere propagate in the form of a cone, known as a Mach cone, as shown in Figure 1.Let us consider a fixed observation point on the ground.The dependence domain of this observation point is the cone region with a shape similar to the Mach cone but with an opposite direction.From the moment the aircraft passes through the cone, observers on the ground begin to feel the propagation of disturbances.
As long as the volume and lift disturbances of the aircraft are concentrated at the intersection of the axis and the Mach surface and their effects can be felt at the observation points, the aircraft can be replaced by an equivalent rotational body, and the crosssectional area distribution of the equivalent rotational body is the equivalent area distribution of the aircraft.
There are two approaches for solving equivalent area distribution, i.e., the forward approach and the backward approach.For the reverse design, the backward approach is a more suitable choice [34], which is implemented through the Abel transformation, as shown below: The research work in references [30][31][32][33][34] only takes the equivalent area as a reference to participate in optimization design, and their goal is still to carry the variation of Burgers equation or manually adjust the shape, which is still not intuitive enough.This paper directly carries out a variation on the distribution of equivalent area to provide more direct guidance on shape optimization and derives a discrete form of the adjoint equation, which is easier to implement in applications compared to the continuous adjoint equation.Furthermore, this paper also develops a reliable approach to provide a reasonable equivalent area distribution as the objective function.

Optimization Method
An adjoint optimization method guided by equivalent area distribution was developed.
Starting from ground sonic boom signals that meet airworthiness requirements, the corresponding near-field sonic boom signal of the target is inverted.Furthermore, the target equivalent area distribution is obtained through the Abel transformation, shown as Equation (2).
We introduce the Euclidean distance vector between the current equivalent area distribution and the target equivalent area distribution as the objective function J: where M is the vector dimension of the equivalent area distribution, which is equal to the number of grid points of the near-field sonic boom signal, Ae is the vector of the current Let us consider a fixed observation point on the ground.The dependence domain of this observation point is the cone region with a shape similar to the Mach cone but with an opposite direction.From the moment the aircraft passes through the cone, observers on the ground begin to feel the propagation of disturbances.
As long as the volume and lift disturbances of the aircraft are concentrated at the intersection of the axis and the Mach surface and their effects can be felt at the observation points, the aircraft can be replaced by an equivalent rotational body, and the cross-sectional area distribution of the equivalent rotational body is the equivalent area distribution of the aircraft.
There are two approaches for solving equivalent area distribution, i.e., the forward approach and the backward approach.For the reverse design, the backward approach is a more suitable choice [34], which is implemented through the Abel transformation, as shown below: The research work in references [30][31][32][33][34] only takes the equivalent area as a reference to participate in optimization design, and their goal is still to carry the variation of Burgers equation or manually adjust the shape, which is still not intuitive enough.This paper directly carries out a variation on the distribution of equivalent area to provide more direct guidance on shape optimization and derives a discrete form of the adjoint equation, which is easier to implement in applications compared to the continuous adjoint equation.Furthermore, this paper also develops a reliable approach to provide a reasonable equivalent area distribution as the objective function.

Optimization Method
An adjoint optimization method guided by equivalent area distribution was developed.Starting from ground sonic boom signals that meet airworthiness requirements, the corresponding near-field sonic boom signal of the target is inverted.Furthermore, the target equivalent area distribution is obtained through the Abel transformation, shown as Equation (2).
We introduce the Euclidean distance vector between the current equivalent area distribution and the target equivalent area distribution as the objective function J: where M is the vector dimension of the equivalent area distribution, which is equal to the number of grid points of the near-field sonic boom signal, Ae is the vector of the current equivalent area distribution, and Ae target is the vector of the target equivalent area distribution.
We parameterize the initial shape and solve the gradient information of the objective function relative to the design variables based on the equivalent area distribution adjoint equation and carry out gradient optimization.
The process of this method is shown in Figure 2. It has a typical architecture of gradient optimization.The program has undergone modular design, such as STU for calculating the logical coordinates of parameterized design variables, PMB3D for solving the flow field, BoomPost for near-field overpressure extraction and far-field prediction, Adjoint for solving the adjoint equation, Get-Grad for solving the gradient, RBF-TFI for achieving grid deformation to obtain a new shape, and SLSQP for optimization algorithm.We parameterize the initial shape and solve the gradient information of the objective function relative to the design variables based on the equivalent area distribution adjoint equation and carry out gradient optimization.
The process of this method is shown in Figure 2. It has a typical architecture of gradient optimization.The program has undergone modular design, such as STU for calculating the logical coordinates of parameterized design variables, PMB3D for solving the flow field, BoomPost for near-field overpressure extraction and far-field prediction, Adjoint for solving the adjoint equation, Get-Grad for solving the gradient, RBF-TFI for achieving grid deformation to obtain a new shape, and SLSQP for optimization algorithm.

Proposal of the Objective Function
A reasonable objective function plays a crucial role in the effectiveness of adjoint optimization.The optimal function should be able to form a sonic boom signal with a minimum perceived noise level in the far field.Therefore, inferring the near field from the target far field is a commonly used method for inverse design [3541].Waveform inversion is a key step [4243]; designers should first design a far-field sonic boom signal from a physical perspective and then inverse the corresponding near field to obtain the target equivalent area distribution.
We combine reverse propagation and detailed inverse propagation to achieve inversion from the far field to the near field.

Proposal of the Objective Function
A reasonable objective function plays a crucial role in the effectiveness of adjoint optimization.The optimal function should be able to form a sonic boom signal with a minimum perceived noise level in the far field.Therefore, inferring the near field from the target far field is a commonly used method for inverse design [35][36][37][38][39][40][41].Waveform inversion is a key step [42,43]; designers should first design a far-field sonic boom signal from a physical perspective and then inverse the corresponding near field to obtain the target equivalent area distribution.
We combine reverse propagation and detailed inverse propagation to achieve inversion from the far field to the near field.
The first step is reverse propagation, which is mainly achieved through the reversed Burgers equation, which takes into account the nonlinear effects, dissipation, non-uniform atmosphere, and molecular relaxation phenomena during sound propagation.It should be noted that we have made some simplifications, such as not considering atmospheric turbulence and assuming that there are no buildings on the ground.In fact, these factors will have some impact on the prediction of the ground sonic boom signal [44][45][46].Equation ( 4) gives the dimensionless expression of the augmented Burgers equation [47], where τ is time (dimensionless), σ is the distance along the ray tube (dimensionless), Γ is the gas dissipation parameter (dimensionless), C υ is the relaxation coefficient, θ ν is the relaxation time (dimensionless), S is the ray tube area, C 0 is the sound velocity, and ρ 0 is the atmospheric density; their values are obtained according to atmospheric parameters at the current altitude.A negative sign is taken for the terms related to time in the augmented Burgers equation, thereby achieving the inverse solution of the near-field signal from the far-field signal.Its expression is as follows: We consider the LM1021 standard model [48] as an example, which was a business civil aircraft proposed at the first and second Sonic Boom Prediction Workshop.The nearfield signal obtained by reverse propagation of the far-field signal is shown in Figure 3a.This result reflects the approximate waveform of the real near-field signal.Due to numerical dissipation in reverse propagation, some local shock wave signals are lost, and further inversion of the waveform details is needed.
The signal obtained by reverse propagation is used as the initial waveform, and the overpressure value corresponding to each grid discrete point is defined as a variable.Based on the sonic boom adjoint equation, the gradient information of the target far-field sonic boom signal relative to the design variable is solved, and the variable is updated.During the descent process of the gradient, the initial waveform gradually approaches the real near-field signal in detail.
The expression of the sonic boom adjoint equation is shown as Equation (6).It reflects the sensitivity of far-field sonic boom signals to near-field sonic boom signals.
where k n is the product of the fourth and fifth terms at the right end of the augmented Burgers Equation ( 4), A 1 , B 1 correspond to the relaxation matrix of nitrogen molecules, A 2 , B 2 correspond to the relaxation matrix of oxygen molecules, A 3 , B 3 are absorption process matrices, and γ 0 , γ 1 , β, λ are all adjoint variables.Figure 3a shows the inverse near-field signal.Although it does not completely match the real near-field waveform, Figure 3b shows that their corresponding far-field sonic boom signals are almost identical, indicating that the inverse near-field wave signal is equivalent to the real near-field signal.
To verify the reliability of inversion, we performed frequency domain analysis of far-field signals, presented in Figure 3c,d as the sound pressure level and loudness distribution of corresponding far-field signals in the frequency domain.It can be seen that the detailed inversion effectively restores the high-frequency signal, which represents the shock wave signals.
From Figure 3c, it can be seen that the equivalent area distribution is highly similar; so, the inversion method can be used to calculate the target equivalent area distribution with good reliability.
To verify the reliability of inversion, we performed frequency domain analysis of farfield signals, presented in Figure 3c,d as the sound pressure level and loudness distribution of corresponding far-field signals in the frequency domain.It can be seen that the detailed inversion effectively restores the high-frequency signal, which represents the shock wave signals.This method may have some limitations in applications.When an aircraft adopts a new concept configuration, it may have a very complex near-field shock wave signal, and this method will bring relatively large errors; furthermore, the optimization objectives proposed may not be able to guide optimization effectively.

Equivalent Area Distribution Adjoint Equation
This paper derives the discrete equivalent area distribution adjoint equation suitable for structured solvers as an efficient tool for solving gradient information.The derivation process is as follows.
The basic ideas of discrete adjoint methods are as follows.Optimization problems minimize the objective function: minJ(W, X) (7) where W is the conservative variable of the flow field and X is the design variable of the shape.When the flow field converges, R(W, X) = 0. Introduce the Lagrange multiplier Λ to construct the objective function in the following form: where J is the objective function of the optimization, as shown in Equation (3).When calculation convergence is achieved, R = 0 and L = J.
For the variation of objective function J, the following can be obtained: From Equation ( 9), it can be seen that if a suitable Λ is found to make the first term at the right end equal to 0, as shown in Equation (10), there is no need to calculate the derivative of the flow field relative to the design variable, saving a lot of CFD calculation time and reducing calculation costs.
Equation ( 10) is the adjoint equation.After solving the adjoint equation to obtain Λ, the gradient information of the objective function relative to the design variable can be solved using Equation (11).
where the two terms at the right end are solved using Equations ( 12) and ( 13), respectively.
The goal of deriving the equivalent area distribution adjoint equation is to construct Equation (10).Equation (10) can be transformed into the following form: Construct the left and right-end terms for the adjoint Equation ( 14) accordingly.First, we provide an analysis of the left-end term.This term mainly deals with the calculation residual R of the objective function J.
Figure 4 gives the calculation process of the objective function J. CFD is used to solve the near-field flow field of the aircraft, extract the near-field pressure distribution, and obtain the equivalent area distribution through the Abel transformation.Then, the Euclidean distance from the target equivalent area distribution is calculated.In this process, the calculation residual R of the objective function J only exists in the near-field CFD calculation.Therefore, the left-end term of the equivalent area distribution adjoint equation is consistent with the left-end term of the flow field adjoint equation.Since detailed derivation has been shown in many studies, there is no need for us to rederive it.The optimization platform used in this paper has also established a modular solution program for residual variations of the flow field based on the central scheme and various upwind schemes.
the near-field flow field of the aircraft, extract the near-field pressure distribution, and obtain the equivalent area distribution through the Abel transformation.Then, the Euclidean distance from the target equivalent area distribution is calculated.In this process, the calculation residual R of the objective function J only exists in the near-field CFD calculation.Therefore, the left-end term of the equivalent area distribution adjoint equation is consistent with the left-end term of the flow field adjoint equation.Since detailed derivation has been shown in many studies, there is no need for us to rederive it.The optimization platform used in this paper has also established a modular solution program for residual variations of the flow field based on the central scheme and various upwind schemes.Then, we provided the derivation of the right-end term, which was refined as Equation ( 14).The right-end term can also be performed using an automatic differentiation program.
Substitute Equation ( 3) into the right-end term of Equation ( 14): It is difficult to solve ( ) directly; so, the primitive variable Q is introduced as the intermediate variable, and the chain derivative rule is applied to handle it in the following form: where , the derivatives of the equivalent area distribution ( ) i Ae relative to ρ , u , v , and w are all zero, and only the derivative to p is not zero.Equation ( 16) is transformed as: Substituting Equation (17) into Equation ( 15), we can obtain the following: Then, we provided the derivation of the right-end term, which was refined as Equation ( 14).The right-end term can also be performed using an automatic differentiation program.
Substitute Equation ( 3) into the right-end term of Equation ( 14): It is difficult to solve ∂Ae(i) ∂W k directly; so, the primitive variable Q is introduced as the intermediate variable, and the chain derivative rule is applied to handle it in the following form: where Q = ρ u v w p T , the derivatives of the equivalent area distribution Ae(i) relative to ρ, u, v, and w are all zero, and only the derivative to p is not zero.Equation ( 16) is transformed as: Substituting Equation (17) into Equation ( 15), we can obtain the following: ∂P k ∂W k is described as follows: where u 0 , v 0 , and w 0 are the components of local velocity, V 0 is the modulus of the velocity vector, which means V 0 = u 2 0 + v 2 0 + w 2 0 , and P k is the local pressure.So the goal of establishing Equation ( 18) is the derivation of the term ∂Ae(i) ∂P k .Discretizing Equation (2) leads to the following formula: where C is the constant coefficient: The variation of Equation ( 20) on P k can be obtained as follows: Only if j = k or k + 1, the sum term is not zero, and the following formula is obtained: Further simplification can lead to: Finally, by substituting Equation ( 24) into Equation ( 15), the right-end term of the equivalent area distribution adjoint equation can be derived, and Equation (25) provides the final expression.

Solution and Verification
Using a self-developed business aircraft shape with the low sonic boom feature, shown in Figure 5 as a calculation example, the adjoint equation is solved, with the gradient information verified.

Solution and Verification
Using a self-developed business aircraft shape with the low sonic boom feature, shown in Figure 5 as a calculation example, the adjoint equation is solved, with the gradient information verified.

A. Assembly of Right-End Terms in Parallel Environment
The prerequisite for correctly solving the adjoint equation is the correct assembly of right term in a parallel computing environment.The following part introduces the assembly method of the right-end term of the equivalent area distribution adjoint equation.
The right-end terms should be assembled at the center of the near-field pressure signal extraction grids, and their positions in physical space are determined by a control box.Since the control box spans multiple grid blocks, in a parallel environment, the grid located in the control box is dispersed across multiple processes.The computational grid and control box for the test case are shown in Figure 6.The assembly steps for the right-end item are as follows: (1) After the flow field solution converges, each process indexes and traverses all the physical coordinates of the grid centers within the task using logical coordinates ( ) , , i j k to determine whether they are located in the control box.If so, we record the process number, grid block, and logical coordinates.
(2) The main process collects information recorded in all processes and restores the grid sorting in the real physical space through the bubble sorting method on flow direction.The process number, grid block, and logical coordinates follow the corresponding sorting of physical coordinates.
(3) According to the sorting results in step 2, the output pressure is converted into a near-field overpressure distribution, and the equivalent area distribution is obtained through the Abel transformation.The right-end term is solved according to Equation ( 25) and the right-end term is recombined with the process number, grid block, and logica coordinates sorted in step 2. The file is stored according to the process number.
(4) When solving the adjoint equation, each sub-process traverses all the grid logica coordinates within the task and assembles the right-end item to the corresponding grid center according to the file information output in step 3.
We developed a program based on the above ideas, which is embedded into the optimization platform as a preprocessing module for solving adjoint equations.The form of the right-end term, as solved by this module in this example, is illustrated in Figure 7.The assembly steps for the right-end item are as follows: (1) After the flow field solution converges, each process indexes and traverses all the physical coordinates of the grid centers within the task using logical coordinates (i, j, k) to determine whether they are located in the control box.If so, we record the process number, grid block, and logical coordinates.
(2) The main process collects information recorded in all processes and restores the grid sorting in the real physical space through the bubble sorting method on flow direction.The process number, grid block, and logical coordinates follow the corresponding sorting of physical coordinates.
(3) According to the sorting results in step 2, the output pressure is converted into a near-field overpressure distribution, and the equivalent area distribution is obtained through the Abel transformation.The right-end term is solved according to Equation (25), and the right-end term is recombined with the process number, grid block, and logical coordinates sorted in step 2. The file is stored according to the process number.
(4) When solving the adjoint equation, each sub-process traverses all the grid logical coordinates within the task and assembles the right-end item to the corresponding grid center according to the file information output in step 3.
We developed a program based on the above ideas, which is embedded into the optimization platform as a preprocessing module for solving adjoint equations.The form of the right-end term, as solved by this module in this example, is illustrated in Figure 7.
Due to the inclusion of terms related to the grid scale in the right-end term of the adjoint equation, the gradient solved based on this adjoint equation has a certain sensitivity to the grid scale.To ensure the accuracy of the gradient, the grid inside the control box should be dense and uniform enough.For structural grids, this can lead to an increase in grid scale and special requirements for grid generation.
(4) When solving the adjoint equation, each sub-process traverses all the grid logical coordinates within the task and assembles the right-end item to the corresponding grid center according to the file information output in step 3.
We developed a program based on the above ideas, which is embedded into the optimization platform as a preprocessing module for solving adjoint equations.The form of the right-end term, as solved by this module in this example, is illustrated in Figure 7.

B. Solution of Adjoint Equation and Gradient Verification
The adjoint equation is solved based on the Vanleer flux vector splitting (FVS) scheme, using a mixed scheme of the first order and second order to improve convergence.The mixing coefficient is taken as 0.5, and the time advancement is LU-SGS.
When the CFL number is taken as 50, it results in the residual convergence curves of the five adjoint variables, as shown in Figure 8. Due to the inclusion of terms related to the grid scale in the right-end term of the adjoint equation, the gradient solved based on this adjoint equation has a certain sensitivity to the grid scale.To ensure the accuracy of the gradient, the grid inside the control box should be dense and uniform enough.For structural grids, this can lead to an increase in grid scale and special requirements for grid generation.

B. Solution of Adjoint Equation and Gradient Verification
The adjoint equation is solved based on the Vanleer flux vector splitting (FVS) scheme, using a mixed scheme of the first order and second order to improve convergence.The mixing coefficient is taken as 0.5, and the time advancement is LU-SGS.
When the CFL number is taken as 50, it results in the residual convergence curves of the five adjoint variables, as shown in Figure 8.The contour of the first adjoint variable is shown in Figure 9.The contour of the first adjoint variable is shown in Figure 9.To verify the accuracy of solving gradient information based on the adjoint equation, any seven design variables were chosen to be tested and compared to the gradient information calculated by the central difference method.Table 1 shows the test results.The direction and magnitude of the adjoint gradient of the selected test variables are consistent with the calculation results of the central difference method.Considering the complexity of the shock/expansion wave system in the flow field, there are also certain errors in the gradient information solved by the central difference method.Therefore, the gradient consistency of some test variables is not perfect, but the comparison results are sufficient to prove the reliability of the adjoint method.The contour of the first adjoint variable is shown in Figure 9.To verify the accuracy of solving gradient information based on the adjoint equation, any seven design variables were chosen to be tested and compared to the gradient information calculated by the central difference method.Table 1 shows the test results.The direction and magnitude of the adjoint gradient of the selected test variables are

Optimization Example
As an example, the adjoint optimization method for low-boom aircraft design using equivalent area is applied for the business aircraft shown in Figure 5.

A. Parameterization and Design Variables
The aerodynamic shape parameterization method adopts the Free Form Deformation (FFD) method, with a total of 306 design variables arranged, as shown in Figure 10.Higherscale design variables were appropriately concentrated at key influencing positions, such as the nose, tail, wings, and flat tail.
Aerospace 2024, 11, x FOR PEER REVIEW 13 of consistent with the calculation results of the central difference method.Considering th complexity of the shock/expansion wave system in the flow field, there are also certai errors in the gradient information solved by the central difference method.Therefore, th gradient consistency of some test variables is not perfect, but the comparison results ar sufficient to prove the reliability of the adjoint method.

Optimization Example
As an example, the adjoint optimization method for low-boom aircraft design usin equivalent area is applied for the business aircraft shown in Figure 5.

A. Parameterization and Design Variables
The aerodynamic shape parameterization method adopts the Free Form Deformatio (FFD) method, with a total of 306 design variables arranged, as shown in Figure 1 Higher-scale design variables were appropriately concentrated at key influencin positions, such as the nose, tail, wings, and flat tail.

B. Optimization Model
We selected the design point shown in Equation ( 26) as the optimized desig condition.

B. Optimization Model
We selected the design point shown in Equation ( 26) as the optimized design condition.
A low sonic boom usually brings lift coefficient loss to the aircraft, but it is very limited for shape modification design.Therefore, this paper adopts a fixed angle of attack to save the cost of calculating fixed lift coefficient constraints and achieve faster design.
The optimized mathematical model is shown below.

C. Objective Function
The optimization objective is proposed based on the near-field sonic boom signal inversion method described in Section 4. The calculation process of the objective function is as follows.
Based on the far-field sonic boom signal of the initial aerodynamic shape, we obtain the target far-field sonic boom signal, as shown in Figure 11a.It is designed with multiple weak shock waves at the head and tail, effectively reducing the maximum overpressure value and rising time, thereby weakening the sonic boom intensity.After evaluation, the perceived noise level of the target far-field sonic boom signal is approximately 81.4PLdB.
Based on the near-field sonic boom signal inversion method described in Section 4, the target near-field sonic boom signal corresponding to the target far-field sonic boom signal was obtained, as shown in Figure 11b.
We then calculated the target equivalent area distribution using the Abel transformation shown as Equation (2).The result is shown in Figure 11c.
Finally, we use Equation (3) to solve the objective function.

D. Optimization Result
The optimization is based on the Sequential Quadratic Programming (SQP) algorithm, and after five searches, the objective function is reduced from 3.227 × 10 2 to 9.193 × 10 −1 , as shown in Figure 12.
The equivalent area distribution of the current aerodynamic shape is compared to the initial one, as shown in Figure 13.The equivalent area distribution of the optimized shape approximates the expected value, indicating that solving the gradient information based on the equivalent area distribution adjoint equation can accurately drive the design variables to update to the optimal design point.
The evaluation of aerodynamic and sonic boom performance parameters of the aircraft before and after optimization is shown in Table 2.After optimization, the far-field perceived noise level of the aircraft is reduced by about 13PLdB, indicating that the sonic boom is significantly reduced.
As shown in Figure 14, in the flow field below, the shock waves at the leading edge of the wing and tail are weakened, making a significant contribution to sonic boom reduction.The drag coefficient also decreased by approximately 23 counts.Figure 15 shows that shockwaves in the near field are weakened, which is beneficial for better aerodynamic performance under this flight condition.On the one hand, in the area marked by the blue ellipse, one strong shock wave is converted into three weak shock waves, which is beneficial for reducing drag force.On the other hand, tail shock waves play a crucial role in the formation of a sonic boom.In the area marked by the red ellipse, the tail shock wave weakens and is quickly dissipated in the near field, effectively reducing the impact on the far-field sonic boom.perceived noise level of the target far-field sonic boom signal is approximately 81.4PLdB.
Based on the near-field sonic boom signal inversion method described in Section 4, the target near-field sonic boom signal corresponding to the target far-field sonic boom signal was obtained, as shown in Figure 11b.
We then calculated the target equivalent area distribution using the Abel transformation shown as Equation ( 2).The result is shown in Figure 11c.The equivalent area distribution of the current aerodynamic shape is compared to the initial one, as shown in Figure 13.The equivalent area distribution of the optimized shape approximates the expected value, indicating that solving the gradient information based on the equivalent area distribution adjoint equation can accurately drive the design The equivalent area distribution of the current aerodynamic shape is compared to the initial one, as shown in Figure 13.The equivalent area distribution of the optimized shape approximates the expected value, indicating that solving the gradient information based on the equivalent area distribution adjoint equation can accurately drive the design variables to update to the optimal design point.The evaluation of aerodynamic and sonic boom performance parameters of the aircraft before and after optimization is shown in Table 2.After optimization, the far-field perceived noise level of the aircraft is reduced by about 13PLdB, indicating that the sonic boom is significantly reduced.As shown in Figure 14, in the flow field below, the shock waves at the leading edge of the wing and tail are weakened, making a significant contribution to sonic boom reduction.The drag coefficient also decreased by approximately 23 counts.Figure 15 shows that shockwaves in the near field are weakened, which is beneficial for better aerodynamic performance under this flight condition.On the one hand, in the area marked by the blue ellipse, one strong shock wave is converted into three weak shock waves, which is beneficial for reducing drag force.On the other hand, tail shock waves play a crucial role in the formation of a sonic boom.In the area marked by the red ellipse, the tail shock wave weakens and is quickly dissipated in the near field, effectively reducing the impact on the far-field sonic boom.It reveals an important advantage of the adjoint optimization method using equivalent area distribution.The low sonic boom design mainly focuses on reducing the shock wave below the aircraft, while the low drag design mainly focuses on reducing the shock wave above the aircraft.There is some uncertainty regarding the positive correlation between the two.The traditional sonic boom adjoint equation needs to be combined with the flow field adjoint equation to develop flow field/sonic boom coupling optimization.However, equivalent area distribution not only determines the sonic boom intensity of the aircraft but also whether the aircraft shape conforms to the design criteria of the supersonic area law.The adjoint optimization guided by the equivalent area distribution can balance the sonic boom and drag reduction, opening up a new approach for the comprehensive design of the supersonic civil aircraft.Compare the aerodynamic shape changes of the aircraft and analyze the profiles of the wing and flat tail.In this paper, the profile position is taken at a distance of 1m from the symmetrical plane in the spanwise direction.Figure 16 shows the comparison results.After optimization, the angle of attack of the wing decreases and its curvature decreases, while the curvature of the flat tail increases, the angle of attack increases, and the installation position moves up.The curvature of the fuselage also varies to a certain extent, with the nose drooping and the tail rising, making the fuselage axis more inclined toward a wavy distribution.It reveals an important advantage of the adjoint optimization method using equivalent area distribution.The low sonic boom design mainly focuses on reducing the shock wave below the aircraft, while the low drag design mainly focuses on reducing the shock wave above the aircraft.There is some uncertainty regarding the positive correlation between the two.The traditional sonic boom adjoint equation needs to be combined with the flow field adjoint equation to develop flow field/sonic boom coupling optimization.However, equivalent area distribution not only determines the sonic boom intensity of the aircraft but also whether the aircraft shape conforms to the design criteria of the supersonic area law.The adjoint optimization guided by the equivalent area distribution can balance the sonic boom and drag reduction, opening up a new approach for the comprehensive design of the supersonic civil aircraft.
Compare the aerodynamic shape changes of the aircraft and analyze the profiles of the wing and flat tail.In this paper, the profile position is taken at a distance of 1m from the symmetrical plane in the spanwise direction.Figure 16 shows the comparison results.After optimization, the angle of attack of the wing decreases and its curvature decreases, while the curvature of the flat tail increases, the angle of attack increases, and the installation position moves up.The curvature of the fuselage also varies to a certain extent, with the nose drooping and the tail rising, making the fuselage axis more inclined toward a wavy distribution.Figure 17 shows the pressure distribution on the upper and lower wing surfaces before and after optimization.After optimization, the shock waves on the upper wing surface significantly weakened, while the shock wave positions on the lower wing surface moved backward and toward the wing tip.In addition, a high-pressure zone appeared on the lower wing surface of the flat tail to increase the lift distribution at the flat tail, indicating that for the configuration of the aircraft with an engine, as shown in this paper, designing a higher lift distribution at the flat tail is beneficial for adjusting the equivalent area distribution of the jet region, thereby reducing the level of the sonic boom.Figure 17 shows the pressure distribution on the upper and lower wing surfaces before and after optimization.After optimization, the shock waves on the upper wing surface significantly weakened, while the shock wave positions on the lower wing surface moved backward and toward the wing tip.In addition, a high-pressure zone appeared on the lower wing surface of the flat tail to increase the lift distribution at the flat tail, indicating that for the configuration of the aircraft with an engine, as shown in this paper, designing a higher lift distribution at the flat tail is beneficial for adjusting the equivalent area distribution of the jet region, thereby reducing the level of the sonic boom.

Conclusions
This paper focuses on the core part of the inverse design method for low sonic boom, which is the adjoint optimization guided by equivalent area distribution.The optimization results showed that the equivalent area adjoint optimization method is capable of a fine optimization design of low sonic boom and has the following outstanding advantages: (1) The magnitude of gradient information on equivalent area distribution is relatively large, indicating that the equivalent area distribution adjoint optimization is a highly sensitive design method.Therefore, it can achieve obvious effects in precise design.
(2) After only five generations of line searches, the sonic boom is reduced by 12.6%, indicating that the equivalent area distribution adjoint optimization has a high convergence speed, which can effectively shorten the detailed design cycle for large-scale design variables and improve optimization efficiency.
(3) The equivalent area distribution determines both the sonic boom performance and supersonic drag performance of the aircraft.With a reasonable equivalent area distribution as optimization guidance, it can simultaneously consider the low sonic boom and low drag design of the aircraft, providing a new approach for the comprehensive design of the supersonic civil aircraft.
In summary, the equivalent area distribution adjoint optimization method has achieved a more detailed design of low sonic boom shapes and could provide support for achieving a multidisciplinary comprehensive optimization design of the supersonic civil aircraft.

Figure 4 .
Figure 4. Residual analysis of objective function J.

Figure 4 .
Figure 4. Residual analysis of objective function J.

Figure 5 .
Figure 5. Aerodynamic shape of calculation example ((a) front view, (b) top view, (c) isometric side view, (d) side view).A. Assembly of Right-End Terms in Parallel Environment The prerequisite for correctly solving the adjoint equation is the correct assembly of the right term in a parallel computing environment.The following part introduces the assembly method of the right-end term of the equivalent area distribution adjoint

Figure 11 .
Figure 11.Objective function and calculation process.

Finally
, we use Equation (3) to solve the objective function.

Figure 11 .
Figure 11.Objective function and calculation process.

Figure 13 .
Figure 13.Comparison of equivalent area distribution.

Figure 13 .
Figure 13.Comparison of equivalent area distribution.
(a) Near-field signals (b) Far-field signals

Figure 17 .
Figure 17.Pressure distribution on the lower surface.

Table 2 .
Evaluation of sonic boom and aerodynamic performance.

Table 2 .
Evaluation of sonic boom and aerodynamic performance.
Aerospace 2024, 11, x FOR PEER REVIEW 17 of 21 B 1 relaxation matrix of nitrogen molecules, dimensionless A 2 , B 2 relaxation matrix of oxygen molecules, dimensionless A 3 , B 3 absorption process matrices, dimensionless γ 0 , γ 1 , β, λ matrix of the median adjoint variable during the solving process Subscripts i, j, k index of the grid