Numerical simulation method for wind tunnel structure under temperature-vibration-wind pressure multiple actions

Abstract

The present application provides a kind of numerical simulation method of wind tunnel body structure under the action of temperature-vibration-wind pressure multiple, belong to wind tunnel technical field, the present application is handled contact boundary condition by using enhanced Lagrangian contact algorithm and setting adaptive penalty parameter and contact stabilization item, modal analysis is carried out to substructure respectively and modal comprehensive coordination matrix is established, the truncated modal residual flexibility correction is introduced and the Pad approximate rational fraction is used to approximate the extrapolation high frequency band dynamic response, when the normal pressure or tangential slip of contact interface exceeds threshold value, the arc length method is used to track the contact state transition path and iteratively update vibration response, the vibration response is input into the contact interface response prediction model based on graph neural network to output stress concentration coefficient and fatigue damage accumulation index, the technical problem of insufficient prediction accuracy of nonlinear response of contact interface under the action of temperature-vibration-wind pressure multiple field coupling of wind tunnel body structure is solved.

Description

Technical Field

[0001] This invention belongs to the field of wind tunnel technology, and more specifically, relates to a numerical simulation method for wind tunnel structure under the multiple effects of temperature, vibration and wind pressure. Background Technology

[0002] As a critical facility for aerospace ground testing, wind tunnel structures endure complex multi-field coupling effects from temperature fields, vibration loads, and wind pressure fluctuations during operation. Traditional numerical simulation methods employ a sequential coupling strategy, first independently calculating temperature stress distribution, then applying it as prestress to vibration analysis, and finally superimposing the wind pressure response. This decoupling method neglects the interaction mechanisms between multiple physics fields. In current wind tunnel structure design, due to the presence of contact interfaces such as bolts and flanges at tunnel connection points, the contact state undergoes strong nonlinear transitions such as opening and closing and slippage under multi-field coupling. Traditional linear superposition methods cannot accurately capture the dynamic response characteristics of the contact interfaces, leading to significant deviations in the prediction of stress concentration locations and fatigue damage accumulation. Existing methods employ the pure penalty function method or the standard Lagrange multiplier method to handle contact nonlinearity. The former suffers from poor convergence and sensitivity to penalty parameters, while the latter has low computational efficiency, making it unsuitable for large-scale engineering problems, and lacks rapid prediction methods for contact interface responses. In other words, existing technologies suffer from insufficient accuracy in predicting the nonlinear response of wind tunnel structures at the contact interfaces under the multi-field coupling effects of temperature, vibration, and wind pressure. Summary of the Invention

[0003] In view of this, the present invention provides a numerical simulation method for wind tunnel structures under multiple effects of temperature, vibration and wind pressure, which can solve the technical problem of insufficient prediction accuracy of nonlinear response of contact interface of wind tunnel structures under multiple field coupling effects of temperature, vibration and wind pressure in the prior art.

[0004] This invention is implemented as follows: It provides a numerical simulation method for wind tunnel structures under multiple effects of temperature, vibration, and wind pressure. This includes collecting historical operating data from similar wind tunnels to establish a multi-field coupled basic database; constructing a three-dimensional finite element model of the wind tunnel and meshing it; setting contact interface elements at the tunnel connection points; decomposing the tunnel structure into main tunnel substructures and connecting segment structures; calculating the thermal stress distribution based on the temperature field distribution data in the multi-field coupled basic database and applying it as the initial stress state; using the enhanced Lagrangian contact algorithm to handle contact boundary conditions; and applying the method to the main tunnel substructures. Modal analysis was performed on the structural components and connecting segments separately. A modal synthesis coordination matrix was established and assembled into the overall structural modal basis. A correction term was introduced to compensate for the contribution of higher-order modes. The Pad approximation rational fraction was used to approximate and extrapolate the high-frequency dynamic response to obtain the full-frequency vibration response time history. When the peak normal pressure or tangential slip at the contact interface exceeds the threshold, the arc length method was used to track the contact state transition path and iteratively update the vibration response. The full-frequency vibration response time history was input into the contact interface response prediction model to output the stress concentration coefficient distribution and fatigue damage accumulation index, generating a multi-field coupling response analysis report of the wind tunnel structure.

[0005] The multi-field coupling basic database includes temperature field distribution data, vibration response data, wind pressure pulsation data, and contact state data of tunnel connection parts.

[0006] Specifically, the process of decomposing the cave structure into a main cave substructure and connecting segment structures involves decomposing it according to stiffness characteristics and refining the connecting segment structures using a high-density mesh.

[0007] The number of modes extracted from the modal analysis of the connecting segment substructure is more than twice that of the main cavity substructure.

[0008] Specifically, the enhanced Lagrange contact algorithm handles contact boundary conditions by setting adaptive penalty parameters and smoothing the friction coefficient, and introducing a contact stabilization term to eliminate high-frequency oscillations at the contact interface.

[0009] The enhanced Lagrange contact algorithm introduces Lagrange multipliers to represent contact pressure at the contact interface, applies contact constraints using a penalty function method, and updates the Lagrange multiplier values ​​in each iteration until the contact conditions are met.

[0010] The adaptive penalty parameter is dynamically adjusted based on the stiffness and penetration of the contact interface unit. A smaller penalty parameter value is used in the initial stage, and the penalty parameter value is gradually increased as the iteration process progresses.

[0011] The friction coefficient smoothing process uses a hyperbolic tangent function to transform the discontinuous friction force-slip velocity relationship of the Coulomb friction model into a smooth continuous function.

[0012] The contact stabilization term is an artificial damping term added to the stiffness matrix of the contact interface element, and the damping coefficient is calculated based on the size and material parameters of the contact interface element.

[0013] Before obtaining the full-band vibration response time history, the wind pressure pulsation data is converted into a frequency domain excitation spectrum, the natural frequencies and mode shapes are extracted, the residual compliance correction coefficient of the truncated mode is calculated, and the wind pressure pulsation frequency domain excitation spectrum is convolved with the corrected frequency response function.

[0014] The modal synthesis coordination matrix establishment process involves introducing interface degrees of freedom at the substructure interfaces, requiring the displacements and forces of each substructure at the interfaces to meet coordination conditions, applying coordination constraints through the Lagrange multiplier method or penalty function method, and assembling the modal coordinate transformation matrices of each substructure into the modal coordinate transformation matrix of the overall structure according to the interface constraint relationship.

[0015] The method for calculating the residual flexibility correction coefficient of the truncated mode is to decompose the total flexibility of the structure into two parts: truncated mode flexibility and residual mode flexibility. The residual mode flexibility is estimated by the difference between the static stiffness matrix of the structure and the truncated mode stiffness matrix. The correction coefficient is the ratio of the residual mode flexibility to the truncated mode flexibility.

[0016] The Pad approximation rational fraction approximation expresses the frequency response function as the ratio of two polynomials, and determines the polynomial coefficients by matching the response values ​​and their derivatives at known frequency points.

[0017] Specifically, the iterative update of the vibration response involves extracting the peak value of the normal pressure and the tangential slip in the vibration response time history, recalculating the contact boundary conditions, until both the peak value of the normal pressure and the tangential slip are lower than the corresponding thresholds.

[0018] The arc length method controls the load step size by introducing an arc length parameter, and automatically adjusts the step size at the extreme points and jump points of the load-displacement curve, so that the load is solved synchronously with the displacement as an unknown quantity.

[0019] The contact pressure threshold is determined based on the yield strength of the material at the connection point of the cavity, and the slip threshold is determined based on the design gap of the connection point. When the contact interface undergoes a state change, a phased loading strategy is adopted, which subdivides the current load step into 5 to 8 sub-steps.

[0020] The contact interface response prediction model is structured as a multi-scale feature extraction architecture based on graph neural networks. The discrete nodes of the contact interface at the connection part of the hole are constructed as a graph structure. The node feature vector includes the node position coordinates, normal pressure, tangential slip, temperature, and historical maximum stress. The edges connect neighboring nodes whose spatial distance is less than 3 times the unit size.

[0021] The graph neural network is designed with three graph convolutional layers. The first graph convolutional layer extracts local features of a single node, the second graph convolutional layer aggregates first-order neighborhood information, and the third graph convolutional layer captures non-local interactions of second-order neighborhoods. A multi-head attention mechanism is introduced between the second and third layers.

