Method for quantitatively analyzing dispersion dissipation characteristics of numerical algorithm

Abstract

The invention relates to a method for quantitatively analyzing dispersion and dissipation characteristics of a numerical algorithm. The method comprises the following steps: carrying out numerical discretization on a linearized convection diffusion equation by using the numerical algorithm; solving a transfer matrix and a plurality of transfer characteristic values; sorting according to the mode length of the complex transfer characteristic value corresponding to each mode to obtain a continuous dispersion main mode curve and a continuous dissipation main mode curve; obtaining a dispersion fitting degree; if a continuous dissipation curve cannot be obtained, projecting the class wave solution to each feature vector to obtain a complex weight, and calculating a complex weight to obtain a joint modal dissipation curve; obtaining the maximum fitting range of the dispersion main mode curve and the dissipation main mode curve / joint mode dissipation curve through a 1% rule; and quantitatively analyzing the dispersion dissipation performance based on the dispersion fitting degree and the maximum fitting range. Compared with the prior art, quantitative analysis can be effectively carried out on dispersion and dissipation, and the problem that the main dispersion dissipation mode is difficult to obtain is solved.

Description

Technical Field

[0001] The invention relates to a characteristic analysis method, in particular to a method for quantitatively analyzing the dispersion characteristics of a numerical algorithm. Background Art

[0002] The application background of eigenanalysis in discrete Galerkin (DG) numerical methods is closely related to their in-depth exploration of stability and accuracy in fluid dynamics and other partial differential equations (PDE) problems. DG methods are widely favored for their ability to handle complex geometric configurations and achieve high-order approximations, but their stability is often limited by the characteristics of the numerical format. In this context, eigenanalysis plays a key role: it not only helps researchers gain insight into the stability of the DG method in the process of time advancement, but also predicts the growth trend of the numerical solution through eigenvalue analysis to ensure the reliability of numerical simulation; it also reveals the propagation speed and attenuation characteristics in wave phenomena, which is crucial for evaluating the performance of numerical methods when dealing with fluctuations in different frequency bands. In addition, eigenanalysis can capture transient growth phenomena in abnormal systems, provide a more comprehensive understanding of flow stability, and reveal the significant amplification effect of transient disturbances even in situations where the flow seems stable. By deeply analyzing the eigenvalues, researchers are able to optimize the parameter configuration of the DG method, thereby improving the efficiency and accuracy of the algorithm. In complex problems involving multi-physics coupling, characteristic analysis is also indispensable, as it helps to identify the interactions and stability of different physical processes. In summary, characteristic analysis provides a powerful analytical tool for DG numerical methods, enabling researchers to more deeply understand and optimize the performance of numerical simulations. However, existing characteristic analysis is difficult to obtain the main dispersion and dissipative modes in some cases, and there is currently a lack of a quantitative characteristic analysis method. Summary of the invention

[0003] The purpose of the present invention is to provide a method for quantitatively analyzing the dispersion characteristics of numerical algorithms by combining joint modal analysis, so as to be applied to the performance research of numerical methods such as DG.

[0004] The purpose of the present invention can be achieved by the following technical solutions:

[0005] A method for quantitatively analyzing the dispersion characteristics of a numerical algorithm comprises the following steps:

[0006] Step S1: numerically discretize the linearized convection-diffusion equation using the numerical algorithm to be analyzed;

[0007] Step S2: replacing the values ​​of adjacent cells in the discrete form by a wave-like solution to obtain a transfer matrix;

[0008] Step S3: using the frequency-solving mode with the specified wave number to solve the complex transfer eigenvalue of the corresponding transfer matrix;

[0009] Step S4: Sort by the modulus length of the complex transfer eigenvalue corresponding to each mode, output in positive order in each interval divided equally by the order, obtain a continuous dispersion main modal curve, and obtain a continuous dissipation main modal curve under the condition that the continuity meets the requirements;

[0010] Step S5: Integrate the absolute deviation between the main mode of dispersion and the reference value at low wavenumber to obtain the dispersion fit;