[0022] The training of the contact interface response prediction model uses mean square error as the loss function for predicting stress concentration coefficient, and binary cross-entropy as the loss function for judging whether the fatigue damage accumulation index exceeds the safety threshold. The total loss function is a weighted sum of two terms.

[0023] Furthermore, it also includes adjusting and recalculating the structural stiffness parameters of the connection when the stress concentration factor exceeds the safety factor threshold.

[0024] This invention achieves accurate simulation of the multi-field coupled response of a wind tunnel structure by constructing a multi-field coupled database including temperature field, vibration response, wind pressure pulsation, and contact state, combined with substructure modal synthesis technology and an enhanced Lagrange contact algorithm. The enhanced Lagrange contact algorithm is used to handle contact boundary conditions, overcoming the poor convergence of the traditional penalty function method through adaptive penalty parameters and friction coefficient smoothing. The introduced contact stabilization term effectively eliminates high-frequency oscillations at the contact interface. Substructure modal synthesis technology refines the connecting section with a high-density mesh and extracts more modes. Combined with truncated modal residual compliance correction and Pad approximation rational fraction approximation of the extrapolated high-frequency response, the problem of stiffness overestimation caused by modal truncation in traditional methods is solved. A contact interface response prediction model based on a graph neural network is established, and historical multi-condition examples are used to train the model to quickly predict stress concentration factor and fatigue damage accumulation index, overcoming the limitation of low computational efficiency in traditional methods. In summary, this invention solves the technical problem mentioned in the background art of insufficient accuracy in predicting the nonlinear response of the contact interface of the wind tunnel structure under the coupled effects of multiple fields of temperature, vibration, and wind pressure. Attached Figure Description

[0025] Figure 1 This is a flowchart of the method of the present invention.

[0026] Figure 2 This is a frequency domain excitation spectrum distribution diagram of wind pressure pulsation.

[0027] Figure 3 Frequency response diagram of the modal residual compliance correction coefficient.

[0028] Figure 4 The graph shows the convergence curve of the loss function during the training process of the contact interface response prediction model.

[0029] Figure 5 This is a spatial distribution diagram of the stress concentration factor in the flange connection section. Detailed Implementation

[0030] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0031] like Figure 1 The diagram shows a flowchart of a numerical simulation method for a wind tunnel structure under multiple effects of temperature, vibration, and wind pressure, provided by this invention. The method includes the following steps:

[0032] S01. Collect historical operating data of similar wind tunnels to establish a multi-field coupling basic database. The multi-field coupling basic database includes temperature field distribution data, vibration response data, wind pressure pulsation data, and contact state data of tunnel connection parts.

[0033] S02. Construct a three-dimensional finite element model of the wind tunnel and divide it into meshes. Set contact interface elements for the connection parts of the tunnel. Decompose the tunnel structure into main tunnel substructure and connection segment structure according to stiffness characteristics. Refine the connection segment structure with high-density mesh.

[0034] S03. Calculate the thermal stress distribution based on the temperature field distribution data in the multi-field coupling basic database, apply the thermal stress as the initial stress state to the three-dimensional finite element model, use the enhanced Lagrangian contact algorithm to process the contact boundary conditions of the cavity connection, set adaptive penalty parameters and smooth the friction coefficient, and introduce a contact stabilization term to eliminate high-frequency oscillations at the contact interface.

[0035] S04. Convert the wind pressure pulsation data into a frequency domain excitation spectrum, perform modal analysis on the main tunnel substructure and the connecting section substructure to extract natural frequencies and mode shapes. The number of modes extracted from the connecting section substructure is more than twice that of the main tunnel substructure. Establish a modal synthesis coordination matrix to assemble the substructure modes into the overall structural modal basis.

[0036] S05. Calculate the residual compliance correction coefficient of the truncated mode and introduce a correction term to compensate for the contribution of higher-order modes. Use the Pad approximation rational fraction to construct the frequency response function and extrapolate the high-frequency dynamic response. Convolve the wind pressure pulsation frequency domain excitation spectrum with the corrected frequency response function to obtain the full-frequency vibration response time history.

[0037] S06. Extract the peak value of normal pressure and tangential slip at the contact interface in the vibration response time history. When the peak value of normal pressure exceeds the contact pressure threshold or the tangential slip exceeds the slip threshold, use the arc length method to trace the contact state transition path and recalculate the contact boundary conditions. Iterate and update the vibration response until both the peak value of normal pressure and the tangential slip are lower than the corresponding threshold.

[0038] S07. Input the full-band vibration response time history into the contact interface response prediction model. The contact interface response prediction model outputs the stress concentration factor distribution and fatigue damage accumulation index of the connection part of the cavity. When the stress concentration factor exceeds the safety factor threshold, adjust the structural stiffness parameters of the connection part and return to step S02 to recalculate.

[0039] S08. Based on the output results of the full-band vibration response time history and contact interface response prediction model, generate a multi-field coupled response analysis report of the wind tunnel structure. The multi-field coupled response analysis report includes temperature stress distribution cloud map, vibration displacement amplitude distribution cloud map, contact pressure distribution cloud map, and stress concentration danger area identification.

[0040] The enhanced Lagrange contact algorithm is a numerical method for handling contact boundary conditions. It introduces Lagrange multipliers to represent contact pressure at the contact interface and applies contact constraints using a penalty function method. The algorithm updates the Lagrange multiplier values ​​in each iteration until the contact conditions are met. Compared to the pure penalty function method, it can achieve stable convergence with a smaller penalty parameter. The adaptive penalty parameter is dynamically adjusted based on the stiffness and penetration of the contact interface elements. A smaller penalty parameter value is used initially to avoid over-constraint, and the penalty parameter value is gradually increased as the iteration progresses to ensure the accuracy of the contact constraints. The friction coefficient smoothing process uses a hyperbolic tangent function to transform the discontinuous friction force-slip velocity relationship of the Coulomb friction model into a smooth, continuous function, eliminating numerical jumps during the stick-slip transition.

[0041] The contact stabilization term is an artificial damping term added to the stiffness matrix of the contact interface element. The damping coefficient is calculated based on the contact interface element size and material parameters. This term suppresses the high-frequency oscillatory response of the contact interface without affecting the low-frequency overall motion, improving iterative stability by dissipating the non-physical high-frequency energy generated during the contact process. The arc-length method is a solution technique for tracing nonlinear equilibrium paths. By introducing an arc-length parameter to control the load step size, the step size is automatically adjusted at extreme points and jump points of the load-displacement curve to ensure path continuity. This method solves for the load as an unknown quantity simultaneously with the displacement, and is suitable for strongly nonlinear problems with multiple equilibrium branches and contact state transitions.

[0042] The contact pressure threshold is determined based on the yield strength of the material at the connection point of the cavity, and 65% of the yield strength is taken as the threshold. The slippage threshold is determined based on the design gap of the connection point, and 15% of the design gap is taken as the threshold. When the contact interface undergoes a state transition, a staged loading strategy is adopted, subdividing the current load step into 5 to 8 sub-steps, with each sub-step applying 12% to 20% of the total load increment, gradually approaching the actual contact state and avoiding iterative divergence.

[0043] The modal synthesis coordination matrix establishment process is as follows: Interface degrees of freedom are introduced at the substructure interfaces, requiring that the displacements and forces of each substructure at the interface satisfy the coordination conditions. Coordination constraints are applied using the Lagrange multiplier method or the penalty function method, and the modal coordinate transformation matrices of each substructure are assembled into the modal coordinate transformation matrix of the overall structure according to the interface constraint relationships. The assembled modal basis contains both the dominant modes of each substructure and satisfies the interface continuity requirements.

[0044] The method for calculating the truncated modal residual flexibility correction coefficient is as follows: the total structural flexibility is decomposed into truncated modal flexibility and residual modal flexibility. The residual modal flexibility is estimated by the difference between the structural static stiffness matrix and the truncated modal stiffness matrix. The correction coefficient is the ratio of the residual modal flexibility to the truncated modal flexibility. In the dynamic response calculation, this coefficient is multiplied by the truncated modal response to obtain an approximate value of the residual modal contribution. The correction term compensates for the static contribution of higher-order modes to the low-frequency response and reduces the overestimation of stiffness caused by modal truncation.

[0045] The Pad approximation rational fraction approximation is a function approximation technique that expresses the frequency response function as the ratio of two polynomials. By matching the response values ​​and their derivatives at known frequency points to determine the polynomial coefficients, the constructed rational fraction is consistent with the original function within the known frequency band and maintains physical plausibility when extrapolating to higher frequencies. Compared to Taylor series expansion, the Pad approximation maintains high accuracy even far from the expansion point and is suitable for broadband response prediction. When extrapolating the high-frequency dynamic response, the high-frequency components exceeding the modal cutoff frequency in the wind pressure pulsation frequency domain excitation spectrum are convolved with the frequency response function constructed using the Pad approximation rational fraction to obtain the high-frequency vibration response, which is then superimposed with the low-frequency response calculated by the modal superposition method.

[0046] The contact interface response prediction model is structured as a multi-scale feature extraction architecture based on graph neural networks. Discrete nodes at the contact interface of the cavity connection are constructed as a graph structure. Node feature vectors include node position coordinates, normal pressure, tangential slip, temperature, and historical maximum stress. Edges connect neighboring nodes within a spatial distance of less than three times the unit size, and edge feature vectors include inter-node distance and relative displacement. A three-layer graph convolutional layer is designed. The first layer extracts local features of a single node, with a kernel size equal to the node degree and 64 output channels. The second layer aggregates first-order neighborhood information, with 128 output channels. The third layer captures non-local interactions of second-order neighborhoods, with 256 output channels. A multi-head attention mechanism is introduced between the second and third layers, with 8 attention heads adaptively assigning weights to different neighboring nodes. Two fully connected layers follow the graph convolutional layers. The first fully connected layer has 512 neurons, and the second fully connected layer outputs the stress concentration coefficient and fatigue damage accumulation index. The activation function is LeakyReLU with a slope parameter of 0.2.

[0047] The steps for establishing the training dataset for the contact interface response prediction model are as follows: Historical calculation examples containing different temperature conditions, wind pressure pulsation intensities, and geometric parameters of different connection parts are extracted from a multi-field coupled basic database. Each example includes the input features of the contact interface nodes and corresponding stress concentration coefficients and fatigue damage accumulation index labels. The node position coordinates are normalized by dividing the coordinate values ​​by the characteristic length of the cavity; the normal pressure and tangential slip are standardized by subtracting the mean and dividing by the standard deviation; the temperature is normalized and divided by the material melting point temperature. The historical calculation examples are divided into a 75% training set, 15% validation set, and 10% test set. The training set is used for model parameter optimization, the validation set for hyperparameter adjustment and early termination judgment, and the test set for final performance evaluation. The dataset contains at least 30% transitional samples showing changes in contact state to ensure the model learns the features of the state transition process.

[0048] The training steps for the contact interface response prediction model are as follows: Mean squared error is used as the loss function for predicting stress concentration coefficients, and binary cross-entropy is used as the loss function for determining whether the fatigue damage accumulation index exceeds the safety threshold. The total loss function is a weighted sum of two terms with weight coefficients of 0.6 and 0.4, respectively. The Adam optimizer is used, with an initial learning rate of 0.001, which is decayed to 0.5 times the original rate every 30 training epochs. The batch size is set to 32, and the maximum number of training epochs is 200. During training, the validation set loss is monitored, and training is terminated when the validation set loss does not decrease for 15 consecutive epochs. Data augmentation is performed on the training set by adding Gaussian noise with a mean of 0 and a standard deviation of 5% of the feature value to the node features to enhance model robustness. After training, the model performance is evaluated on the test set, requiring an average relative error of less than 8% for stress concentration coefficient prediction and an accuracy of over 92% for determining whether the fatigue damage accumulation index exceeds the safety threshold.

[0049] The stress concentration factor is the ratio of the maximum principal stress at the contact interface node to the nominal stress, where the nominal stress is the average stress of the structure far from the contact area. The stress concentration factor reflects the degree of local stress amplification at the contact interface and is a key parameter for assessing structural strength. The fatigue damage accumulation index is calculated using the rainflow counting method to statistically analyze vibration stress cycles and Miner's linear cumulative damage theory. An index of 1.0 indicates fatigue failure of the structure. The safety factor threshold is determined according to structural design specifications, and is taken as 1.8 to 2.2 for wind tunnel connection areas.

[0050] The variational multiscale turbulent subgrid stress closure method decomposes the velocity field into solvable scale components and subgrid scale components. An orthogonal decomposition is constructed using a variational projection operator to ensure the uniqueness of scale separation. A local boundary value problem is established for the subgrid stress tensor, and an analytical expression including residual basis functions is derived. The residual basis functions characterize small-scale eddy motion features that the mesh cannot analyze. Based on Kolmogorov's energy cascade theory, a scale-adaptive filter width is designed that varies with the local mesh Reynolds number, with the filter width being the geometric mean of the element size and the Kolmogorov scale. The subgrid model coefficients are dynamically adjusted based on the ratio of solvable scale eddy viscosity to total eddy viscosity. The model coefficients are larger in the fully developed turbulent region and close to zero in the laminar or transition regions. Compared to the Smagorinsky model, this method can more accurately predict the Reynolds stress distribution and transition location in separated flows.

[0051] The fast multipole acceleration method for cyclic boundaries is used to handle periodically connected structures. It requires only a finite element model of a single repeating element, applying Floquet-Bloch periodic conditions at the element boundaries to simulate the effect of an infinite periodic array. The Floquet-Bloch conditions require that the physical quantities at the boundaries of adjacent elements satisfy a phase relationship, with the phase difference determined by the Bloch wave vector. The fast multipole method is used to calculate the effect of far-field elements on the current element in a periodic array. The interaction potential function is expanded according to the multipole moment, and the far-field contribution is aggregated through a hierarchical tree structure, reducing the computational complexity from that of traditional methods. Reduce to , The total number of degrees of freedom is represented by . Different Bloch wave vectors are scanned within the Brillouin region, and the corresponding eigenvalue problems are solved to obtain dispersion relation curves and identify the bandgap frequency range. A preconditioned conjugate gradient iterative method is introduced to solve the dense complex coefficient matrix equation generated by the fast multipole method. The preconditioned matrix is ​​constructed using incomplete LU decomposition, accelerating iterative convergence and improving the simulation efficiency of large-scale periodic structures. This method is applicable to the wave propagation characteristic analysis of phononic crystals and metamaterial structures.

[0052] As a preferred embodiment, the present invention also provides a computer-based multi-field coupled numerical simulation system for wind tunnel structures. The computer is equipped with a readable storage medium that stores program instructions. When the program instructions are run on the computer, they execute the above-mentioned numerical simulation method for wind tunnel structures under multiple effects of temperature, vibration, and wind pressure.

[0053] The specific implementation methods of the above steps are described in detail below.

[0054] The specific implementation of step S01 is to first collect multi-field coupling data from the historical operation records of similar wind tunnels and establish a basic database. This process achieves data accumulation through long-term monitoring of the actual operating conditions of the wind tunnel. The collected temperature field distribution data includes the temperature time history curves and spatial distribution characteristics at various locations on the tunnel wall; vibration response data includes the acceleration time-domain signal and spectral characteristics at key locations in the tunnel; wind pressure pulsation data includes the pressure fluctuation amplitude and frequency components at different cross sections inside the tunnel; and contact state data at the tunnel connection points includes the contact pressure distribution and historical slip displacement records. The data acquisition frequency is determined according to the characteristic frequencies of each physical quantity: the temperature field sampling frequency is not less than 0.1Hz, the vibration response sampling frequency is not less than 500Hz, the wind pressure pulsation sampling frequency is not less than 1000Hz, and the contact state monitoring frequency is not less than 100Hz. The established multi-field coupling basic database provides measured benchmark data for subsequent finite element model verification and machine learning model training, ensuring the consistency between simulation results and actual engineering conditions.