[0011] Step S6: If a continuous dissipation curve cannot be extracted in step S4, the wavelike solution is projected onto each eigenvector to obtain the complex weight;

[0012] Step S7: perform eigendecomposition on the wavelike solution according to the complex weights and then perform multi-integration to obtain the complex weight of each mode in the two-norm energy;

[0013] Step S8: After weighting according to the complex weight, each mode is accumulated to obtain a joint modal dissipation curve;

[0014] Step S9: Obtaining the maximum fitting range between the dispersion main mode curve and the dissipation main mode curve or the joint mode dissipation curve by using the 1% rule;

[0015] Step S10: quantitatively analyzing the performance of the numerical algorithm in dispersion and dissipation based on the dispersion fitting degree and the maximum fitting range, wherein the greater the dispersion fitting degree, the better the dispersion performance, and the larger the maximum fitting range, the better the robustness and accuracy of the numerical algorithm in solving complex flows.

[0016] Among them, steps S1, S2, S3, and S4 are the feature analysis part, which extracts the main mode through the periodicity of the dispersion dissipation curve. Steps S6, S7, and S8 are the joint modal analysis part, which can avoid the difficulty of determining the main mode when discontinuous ghost modes appear. Steps S5 and S9 calculate the quantitative evaluation index from the fit degree and fitting interval with the theoretical value respectively. Step S10 performs analysis based on the quantitative evaluation index obtained in steps S5 and S9.

[0017] In step S1, control parameters of the numerical algorithm are determined according to the numerical algorithm and its application scenario, and the control parameters at least include the polynomial order p, node distribution and the windward coefficient λ of the numerical flux; if the problem to be solved according to the application scenario includes viscosity, the control parameters also include the grid Peclet number Pe and the auxiliary equation numerical flux λ used in calculating the viscosity term a ; When the algorithm includes a stabilization method, the control parameters also include specific parameters corresponding to the stabilization method.

[0018] The step S2 is specifically as follows:

[0019] The values ​​of adjacent cells in the discrete form are replaced by wave-like solutions, that is, u represents any parameter in the solution vector u of the numerical algorithm, I j is the current discretized finite element, i is the imaginary unit, k is the spatial frequency, h represents the actual grid size, and kh represents the wave number. The discrete form is then organized into In this form, the transfer matrix H is obtained.

[0020] The step S3 is specifically as follows: traverse the wave number kh∈[0,(p+1)π], take points at preset intervals, determine a unique transfer matrix according to the value of each kh, and solve the eigenvalue Λ of the corresponding transfer matrix H, where p is the polynomial order.

[0021] The step S4 is specifically as follows:

[0022] According to the eigenvalue Λ corresponding to each mode i i The modulus lengths are sorted and output in positive order in the p+1 equally divided interval on [0,π], that is, the ith wavenumber interval outputs the ith eigenvalue of the current modulus length sorting. When the continuity requirements are met, continuous dispersion main mode curves and dissipation main mode curves are obtained, and the horizontal coordinates of both are standardized wavenumbers. The ordinate of the main mode curve of dispersion is the dispersion rate, which is calculated as follows: The ordinate of the dissipation main mode curve is the dissipation rate, which is calculated as follows: Among them, ω is the corresponding spatial frequency, h represents the actual grid size, ωh is used in combination to reflect the dispersion rate and dissipation rate of the semi-discrete format, k is the spatial frequency, kh represents the wave number, and p is the polynomial order. I corresponds to the real and imaginary parts of the complex number.

[0023] The step S5 is specifically as follows:

[0024] The absolute deviation between the main mode of dispersion and the reference value is integrated at low wavenumbers to obtain the dispersion fit The reference value is given by the equation Given, Pe is the grid Peclet number, and i is the imaginary unit.

[0025] The step S6 is specifically as follows: if a continuous dissipation curve cannot be extracted in step S4, the wave-like solution Project it onto each eigenvector of the transfer matrix H, that is, u=Vα, and use the inverse of the eigenvector matrix V to obtain the complex weight α, where x is the node coordinate of the finite element solution.

[0026] The step S7 is specifically as follows:

[0027] Each component Z of the wavelike solution after eigendecomposition according to the complex weights m =α m v m Perform the double integral, v m is the mth eigenvector, α m is the mth complex weight, and the complex weight of each mode in the two-norm energy is obtained Represents Z m is the mean of , p is the polynomial order, where the integration over the reference unit is done using Gaussian integration.

[0028] The step S8 is specifically as follows:

[0029] For complex weight c m Reweight the real part of Accumulate each mode to obtain That is, the joint modal dissipation curve, where That is, the dissipation value corresponding to the mth eigenvalue in step S4.

[0030] The step S9 is specifically as follows:

[0031] Through the 1% rule, That is, the maximum k value that continuously satisfies the above deviation formula is k 1％ , obtain the maximum fitting range of the dispersion main mode and the dissipation main mode / dissipation joint curve, where, for dispersion, For dissipation,

[0032] Compared with the prior art, the present invention avoids the difficulty of obtaining the main dispersive dissipation mode in some cases by combining modal dissipation. Moreover, the two quantitative analyses adopted in the analysis method can effectively compare the performance of numerical methods and stabilization methods in dispersion and dissipation. Finally, the reference values ​​obtained by this method provide researchers with the possibility of multi-parameter joint optimization of algorithm parameters or stabilization parameters through convex optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 The present invention is a flow chart of the method.

[0034] Figure 2 This is a schematic diagram of the influence of the auxiliary upwind coefficient on the dispersion dissipation characteristics and dispersion fitting degree of the numerical method obtained based on the high-order DG method of the present invention.

[0035] Figure 3 This is the joint modal dissipation curve under different core vector parameter settings obtained by the high-order DG method with exponential filter of the present invention.

[0036] Figure 4 This is the 1% rule maximum fitting interval curve obtained under different artificial viscosity coefficients and core vector parameters based on the high-order DG method with spectral viscosity elimination method of the present invention.

[0037] Figure 5 This is a schematic diagram of the influence of the mixing coefficient on the dispersion dissipation characteristics and the dispersion fitting degree obtained by the high-order DG method based on the positive-preserving limiter of the present invention. DETAILED DESCRIPTION

[0038] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0039] This embodiment provides a method for quantitatively analyzing the dispersion characteristics of numerical algorithms, such as Figure 1 As shown, the following steps are included:

[0040] Step S1: numerically discretize the linearized convection-diffusion equation using the numerical algorithm to be analyzed.

[0041] The high-order finite element method requires the clarification of basic parameters, namely the polynomial order p, node distribution (mainly affecting the core mathematical matrix of the algorithm), and the windward coefficient λ of the numerical flux. If the solution problem determined according to the application scenario includes viscosity, the control parameters also include the grid Peclet number Pe and the auxiliary equation numerical flux λ used to calculate the viscosity term a When the algorithm includes a stabilization method, the control parameters also include specific parameters corresponding to the stabilization method.

[0042] Use the method to be analyzed to numerically discretize the linearized convection (including diffusion) equations, taking the DG method as an example, Where h is the reference value of the discrete unit length, representing the actual grid size, and is usually defined together with the wave number, so there is no need to give a separate value in the calculation. is the mass matrix, is the stiffness matrix, E is the projection matrix from the node to the boundary interpolation, both are determined by the polynomial order p and the node distribution, B = diag (-1, 1), u, f refer to the solution and flux, is the numerical flux. When diffusion (viscosity term) is not considered, f = au. When diffusion is considered, f = au-μq, where a is the convection coefficient, which defaults to 1, μ is the diffusion coefficient, which is determined according to the Peclet number, and q is the first-order derivative of the solution. In DG, it is based on Given, where u+ 、u - Identifies the values ​​on both sides of the finite element boundary, + represents outside the element and - represents inside the element.

[0043] Step S2: The transfer matrix is ​​obtained by replacing the values ​​of adjacent cells in the discrete form with the wavelike solution.

[0044] The values ​​of adjacent cells in the discrete form are replaced by wave-like solutions, that is, u represents any parameter in the solution vector u of the numerical algorithm, I j is the current discretized finite element, i is the imaginary unit, k is the spatial frequency, h represents the actual grid size, and kh represents the wave number. The discrete form is then organized into In this form, the transfer matrix H is obtained.