[0055] The specific implementation of step S02 involves constructing a three-dimensional finite element model and meshing it based on the geometric dimensions and material properties of the wind tunnel. During the modeling process, the tunnel structure is decomposed into two parts according to stiffness characteristics: the main tunnel substructure and the connecting section substructure. The main tunnel substructure includes parts with high stiffness such as the tunnel wall and supporting frame, while the connecting section substructure includes parts with relatively low stiffness and complex contact behavior, such as flange connections, bolt connections, and sealing structures. The main tunnel substructure is meshed using conventional density hexahedral or tetrahedral elements, with the element feature size controlled within 2% to 5% of the tunnel's feature length. The connecting section structure is refined using high-density meshes, with the element feature size controlled within 0.5% to 1% of the connecting section's feature length to capture stress gradient changes at the contact interface. Contact interface elements are set at the tunnel connection points to simulate contact boundary conditions. The contact interface element type is selected as face-to-face contact elements to adapt to changes in contact state under large deformation conditions. This substructure decomposition and mesh refinement strategy ensures the computational accuracy of the contact area while avoiding the problem of reduced computational efficiency caused by an excessively large number of meshes in the overall model.

[0056] The specific implementation of step S03 is to calculate the thermal stress distribution of the wind tunnel body based on the temperature field distribution data in the multi-field coupling basic database and use it as the initial stress state. During the calculation, the temperature field data is applied as a thermal load to the finite element model. The thermal strain and thermal stress caused by the temperature gradient are calculated through the thermoelastic constitutive relation. The calculation of thermal stress considers the variation characteristics of the material's thermal expansion coefficient with temperature and the influence of constraint conditions on thermal deformation. The calculated thermal stress field is used as the initial stress state and applied to the subsequent vibration response analysis to reflect the influence of temperature on the structural stiffness and stress state. When dealing with the contact boundary conditions at the connection parts of the tunnel body, the enhanced Lagrange contact algorithm is adopted. This algorithm introduces Lagrange multipliers to characterize the contact pressure at the contact interface and combines the penalty function method to apply contact constraints. Compared with the pure penalty function method, it can use a smaller penalty parameter to obtain the results. To achieve stable convergence and avoid over-constraint, the adaptive penalty parameter is dynamically adjusted based on the stiffness and penetration of the contact interface element. In the initial stage, a small penalty parameter value of approximately 0.01 times the stiffness of the contact element is used to avoid over-constraint. As the iteration progresses, the penalty parameter value is gradually increased to 0.1 to 1 times the stiffness of the contact element to ensure the accuracy of the contact constraint. The friction coefficient is smoothed by using a hyperbolic tangent function to transform the discontinuous friction force-slip velocity relationship of the Coulomb friction model into a smooth continuous function, eliminating the numerical jump problem during the stick-slip transition and improving the iteration stability. The introduced contact stabilization term is an artificial damping term added to the stiffness matrix of the contact interface element. The damping coefficient is calculated based on the size and material parameters of the contact interface element. This term suppresses the high-frequency oscillation response of the contact interface by dissipating the non-physical high-frequency energy generated during the contact process, but does not affect the low-frequency overall motion.

[0057] The specific implementation of step S04 involves converting wind pressure pulsation data into a frequency domain excitation spectrum using a fast Fourier transform to obtain the excitation amplitude of each frequency component. Modal analysis is then performed on the main tunnel substructure and the connecting section substructure to extract natural frequencies and mode shapes. During the modal analysis, the eigenvalue problem is solved to obtain the natural frequencies and corresponding mode shape vectors of the structure. The number of modes extracted from the connecting section substructure is more than twice that of the main tunnel substructure because the connecting section structure has more complex local vibration characteristics and requires more modes to describe its dynamic behavior. The process of establishing the modal synthesis coordination matrix involves introducing interface degrees of freedom at the substructure interfaces and requiring the displacements and forces of each substructure at the interfaces to meet coordination conditions. Coordination constraints are applied using the Lagrange multiplier method or penalty function method, and the modal coordinate transformation matrices of each substructure are assembled into the modal coordinate transformation matrix of the overall structure according to the interface constraint relationship. The assembled modal basis contains both the dominant modes of each substructure and meets the interface continuity requirements. Compared with directly performing modal analysis on the overall structure, this substructure modal synthesis method can significantly reduce the computational scale while maintaining sufficient accuracy.

[0058] The specific implementation of step S05 involves calculating the truncated modal residual compliance correction coefficient and introducing a correction term to compensate for higher-order modal contributions to improve the accuracy of vibration response calculation. The calculation method decomposes the total structural compliance into truncated modal compliance and residual modal compliance. The residual modal compliance is estimated by the difference between the structural static stiffness matrix and the truncated modal stiffness matrix. The correction coefficient is the ratio of the residual modal compliance to the truncated modal compliance. In the dynamic response calculation, this coefficient is multiplied by the truncated modal response to obtain an approximate value of the residual modal contribution. The correction term compensates for the static contribution of higher-order modes to the low-frequency response, reducing the overestimation of stiffness caused by modal truncation. The Pad approximation rational fraction approximation is used to construct the frequency response function and extrapolate the high-frequency dynamic response. The Pad approximation expresses the frequency response function... The expression is presented as a ratio of two polynomials. By matching the response values ​​and derivatives at known frequency points, the polynomial coefficients are determined. The constructed rational fraction is consistent with the original function within the known frequency band and maintains physical rationality when extrapolating in the high-frequency band. Compared with the Taylor series expansion, the Pad approximation still has high accuracy when far from the expansion point. When extrapolating the high-frequency dynamic response, the high-frequency components exceeding the modal cutoff frequency in the wind pressure pulsation frequency domain excitation spectrum are convolved with the frequency response function constructed by the Pad approximation rational fraction to obtain the high-frequency vibration response. Finally, the high-frequency vibration response is superimposed with the low-frequency response calculated by the modal superposition method to obtain the full-frequency vibration response time history. This combined method utilizes the efficiency of the modal superposition method in the low-frequency band and ensures the accuracy of the high-frequency response through the Pad approximation.

[0059] The specific implementation of step S06 involves extracting the peak normal pressure and tangential slip at the contact interface from the full-frequency vibration response time history to determine whether the contact state has changed. The contact pressure threshold is determined based on the yield strength of the material at the connection point of the cavity, taking 65% of the yield strength as the threshold. The slip threshold is determined based on the design gap of the connection point, taking 15% of the design gap as the threshold. When the peak normal pressure exceeds the contact pressure threshold or the tangential slip exceeds the slip threshold, it indicates that the contact state may change from adhesion to slip or from contact to separation. At this time, the arc length method is used to trace the contact state change path and recalculate the contact boundary conditions. The arc length method is a solution technique for tracing nonlinear equilibrium paths. By introducing an arc length parameter to control the load step size, the step size is automatically adjusted at the extreme points and jump points of the load-displacement curve to ensure path continuity. This method treats the load as an unknown quantity and solves it synchronously with the displacement. It is suitable for strongly nonlinear problems with multiple equilibrium branches and contact state transitions. When the contact interface undergoes a state transition, a staged loading strategy is adopted to subdivide the current load step into 5 to 8 sub-steps. Each sub-step applies 12% to 20% of the total load increment to gradually approach the true contact state and avoid iterative divergence. The vibration response is iteratively updated until the peak normal pressure and tangential slip are both below the corresponding thresholds, indicating that the contact state is stable. This iterative correction strategy ensures that the influence of contact nonlinearity on the vibration response is accurately considered.

[0060] The specific implementation of step S07 involves inputting the full-band vibration response time history into the contact interface response prediction model to quickly predict the stress concentration coefficient distribution and fatigue damage accumulation index of the tunnel connection. This prediction model is constructed using a multi-scale feature extraction architecture based on graph neural networks. The discrete nodes of the contact interface at the tunnel connection are constructed as a graph structure. The node feature vector contains five types of features: node position coordinates, normal pressure, tangential slip, temperature, and historical maximum stress. Edges connect neighboring nodes with a spatial distance less than three times the unit size, and the edge feature vector contains the distance between nodes and their relative displacement. Three graph convolutional layers are designed to extract features at different scales. The first graph convolutional layer extracts local features of a single node with 64 output channels. The second graph convolutional layer aggregates first-order neighborhood information with 128 output channels. The third graph convolutional layer captures non-local interactions of the second-order neighborhood with 256 output channels. A multi-head attention mechanism with 8 attention heads is introduced between the second and third layers for adaptive allocation. The weights of different neighboring nodes are assigned. After the graph convolutional layer, two fully connected layers are connected. The first layer has 512 neurons. The second layer outputs the stress concentration factor and fatigue damage accumulation index. The activation function uses LeakyReLU with a slope parameter of 0.2 to avoid neuron death. The stress concentration factor is the ratio of the maximum principal stress to the nominal stress at the contact interface node, reflecting the degree of local stress amplification at the contact interface. The fatigue damage accumulation index is calculated using the Miner linear cumulative damage theory based on the rainflow counting method and vibration stress cycles. When the index reaches 1.0, it indicates that the structure has experienced fatigue failure. The safety factor threshold is determined according to the structural design specifications. For the connection parts of the wind tunnel, it is taken as 1.8 to 2.2. When the stress concentration factor exceeds the safety factor threshold, it indicates that the current connection part of the structural design has insufficient strength and the structural stiffness parameters of the connection part need to be adjusted, such as increasing the number of bolts or changing the flange thickness. Then, return to step S02 to rebuild the finite element model for calculation until the stress concentration factor meets the safety requirements.

[0061] The specific implementation of step S08 is to generate a multi-field coupled response analysis report of the wind tunnel structure based on the output results of the full-band vibration response time history and contact interface response prediction model. The report includes a temperature stress distribution cloud map to visually display the distribution characteristics and maximum value location of thermal stress in the tunnel structure, a vibration displacement amplitude distribution cloud map to display the vibration intensity distribution of various parts of the tunnel under wind pressure pulsation excitation, a contact pressure distribution cloud map to display the pressure concentration area and pressure gradient change at the contact interface of the connection part, and a stress concentration danger area marker to highlight the area where the stress concentration coefficient exceeds or is close to the safety threshold, providing a clear target for structural optimization design. During the report generation process, a unified color scale range is used for each cloud map to facilitate the comparison and analysis of results under different working conditions. Dangerous areas are highlighted in red and marked with specific stress concentration coefficient values ​​and fatigue damage accumulation index values. This analysis report provides comprehensive technical support for the safety assessment and optimization design of the wind tunnel structure.

[0062] It should be noted that the key technical ideas of this invention include the following aspects. The first key technology is to use a substructure modal synthesis method combined with modal truncation residual compliance correction and Pad approximation rational fraction extrapolation to achieve efficient calculation of the full-frequency vibration response. Traditional methods directly perform modal analysis on the overall structure, resulting in a large computational scale and stiffness overestimation caused by high-order mode truncation. However, this invention decomposes the structure into main tunnel substructures and connecting section substructures, extracts modes separately, and then assembles them through a coordination matrix, which can significantly reduce the number of degrees of freedom in modal analysis. The introduction of residual compliance correction terms compensates for the static contribution of truncated high-order modes to the low-frequency response. The use of Pad approximation rational fraction to approximate the extrapolation of the high-frequency response avoids the problem of insufficient accuracy in the high-frequency band of traditional methods. This combined strategy significantly reduces the computational cost while ensuring computational accuracy, making broadband dynamic response analysis of large and complex wind tunnel structures feasible. The second key technology is to establish a contact interface response prediction model based on graph neural networks to achieve rapid assessment of stress concentration and fatigue damage. Traditional finite element post-processing methods require calculating stress and fatigue indices for each node one by one, which is computationally intensive and makes it difficult to identify the spatial distribution pattern of stress concentration. In contrast, this invention constructs the contact interface nodes as a graph structure and uses graph convolutional layers to extract multi-scale geometric and physical features. Through a multi-head attention mechanism, it adaptively captures the interaction relationship between different nodes. The trained model can directly output the stress concentration coefficient and fatigue damage index based on the vibration response time history. Compared with the traditional node-by-node calculation method, the prediction speed is increased by tens of times and it can learn the potential laws of stress concentration to provide intelligent decision support for structural optimization. The synergistic effect of these two key technologies is reflected in the fact that the substructure modal synthesis method efficiently obtains the full-frequency vibration response, providing high-quality input data for the graph neural network prediction model, while the graph neural network prediction model quickly identifies stress concentration danger areas, providing feedback guidance for substructure partitioning and mesh refinement strategies. The combination of the two forms a closed-loop analysis process from load excitation to structural response and then to safety assessment, which makes the numerical simulation of wind tunnel structure under the multi-field coupling of temperature, vibration and wind pressure both accurate and efficient, breaking through the bottleneck of traditional methods in which it is difficult to balance computational efficiency and accuracy in complex multi-field coupling problems.

[0063] It should be noted that this invention also solves the following technical problem: Traditional global modal analysis methods suffer from large computational scale and difficulty in capturing local high-order vibration characteristics of connecting sections when dealing with wind tunnel structures with local connections. This invention decomposes the tunnel structure into main tunnel substructures and connecting section substructures according to stiffness characteristics. The connecting section substructures are refined with high-density meshes, and more than twice the number of modes as the main tunnel substructures are extracted. Then, a modal synthesis coordination matrix is ​​established to assemble the substructure modes into a global structural modal basis. This substructure modal synthesis technique significantly reduces the total number of computational degrees of freedom and accurately captures local high-order vibration characteristics of the contact interface by extracting more modes from the connecting section. Simultaneously, introducing interface degrees of freedom at the substructure interfaces and applying coordination constraints ensures that the assembled modal basis meets the interface continuity requirements, thereby improving the prediction accuracy of the vibration response of the connecting section while maintaining computational efficiency.

[0064] Furthermore, this invention addresses the technical problem of underestimation of the response due to the loss of high-order mode contributions when predicting the dynamic response of structures under broadband wind pressure excitation using the modal truncation method. This invention decomposes the total structural flexibility into truncation modal flexibility and residual modal flexibility by calculating the residual flexibility correction coefficient of the truncation mode. The residual modal flexibility is estimated using the difference between the static stiffness matrix and the truncation modal stiffness matrix. A correction term is introduced into the dynamic response calculation to compensate for the static contribution of high-order modes to the low-frequency response. Simultaneously, a Pad approximation rational fraction approximation technique is used to construct the frequency response function. The rational fraction coefficients are determined by matching the response values ​​and derivatives at known frequency points. The high-frequency components exceeding the modal truncation frequency in the wind pressure pulsation frequency domain excitation spectrum are convolved with the extrapolated frequency response function to obtain the high-frequency vibration response. Finally, this is superimposed with the low-frequency response calculated by the modal superposition method to obtain the full-frequency vibration response time history. This strategy of combining residual flexibility correction and high-frequency extrapolation effectively reduces stiffness overestimation caused by modal truncation, ensuring the completeness and accuracy of the dynamic response prediction under broadband excitation.

[0065] Specifically, the principle of this invention is as follows: This invention can solve the technical problem of insufficient prediction accuracy of nonlinear response at contact interfaces. Its principle lies in establishing an asymptotic solution framework involving multi-field coupling of temperature, vibration, and wind pressure. First, the thermal stress caused by the temperature field is applied as the initial stress state, changing the initial contact pressure distribution at the contact interface, making the contact boundary conditions in subsequent vibration analysis more consistent with actual working conditions. The enhanced Lagrange contact algorithm maintains iterative stability while ensuring the accuracy of contact constraints by iteratively updating the Lagrange multiplier values ​​and adaptively adjusting the penalty parameters. Friction coefficient smoothing eliminates numerical jumps during stick-slip transitions, and the contact stabilization term dissipates non-physical high-frequency energy during the contact process. The combination of these three significantly improves the convergence of the contact nonlinear problem. Substructure modal synthesis technology extracts more modes for the connecting substructure, capturing the local high-order vibration characteristics of the contact interface. Truncation mode residual compliance correction compensates for the static contribution of high-order modes to the low-frequency response. Pad approximation rational fraction approximation extrapolates the high-frequency response using known frequency band information. Both of these together reduce modal truncation errors. When the contact state changes, the arc-length method is used to trace the equilibrium path and apply load in stages, avoiding iterative divergence and ensuring the continuity of the state transition process. The graph neural network contact interface response prediction model extracts local and non-local features of the contact interface nodes through multiple graph convolutional layers, learns the mapping relationship between contact state transitions and stress concentrations in historical examples, and achieves fast and accurate response prediction.