[0045] Taking the classical DG method to discretize the pure convection equation as an example, For the convection-diffusion equation, For DG including stabilization method, stabilization operator is added on the basis of the above two transfer matrices. In common stabilization method, the filter form is H′=HV q FP q , V q is the matrix for projection from node space to modal space, P q is its inverse, F is the diagonal matrix of the core filter vector with the first term being 1, and each value of which is in the interval [0,1]. The form of spectral viscosity elimination is F SVV It has similar properties to the core filter vector in the filter, but its first term is required to be 0.

[0046] Step S3: Using the frequency-solving mode with the specified wave number, solve the complex transfer eigenvalues ​​of the corresponding transfer matrix.

[0047] Traverse the wave number kh∈[0,(p+1)π], select points with a sufficiently dense preset spacing, determine the unique transfer matrix according to the value of each kh, and solve the eigenvalue Λ of the corresponding transfer matrix H.

[0048] Step S4: Sort by the modulus length of the complex transfer eigenvalue corresponding to each mode, output in positive order in each interval equally divided by the order, obtain a continuous dispersion main modal curve, and obtain a continuous dissipation main modal curve under the condition that the continuity requirements are met.

[0049] According to the eigenvalue Λ corresponding to each mode i iThe modulus lengths are sorted and output in positive order in the p+1 equally divided interval on [0,π], that is, the ith wavenumber interval outputs the ith eigenvalue of the current modulus length sorting. When the continuity requirements are met, continuous dispersion main mode curves and dissipation main mode curves are obtained, and the horizontal coordinates of both are standardized wavenumbers. The ordinate of the main mode curve of dispersion is the dispersion rate, which is calculated as follows: The ordinate of the dissipation main mode curve is the dissipation rate, which is calculated as follows: Among them, ω is the corresponding spatial frequency, h represents the actual grid size, ωh is used in combination to reflect the dispersion rate and dissipation rate of the semi-discrete format, k is the spatial frequency, kh represents the wave number, and p is the polynomial order. I corresponds to the real and imaginary parts of the complex number.

[0050] Step S5: Integrate the absolute deviation of the main mode of dispersion from the reference value at low wavenumbers to obtain the dispersion fit. Integrate the absolute deviation of the main mode of dispersion from the reference value at low wavenumbers to obtain the dispersion fit. The reference value is given by the equation Given, Pe is the grid Peclet number, and i is the imaginary unit.

[0051] Step S6: If a continuous dissipation curve cannot be extracted in step S4, the wavelike solution is projected onto each eigenvector to obtain the complex weight.

[0052] If a continuous dissipation curve cannot be extracted in step S4, the wavelike solution Project it onto each eigenvector of the transfer matrix H, that is, u=Vα, and use the inverse of the eigenvector matrix V to obtain the complex weight α, where x is the node coordinate of the finite element solution.

[0053] Step S7: The wave-like solution is eigen-decomposed according to the complex weights and then multi-integrated to obtain the complex weight of each mode in the two-norm energy.

[0054] Each component Z of the wavelike solution after eigendecomposition according to the complex weights m =α m v m Perform the double integral, v m is the mth eigenvector, α m is the mth complex weight, and the complex weight of each mode in the two-norm energy is obtained Represents Z m where the integration over the reference cell is done using the Gaussian integral, i.e. the mass matrix W.

[0055] Step S8: After weighting according to the complex weight, each mode is accumulated to obtain a joint modal dissipation curve.

[0056] For complex weight c m Reweight the real part of Accumulate each mode to obtain That is, the joint modal dissipation curve, where That is, the dissipation value corresponding to the mth eigenvalue in step S4.

[0057] Step S9: Obtain the maximum fitting range between the dispersion main mode curve and the dissipation main mode curve or the joint mode dissipation curve by using the 1% rule.

[0058] The step S9 is specifically as follows:

[0059] Through the 1% rule, That is, the maximum k value that continuously satisfies the above deviation formula is k 1％ , obtain the maximum fitting range of the dispersion main mode and the dissipation main mode / dissipation joint curve, where, for dispersion, For dissipation,

[0060] When analyzing the fitting range of the dissipation curve of the pure convection equation, since the theoretical dissipation is 0, a large Peclet number (corresponding to extremely small viscosity) can be used to continue to apply the 1% rule.

[0061] Step S10: quantitatively analyzing the performance of the numerical algorithm in dispersion and dissipation based on the dispersion fitting degree and the maximum fitting range, wherein the greater the dispersion fitting degree, the better the dispersion performance, and the larger the maximum fitting range, the better the robustness and accuracy of the numerical algorithm in solving complex flows.

[0062] When certain parameters of the algorithm need to be compared, S2-S9 can be repeated for the algorithms under different parameters to obtain the dispersion dissipation characteristic analysis under different parameters.

[0063] Based on the high-order DG method, the influence of the auxiliary upwind coefficient on the dispersion characteristics of the numerical method and the dispersion fitting degree is obtained, such as Figure 2 As shown in the figure, the results show that when the Peclet number is large, that is, the viscosity is small, the absolute value of the headwind coefficient of the viscous term flux |λ a |The dispersion is smaller when it is close to 1. When the Peclet number is small, that is, the viscosity is large, the absolute value of the windward coefficient of the viscous term flux |λ a |The dispersion is smaller when it is close to 0.

[0064] High-order DG based on exponential filter, namely core vector The joint modal dissipation curves obtained under different core vector parameter η and s settings are as follows: Figure 3As shown in Figure 2, the results show that when η remains unchanged and s changes, the excessive dissipation of the numerical method at low wavenumbers will be quite obvious, while when η changes and s remains unchanged, the dissipation of the numerical method can be appropriately adjusted while avoiding significant excessive dissipation. That is, for the filtering method, η is more suitable to be adjusted as an adaptive value.

[0065] Based on the high-order DG method with spectral viscosity elimination, the maximum fitting interval curve of the 1% rule under different artificial viscosity coefficients and kernel vector parameters is obtained, such as Figure 4 As shown in the figure, the results show that no matter what kind of kernel vector, the spectral viscosity elimination method can quickly converge to the maximum dispersion wave number, but cannot converge to the maximum dissipation wave number, that is, there is an additional dissipation that is difficult to suppress.

[0066] Based on the high-order DG method of mixing low-order positive-preserving limiters, the influence of the obtained mixing coefficient on the dispersion dissipation characteristics and the dispersion fitting degree, such as Figure 5 As shown, the results show that for high-order methods, appropriately mixing a very small amount of low-order methods can effectively improve the algorithm dispersion performance.

[0067] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning, or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A method for quantitatively analyzing the dispersion characteristics of numerical algorithms, characterized in that: The following steps are involved: Step S1: numerically discretize the linearized convection-diffusion equation using the numerical algorithm to be analyzed; Step S2: replacing the values ​​of adjacent cells in the discrete form by a wave-like solution to obtain a transfer matrix; Step S3: using the specified wave number to solve the frequency mode, solving the complex transfer eigenvalue of the corresponding transfer matrix; Step S4: Sort by the modulus length of the complex transfer eigenvalue corresponding to each mode, output in positive order in each interval divided equally by the order, obtain a continuous dispersion main modal curve, and obtain a continuous dissipation main modal curve under the condition that the continuity meets the requirements; Step S5: Integrate the absolute deviation between the main mode of dispersion and the reference value at low wavenumber to obtain the dispersion fit; Step S6: If a continuous dissipation curve cannot be extracted in step S4, the wavelike solution is projected onto each eigenvector to obtain the complex weight; Step S7: perform eigendecomposition on the wavelike solution according to the complex weights and then perform multi-integration to obtain the complex weight of each mode in the two-norm energy; Step S8: After weighting according to the complex weight, each mode is accumulated to obtain a joint modal dissipation curve; Step S9: Obtaining the maximum fitting range between the dispersion main mode curve and the dissipation main mode curve or the joint mode dissipation curve by using the 1% rule; Step S10: quantitatively analyzing the performance of the numerical algorithm in dispersion and dissipation based on the dispersion fitting degree and the maximum fitting range, wherein the greater the dispersion fitting degree, the better the dispersion performance, and the larger the maximum fitting range, the better the robustness and accuracy of the numerical algorithm in solving complex flows.