[0066] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0067] The specific implementation methods of steps S01, S02 and S08 are the same as those described above, and will not be repeated in detail here.

[0068] The specific implementation of step S03 is as follows: Calculate the thermal stress distribution based on the temperature field distribution data in the multi-field coupling basic database. The formula for the thermal stress distribution is as follows:

[0069] ;

[0070] In the formula, Thermal stress, unit: ; Reference stress, unit: The value is usually 100; This is the elastic modulus, with units of 1. ; For reference elastic modulus, the unit is . The value is usually 200; The coefficient of linear expansion is given by . ; The reference linear expansion coefficient is expressed in units of 1000 ppm. The value is usually taken as ; This is the change in temperature, in units of... ; Poisson's ratio is dimensionless. For reference temperature, the unit is The value is typically taken as 293. The contact constraint equations of the enhanced Lagrange contact algorithm are expressed as follows:

[0071] ;

[0072] In the formula, The normal gap at the contact interface, in units of ; For reference length, the unit is... Take the average size of the contact unit; These are Lagrange multipliers, with units of 1. ; The penalty parameter is in units of The formula for calculating the adaptive penalty parameter is as follows:

[0073] ;

[0074] In the formula, For the first The penalty parameter for the next iteration, in units of ; The initial penalty parameter, in units of experience value ; For the first The penalty parameter for the next iteration, in units of ; The adjustment coefficient is dimensionless and has an empirical value of 0.3. For the first The normal gap of the next iteration, in units of ; These are the characteristic dimensions of the contact element, in units of... The friction coefficient smoothing process uses the hyperbolic tangent function, and the formula is expressed as follows:

[0075] ;

[0076] In the formula, The coefficient of friction for smooth surfaces is dimensionless. The static friction coefficient is dimensionless and has an empirical value of 0.15. Slip velocity, unit: ; For reference speed, the unit is... It is usually 0.01; Let be the hyperbolic tangent function, defined as ,in The base of the natural logarithm. The contact stabilization term is expressed in the stiffness matrix of the contact interface element as follows:

[0077] ;

[0078] In the formula, For stabilizing stiffness, the unit is... ; For reference stiffness, the unit is . The value is usually taken as ; This is the stabilization coefficient, dimensionless, with a default value of 0.05. The equivalent elastic modulus of the contact interface material, in units of ; For reference elastic modulus, the unit is . The value is usually 200; These are the characteristic dimensions of the contact element, in units of... ; The characteristic length of the cave is expressed in units of 1. .

[0079] The interpretation of the enhanced Lagrange contact algorithm involves the Lagrange multiplier update equation, which is expressed as follows:

[0080] ;

[0081] In the formula, For the first The Lagrange multipliers of the next iteration, in units of ; Reference stress, unit: The yield strength of the material is usually taken as the reference value. For the first The Lagrange multipliers of the next iteration, in units of ; The penalty parameter is in units of ; For the first The normal gap of the next iteration, in units of ; For reference length, the unit is... The average size of the contact element is taken. The formula for calculating the damping coefficient of the contact stabilization term is as follows:

[0082] ;

[0083] In the formula, The contact stabilization damping coefficient is expressed in units of 1000 ppm. ; The reference damping coefficient is expressed in units of [missing information]. It usually takes the value 1; The damping ratio is dimensionless and defaults to 0.02. For stabilizing stiffness, the unit is... ; For reference stiffness, the unit is . The value is usually taken as ; Contact unit mass, unit: ; For reference density, the unit is... The value is usually taken as ; The characteristic length of the cave is expressed in units of 1. .

[0084] The specific implementation of step S04 is as follows: The wind pressure fluctuation data is converted into a frequency domain excitation spectrum, and a fast Fourier transform is used to perform modal analysis on the main tunnel substructure and the connecting section substructure to extract natural frequencies and mode shapes. The characteristic equation of the modal analysis is expressed as follows:

[0085] ;

[0086] In the formula, This is the structural stiffness matrix, in units of 1. ; For the first First-order natural angular frequency, in units of ; This is the structural mass matrix, in units of... ; For the first Mode shape vector, dimensionless normalized vector; Let be the modal order. The interface coordination conditions in the process of establishing the modal synthesis coordination matrix are described as follows:

[0087] ;

[0088] In the formula, The number of substructures; For the first The modal matrix of each substructure is dimensionless. For the first The quality matrix of each substructure, in units of ; For the first Modal coordinate vectors of each substructure, in units of ; The overall structural mode matrix is ​​dimensionless. This is the overall structural mass matrix, in units of... ; The overall structural modal coordinate vector, in units of ; For reference quality, the unit is... The total mass of the structure is usually taken as the total mass of the structure. For reference length, the unit is... The characteristic length of the cave is usually taken; Number the substructures.

[0089] The specific implementation of step S05 is as follows: The formula for calculating the residual compliance correction coefficient of the truncated mode is expressed as follows:

[0090] ;

[0091] In the formula, This is the residual compliance correction factor, which is dimensionless. This is the residual modal compliance matrix, in units of ; To truncate the modal compliance matrix, in units of ; This is the static compliance matrix of the structure, in units of . ; To truncate the number of modes; For the first Mode shape vector, dimensionless normalized vector; For the first First-order natural angular frequency, in units of ; Let be the modal order. The frequency response function constructed by the Pad rational fraction approximation is expressed as follows:

[0092] ;

[0093] In the formula, This is the frequency response function, in units of... ; The reference frequency response amplitude is expressed in units of... Take the response amplitude at the cutoff frequency; To excite the angular frequency, the unit is . ; Reference angular frequency, unit: The cutoff frequency is usually chosen; The order of the numerator polynomial is usually taken as 3 to 5; The order of the denominator polynomial is usually taken as 3 to 5; The coefficients of the numerator polynomial are dimensionless. ; The coefficients of the denominator polynomial are dimensionless. ,in ; The index of the numerator polynomial; The index is the denominator polynomial. The convolution operation of the full-frequency vibration response time history is expressed as follows:

[0094] ;

[0095] In the formula, The vibration response displacement time history is given in units of 1. ; The characteristic length of the cave is expressed in units of 1. ; For time, the unit is ; This is the corrected frequency response function, in units of... ; The reference frequency response amplitude is expressed in units of... ; The frequency domain excitation spectrum of wind pressure fluctuations, in units of ; For reference pressure, the unit is... Dynamic pressure is usually taken; The imaginary unit; To excite the angular frequency, the unit is . ; Reference angular frequency, unit: ; is the base of the natural logarithm.

[0096] The specific implementation of step S06 is as follows: extract the peak value of the normal pressure and the tangential slip in the vibration response time history. The contact pressure threshold is determined based on the yield strength of the material at the connection part of the cavity, and the calculation formula is expressed as follows:

[0097] ;

[0098] In the formula, Contact pressure threshold, in units of ; The yield strength of the material, in units of... ; The safety factor is dimensionless and typically takes the value of 1.5. The slip threshold is determined based on the design clearance of the connection part, and the calculation formula is expressed as follows:

[0099] ;

[0100] In the formula, The sliding threshold is dimensionless. For design gaps, the unit is... ; The characteristic length of the cave is expressed in units of 1. The governing equations for tracing the contact state transition path using the arc-length method are expressed as follows:

[0101] ;

[0102] In the formula, This is the arc length increment, in units of ; The characteristic length of the cave is expressed in units of 1. ; This is a displacement increment vector, in units of ; This is a scaling factor, dimensionless, and its empirical value is the ratio of the load to the displacement magnitude. This represents the increment of the load parameters, which is dimensionless.

[0103] The specific implementation of step S07 is as follows: The contact interface response prediction model outputs the stress concentration factor distribution at the connection part of the cavity, and the formula for calculating the stress concentration factor is expressed as follows:

[0104] ;

[0105] In the formula, The stress concentration factor is dimensionless. The maximum principal stress at the contact interface node, in units of ; Nominal stress, unit: ; Reference stress, unit: The fatigue damage accumulation index is typically taken as the material's yield strength. It is calculated using Miner's linear cumulative damage theory, and the formula is as follows:

[0106] ;

[0107] In the formula, The fatigue damage cumulative index is dimensionless. Number of stress cycle types; For the first The number of cycles for each stress amplitude; For the first Fatigue life corresponding to a certain stress amplitude, through material Curve acquisition; This is the number for the stress cycle type. The safety factor threshold is determined according to the structural design specifications, and is taken as 1.8 to 2.2 for the connection parts of the wind tunnel body.

[0108] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0109] The technical team first collected nearly five years of operational data from this wind tunnel and two similar wind tunnels to establish a multi-field coupling basic database. The database contains 86 sets of temperature field distribution data under different operating conditions, covering a temperature range of -15 to 65℃, and records the temperature time history of 217 measuring points on the tunnel surface. Vibration response data comes from accelerometers deployed at eight key locations within the tunnel, with a sampling frequency of 2000Hz and a cumulative recording time exceeding 1200 hours. Wind pressure pulsation data is acquired through 32 dynamic pressure sensors, with a measurement frequency range of 0.1 to 500Hz and a pulsating pressure amplitude range of 0.05 to 8.5 kPa. Contact state data at tunnel connection points includes the contact pressure distribution and relative slippage of 15 flange connection sections, obtained through strain gauges and displacement sensors. After screening and preprocessing, this data forms a basic database containing 1850 valid samples, as shown in Table 1.

[0110] Table 1. Statistical information of the multi-field coupling basic database

[0111]

[0112] The technical team constructed a three-dimensional finite element model of the wind tunnel, comprising four main parts: the main tunnel body, the diffuser section, the contraction section, and the test section. Eight-node hexahedral elements were used for mesh generation, with the mesh size for the main tunnel body controlled at 150mm, resulting in a total of 236,000 elements. The tunnel structure was decomposed according to stiffness characteristics. The main tunnel body substructure includes concrete sections and stable sections, while the connecting section substructure includes 15 flange connection sections and 8 expansion joints. High-density mesh refinement was applied to the connecting section structure, with the mesh size refined to 30mm, and further refined to 8mm in bolt holes and stress concentration areas, resulting in a total of 94,000 connecting section structural elements. Surface-to-surface contact interface elements were set at the flange contact surfaces, with a total contact area of ​​23.6. Each flange face contains an average of 1260 contact points.

[0113] Based on temperature field distribution data from a multi-field coupling database, the technical team calculated the thermal stress distribution under typical high-temperature conditions. The most unfavorable summer condition was selected: an ambient temperature of 38℃, and after 3 hours of operation, the inner wall temperature of the tunnel rose to 56℃. The maximum thermal stress caused by the temperature gradient occurred at the connection between the concrete and steel structure, reaching 17.8 MPa. This thermal stress was applied as the initial stress state to the finite element model, and the enhanced Lagrangian contact algorithm was used to handle the contact boundary conditions of the flange connection surface. The initial value of the adaptive penalty parameter was set to 0.08 times the stiffness of the contact interface element, increasing to 0.35 times with the iteration process. After 15 iterations, the contact penetration converged to 0.012 mm. The friction coefficient was smoothed using a hyperbolic tangent function, smoothly transitioning the discontinuous transition between the static friction coefficient of 0.45 and the dynamic friction coefficient of 0.32. The slip velocity range in the transition zone was set to 0.01 to 0.1 mm / s. A contact stabilization term is introduced into the stiffness matrix of the contact interface element. The damping coefficient is calculated based on the element size of 30 mm and the elastic modulus of steel of 206 GPa. N·s / m, successfully eliminating high-frequency oscillations in the 0.8 to 1.5 kHz frequency band at the contact interface.

[0114] The technical team converted the wind pressure pulsation data into a frequency domain excitation spectrum using a fast Fourier transform. For example... Figure 2 As shown, the wind pressure pulsation energy is mainly concentrated in the 2 to 120 Hz frequency band, with two obvious peaks at 15 Hz and 45 Hz, corresponding to the vortex shedding frequency of the diffusion section and the standing wave frequency of the test section, respectively. Modal analysis was performed on the main tunnel substructure, extracting the first 80 natural frequencies and mode shapes, covering a frequency range of 3.2 to 186 Hz. The first 180 modes were extracted from the connecting section substructure, extending the frequency range to 298 Hz to ensure the capture of the local vibration characteristics of the connecting section. A modal synthesis compatibility matrix was established, introducing 468 interface degrees of freedom at the substructure interfaces. The Lagrange multiplier method was used to apply displacement compatibility conditions, assembling the substructure modes into a global structural modal basis containing 260 modes.

[0115] like Figure 3As shown, the correction coefficients for the truncated modal residual compliance are calculated. The total structural compliance is decomposed into truncated modal compliance and residual modal compliance. The residual modal compliance is estimated using the static stiffness matrix, and the correction coefficients range from 0.12 to 0.28 in the 5-50 Hz frequency band. A correction term is introduced to compensate for the static contribution of higher-order modes to the low-frequency response. The corrected frequency response function shows an 18% increase in amplitude at 20 Hz, effectively reducing the overestimation of stiffness caused by modal truncation. The Pad approximation rational fraction is used to construct the frequency response function and extrapolate the high-frequency dynamic response. The order of the numerator polynomial is set to 12, and the order of the denominator polynomial is set to 10. The polynomial coefficients are determined by matching eight known response points in the 150-200 Hz frequency band. The high-frequency components of the wind pressure pulsation frequency domain excitation spectrum from 200 to 500 Hz are convolved with the Pad approximate rational fraction to obtain the high-frequency vibration response. This high-frequency response is then superimposed with the low-frequency response calculated by the modal superposition method to obtain the full-frequency vibration response time history. The simulation duration is 60 s with a time step of 0.0005 s.

[0116] The peak normal pressure and tangential slip at 15 flange contact interfaces were extracted during the vibration response time history. At flange #3, the peak normal pressure reached 32.5 MPa, exceeding the contact pressure threshold of 32.5 MPa corresponding to 65% of the material's yield strength of 500 MPa. At flange #7, the tangential slip reached 0.048 mm, exceeding the slip threshold of 0.045 mm corresponding to 15% of the design clearance of 0.3 mm. The technical team used the arc-length method to trace the contact state transition path at these two locations, subdividing the current load step into 6 sub-steps, applying 16% of the total load increment to each sub-step, and introducing an arc-length parameter to control the iteration step size range from 0.08 to 0.15. After 8 iterations at flange #3, the peak normal pressure decreased to 31.2 MPa; at flange #7, after 11 iterations, the tangential slip decreased to 0.042 mm, both below the corresponding thresholds, indicating that the contact boundary conditions had returned to equilibrium.

[0117] The technical team input the full-frequency vibration response time history into the contact interface response prediction model. This model, based on a graph neural network architecture, constructs a graph structure from 18,900 contact nodes across 15 flange connection sections. The node feature vector contains five components: node position coordinates, normal pressure, tangential slip, temperature, and historical maximum stress. Edges connect neighboring nodes within a spatial distance of less than 90 mm, generating 156,000 edges, whose feature vectors include inter-node distance and relative displacement. The model comprises three graph convolutional layers: the first layer has 64 output channels, extracting local features of a single node; the second layer has 128 output channels, aggregating first-order neighborhood information; and the third layer has 256 output channels, capturing non-local interactions of second-order neighborhoods. An 8-head attention mechanism is introduced between the second and third layers to adaptively allocate weights to neighboring nodes. Following the graph convolutional layers are two fully connected layers with 512 and 2 neurons respectively, outputting the stress concentration coefficient and fatigue damage accumulation index.

[0118] The training dataset for the contact interface response prediction model extracts 1850 historical examples from a multi-field coupled database, covering different combinations of working conditions: ambient temperature -15 to 65℃, wind speed 80 to 280 m / s, and flange bolt preload 150 to 420 kN. Node position coordinates are normalized by dividing by the tunnel's characteristic length of 45 m; normal pressure and tangential slip are normalized by subtracting the mean of 15.6 MPa and 0.018 mm and then dividing by the standard deviation of 8.2 MPa and 0.013 mm; temperature is normalized by dividing by the material's melting point of 1450℃. The dataset is divided into a 75% training set, 15% validation set, and 10% test set, with 1388 samples in the training set, 278 samples in the validation set, and 184 samples in the test set. Transitional samples representing changes in contact state account for 32% of the dataset, ensuring the model fully learns the characteristics of state transitions.