2. A method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 1, characterized in that: In step S1, control parameters of the numerical algorithm are determined according to the numerical algorithm and its application scenario, and the control parameters at least include the polynomial order p, node distribution and the windward coefficient λ of the numerical flux; if the problem to be solved according to the application scenario includes viscosity, the control parameters also include the grid Peclet number Pe and the auxiliary equation numerical flux λ used in calculating the viscosity term a ; When the algorithm includes a stabilization method, the control parameters also include specific parameters corresponding to the stabilization method.

3. The method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 1, characterized in that: The step S2 is specifically as follows: The values ​​of adjacent cells in the discrete form are replaced by wave-like solutions, that is, u represents any parameter in the solution vector u of the numerical algorithm, I j is the current discretized finite element, i is the imaginary unit, k is the spatial frequency, h represents the actual grid size, and kh represents the wave number. The discrete form is then organized into In this form, the transfer matrix H is obtained.

4. The method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 1, characterized in that: The step S3 is specifically as follows: traverse the wave number kh∈[0,(p+1)π], take points at preset intervals, determine a unique transfer matrix according to the value of each kh, and solve the eigenvalue Λ of the corresponding transfer matrix H, where p is the polynomial order.

5. The method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 1, characterized in that: The step S4 is specifically as follows: According to the eigenvalue Λ corresponding to each mode i i The modulus lengths are sorted and output in positive order in the p+1 equally divided interval on [0,π], that is, the ith wavenumber interval outputs the ith eigenvalue of the current modulus length sorting. When the continuity requirements are met, continuous dispersion main mode curves and dissipation main mode curves are obtained, and the horizontal coordinates of both are standardized wavenumbers. The ordinate of the main mode curve of dispersion is the dispersion rate, which is calculated as follows: The ordinate of the dissipation main mode curve is the dissipation rate, which is calculated as follows: Among them, ω is the corresponding spatial frequency, h represents the actual grid size, ωh is used in combination to reflect the dispersion rate and dissipation rate of the semi-discrete format, k is the spatial frequency, kh represents the wave number, p is the polynomial order, R, I correspond to the real part and imaginary part of the complex number.

6. A method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 5, characterized in that: The step S5 is specifically as follows: The absolute deviation between the main mode of dispersion and the reference value is integrated at low wavenumbers to obtain the dispersion fit The reference value is given by the equation Given, Pe is the grid Peclet number, and i is the imaginary unit.

7. The method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 5, characterized in that: The step S6 is specifically as follows: if a continuous dissipation curve cannot be extracted in step S4, the wave-like solution Project it onto each eigenvector of the transfer matrix H, that is, u=Vα, and use the inverse of the eigenvector matrix V to obtain the complex weight α, where x is the node coordinate of the finite element solution.

8. A method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 7, characterized in that: The step S7 is specifically as follows: Each component Z of the wavelike solution after eigendecomposition according to the complex weights m =α m v m Perform the double integral, v m is the mth eigenvector, α m is the mth complex weight, and the complex weight of each mode in the two-norm energy is obtained Represents Z m is the mean of , p is the polynomial order, where the integration over the reference unit is done using Gaussian integration.

9. A method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 8, characterized in that: The step S8 is specifically as follows: For complex weight c m Reweight the real part of Accumulate each mode to obtain That is, the joint modal dissipation curve, where That is, the dissipation value corresponding to the mth eigenvalue in step S4.

10. The method for quantitatively analyzing the dispersion characteristics of numerical algorithms according to claim 5, characterized in that: The step S9 is specifically as follows: Through the 1% rule, That is, the maximum k value that continuously satisfies the above deviation formula is k 1％ , obtain the maximum fitting range of the dispersion main mode and the dissipation main mode / dissipation joint curve, where, for dispersion, For dissipation,

Images