[0119] Mean squared error was used as the loss function for predicting stress concentration coefficients, and binary cross-entropy was used as the loss function for judging fatigue damage accumulation index. The weight coefficients of the total loss function were 0.6 and 0.4, respectively. The Adam optimizer was used with an initial learning rate of 0.001, which decayed to 0.0005 every 30 training epochs. The batch size was 32, and the maximum number of training epochs was 200. Figure 4 As shown, the validation set loss reached its lowest value of 0.0234 in the 112th epoch during training, and then remained unchanged for 15 consecutive epochs, triggering the early stopping mechanism to terminate training. Gaussian noise with a mean of 0 and a standard deviation of 5% of the eigenvalues ​​was added to the training set for data augmentation. Evaluation on the test set showed that the average relative error of the stress concentration factor prediction was 6.8%, and the accuracy of identifying fatigue damage accumulation index exceeding the safety threshold was 94.3%, both meeting performance requirements.

[0120] The model output shows that the maximum stress concentration factor around the bolt holes of flange #3 is 2.35, exceeding the safety factor threshold of 2.2. The fatigue damage accumulation index at flange #7 is 0.87, and it is expected to reach 1.0 and fail after 8600 hours of operation under the current load conditions. Figure 5 As shown, the stress concentration danger areas are mainly distributed in the annular area with a radius of 30 to 50 mm around the flange bolt holes and at the connection between the outer edge of the flange and the cavity wall. The technical team adjusted the structural stiffness parameters of flanges No. 3 and No. 7, increasing the flange thickness from 38 mm to 46 mm, the bolt diameter from M24 to M27, and the number of bolts from 16 to 20. They then returned to the finite element modeling step, updated the substructure mesh of the connection segment, and recalculated. After three iterations of optimization, the stress concentration factor at all flanges decreased to below 2.0, the maximum fatigue damage accumulation index decreased to 0.63, and the expected fatigue life was extended to 13,800 hours.

[0121] Based on the full-band vibration response time history and contact interface response prediction model outputs, the technical team generated a multi-field coupled response analysis report for the wind tunnel structure. The report includes a temperature stress distribution cloud map, showing the range of temperature stress concentration at the concrete-steel structure connection; a vibration displacement amplitude distribution cloud map, showing that the maximum vibration displacement amplitude at the center of the test section reaches 0.82 mm; a contact pressure distribution cloud map, showing the uniformity of pressure distribution on each flange contact surface; and stress concentration hazard area markers, indicating 12 key areas requiring enhanced monitoring. The report provides a scientific basis for subsequent structural reinforcement and maintenance.

[0122] The technological advancements of this invention compared to traditional wind tunnel structural analysis methods are reflected in multiple aspects. Traditional methods typically apply temperature field, vibration response, and wind pressure load as independent boundary conditions, neglecting the coupling effects between multiple physics fields, leading to significant deviations in the predicted stress state of the contact interface. This invention establishes a multi-field coupling database and introduces thermal stress as the initial stress state into vibration analysis, realistically reflecting the influence of the temperature field on the structural stiffness distribution and contact state. Traditional contact analysis uses a pure penalty function method to handle contact constraints, which is prone to numerical oscillations and convergence difficulties when contact states frequently change. This invention employs an enhanced Lagrangian contact algorithm combined with adaptive penalty parameters and friction coefficient smoothing, ensuring both contact constraint accuracy and improved iterative stability, effectively suppressing high-frequency non-physical oscillations through contact stabilization terms. Traditional modal analysis methods extract overall modes from complex connected structures, resulting in low computational efficiency and difficulty in capturing local vibration characteristics. This invention decomposes the structure into substructures based on stiffness characteristics and assembles the overall modal basis using modal synthesis. It extracts more modal orders from the connecting substructures to ensure accurate characterization of local features. Full-band response prediction is achieved through truncated modal residual compliance correction and Padre approximation rational fraction extrapolation, avoiding errors caused by high-order mode truncation. Traditional fatigue assessment relies on empirical formulas and simplifying assumptions, failing to accurately predict fatigue damage evolution under multi-field coupling. This invention constructs a contact interface response prediction model based on graph neural networks, fully utilizing historical operating data to learn the nonlinear characteristics of the contact state transition process. Through multi-scale feature extraction and an attention mechanism, it adaptively captures the interactions of different neighboring nodes, achieving accurate prediction of stress concentration coefficients and fatigue damage accumulation exponents. These technological improvements make numerical simulation results closer to actual engineering conditions, providing a reliable technical means for wind tunnel structural safety assessment and optimization design.

[0123] It should be noted that the variables involved in this invention are explained in detail in Tables 2 and 3.

[0124] Table 2. Variable Explanation Table (Part 1)

[0125]

[0126] Table 3. Variable Explanation Table (Part Two)

[0127]

[0128] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A numerical simulation method for wind tunnel structure under multiple effects of temperature, vibration, and wind pressure, characterized in that, This process includes collecting historical operational data from similar wind tunnels to establish a multi-field coupling foundation database; constructing a three-dimensional finite element model of the wind tunnel and meshing it; setting contact interface elements at the tunnel connection points; decomposing the tunnel structure into main tunnel substructures and connecting section structures; calculating the thermal stress distribution based on the temperature field distribution data in the multi-field coupling foundation database and applying it as the initial stress state; using the enhanced Lagrange contact algorithm to handle contact boundary conditions; performing modal analysis on the main tunnel substructure and connecting section structures respectively; establishing a modal synthesis coordination matrix and assembling it into the overall structural modal basis; introducing correction terms to compensate for higher-order modal contributions; using the Pad approximation rational fraction to approximate and extrapolate the high-frequency dynamic response to obtain the full-frequency vibration response time history; when the peak normal pressure or tangential slip at the contact interface exceeds a threshold, using the arc length method to track the contact state transition path and iteratively update the vibration response; inputting the full-frequency vibration response time history into the contact interface response prediction model to output the stress concentration coefficient distribution and fatigue damage accumulation index; and generating a multi-field coupling response analysis report for the wind tunnel structure.

2. The method according to claim 1, characterized in that, The multi-field coupling basic database includes temperature field distribution data, vibration response data, wind pressure pulsation data, and contact state data of the tunnel connection parts.

3. The method according to claim 2, characterized in that, The process of decomposing the cave structure into a main cave substructure and connecting segment structures is specifically carried out according to stiffness characteristics, and the connecting segment structures are refined using a high-density mesh.

4. The method according to claim 3, characterized in that, The number of modes extracted from the modal analysis of the connecting segment substructure is more than twice that of the main cavity substructure.

5. The method according to claim 4, characterized in that, The enhanced Lagrange contact algorithm handles contact boundary conditions by setting adaptive penalty parameters and smoothing the friction coefficient, and introducing a contact stabilization term to eliminate high-frequency oscillations at the contact interface.

6. The method according to claim 5, characterized in that, The enhanced Lagrange contact algorithm introduces Lagrange multipliers to represent contact pressure at the contact interface, applies contact constraints in conjunction with the penalty function method, and updates the Lagrange multiplier values ​​in each iteration until the contact conditions are met.

7. The method according to claim 6, characterized in that, The adaptive penalty parameter is dynamically adjusted based on the stiffness and penetration of the contact interface unit. A smaller penalty parameter value is used in the initial stage, and the penalty parameter value is gradually increased as the iteration process progresses.

8. The method according to claim 7, characterized in that, The friction coefficient smoothing process uses a hyperbolic tangent function to transform the discontinuous friction force-slip velocity relationship of the Coulomb friction model into a smooth, continuous function.

9. The method according to claim 8, characterized in that, The contact stabilization term is an artificial damping term added to the stiffness matrix of the contact interface element. The damping coefficient is calculated based on the size and material parameters of the contact interface element.

10. The method according to claim 9, characterized in that, Before obtaining the full-band vibration response time history, the wind pressure pulsation data is converted into a frequency domain excitation spectrum, the natural frequency and mode shape are extracted, the residual compliance correction coefficient of the truncated mode is calculated, and the wind pressure pulsation frequency domain excitation spectrum is convolved with the corrected frequency response function.

Images