Topological optimisation technology of gravity dam section structure based on ANSYS partial differential equation operation

Xin Guan; Peng Yao; Reem Alotaibi; Hamdy M. Mohamed

Dziennik i szczegóły wydania

Applied Mathematics and Nonlinear Sciences

Informacje o czasopiśmie

Format: Czasopismo
eISSN: 2444-8656
Pierwsze wydanie: 01 Jan 2016
Częstotliwość wydawania: 2 razy w roku
Języki: Angielski

This paper uses the finite element method to explain the specific nature of numerical instability such as network dependence in the topology optimisation of engineering structures from the perspective of partial differential equations. Gaussian function filtering method reduces the global impact of local extremum on topology optimisation. Finally, the method is introduced into the topology optimisation of concrete gravity dams in hydraulic engineering, and the topology optimisation program is developed in conjunction with ANSYS software language to achieve the topology optimisation of building structures in hydraulic engineering from a technical perspective.

Keywords

ANSYS software
partial differential equations
hydraulic engineering
gravity dam
building structure
topology optimisation

MSC 2010

12H99

Artykuł

Introduction

Structural topology optimisation is the development of structural shape optimisation and an aspect of layout optimisation. As the shape optimisation gradually matured, the new concept of structural topology optimisation began to develop. Now topology optimisation has become a new hot research topic in the field of international structural optimisation. In the shape optimisation process, the initial structure and the final structure have the same topology. For example, after the shape optimisation of a plate-like structure with three openings, the only change is the boundary shape of the openings. The number of openings has not increased or decreased but is still three. However, there may be a situation where the change in the number of openings may be more effective in reducing the weight of the plate than the change in the shape of the openings under the same design constraints. This is the original research purpose of topology optimisation. Although both finite element and boundary element can be used to automatically divide the mesh, there is no research on the data processing of a model with a topology change, so it is very difficult to automatically generate openings in the design area. To overcome this limitation, consider using a ‘fixed’ finite element model in which the less stressed elements are artificially designated with very soft materials to approximate the openings.

Structural topology optimisation discusses the interconnection methods of structural members, including the topological forms such as the presence or absence of voids and the number and location of holes in the structure, so that the structure can transfer external loads to the support under the constraints of equilibrium, stress and displacement. Seat, and at the same time make certain morphological indicators of the structure optimal. The main difficulty of topology optimisation is that the structural topology that meets certain functional requirements has many or even infinite forms, and this topology form of the target structure is difficult to quantitatively describe or parameterise, and because the area to be designed is unknown in advance, more added difficulty in solving problems. Strictly speaking, logical variables (presence or absence) or integer variables should be used to indicate whether there are rods connected between nodes, but this will complicate the optimisation model; so people are used to starting from the section optimisation model for topology. Optimisation, in simple terms, uses the theory of finite elements to convert topology optimisation into a form of size optimisation for calculation [1].

Based on the above research background, this paper starts from the mathematical theory of topology and in-depth studies the theory of structural topology optimisation and applies it to the practice of hydraulic engineering. The main work is as follows: First, the thesis explains the numerical instability phenomena such as grid dependence in structural topology optimisation from the perspective of the inverse problem of partial differential equations in the proposed structure according to the specific situation of the structural design of hydraulic buildings. In essence, combining the Gaussian function filtering method in digital signal processing and the multigrid method in fluid mechanics, using the convolution operation of the Gaussian function and the nature of the softenable kernel function, comprehensive use of the filtering method and multigrid method is proposed. The numerical instability is solved, and the filtering radius extension method is used to reduce the influence of local extreme values on the global optimisation of the structure [2]. The calculations show that the choice of the scheme is reasonable and effective, saving a lot of calculation time. At the same time, this paper chooses ANSYS software to realise the calculation of structural topology optimisation, and from the technical point of view, it realises the structural topology optimisation of hydraulic buildings.

Structural topology optimisation based on homogenisation method theory

Topology describes a special type of graphics, the so-called ‘topological properties’. However, although the topological property is a basic property of graphics, it also has strong geometric intuition; but it is difficult to accurately describe it in simple and popular language. Its exact definition is described in abstract language. Here we try to describe topology in terms of topological properties. Topological properties reflect the characteristics of the overall structure of the graphic. The graphic can be deformed (such as squeezed, stretched or twisted) at will; if it is not torn, no adhesion occurs, and the overall structure is not damaged. Nature remains the same. The above-mentioned deformations are called ‘topological transformations’ of the graphics, and then the topological properties are the properties that the geometric figures remain unchanged during the topological transformation. If the topology is described in terms of sets and mappings [3], topology is the study of the invariant or invariant properties of space under a topological transformation (or homeomorphism). The so-called homeomorphic space X and Y means that there is a two-way continuous (that is, reciprocal and continuous) correspondence between X and Y, which is figuratively that plasticine X can be kneaded into Y without allowing partitioning. Topological properties are geometric properties that are common to homeomorphic figures, so the two spaces X and Y of homeomorphic cannot be distinguished. Because the nature of topological research does not change when graphics are elastically deformed, topology is figuratively called ‘rubber geometry.’

The homogenisation method was proposed by Bendsoe and Kikuchi. It is the most widely used method for topology optimisation of continuum structures and belongs to the material description method. The basic idea is to introduce the microstructure (unit cell) into the material of the topological structure [4]. The form and size parameters of the microstructure determine the elastic properties and density of the micromaterial at this point. The optimisation process takes the size of the unit cell of the microstructure as topological design variables, the addition and deletion of microstructures are achieved by the growth and decline of the unit cell size, and composite materials composed of intermediate unit cells are generated to expand the design space and realise the unification and continuity of the structural topology optimisation model and the size optimisation model. The work on the homogenisation method mainly includes the research of the microstructure model theory and the practical application of the homogenisation model. The research on the theory of microstructure model mainly including proposed models of square guarding microstructure, two-level hierarchical microstructure, rectangular hollow microstructure, three-dimensional hierarchical microstructure, etc. pointed out that the orthogonal microstructure hypothesis will inevitably lead to incorrect results. Research on the application of the homogenisation model in the optimisation design of continuum topology. The research scope covers multi-load plane problems, 3D continuum problems, vibration problems, elasticity problems, buckling problems, 3D shell problems, thin shell structure problems, composite topology, optimisation issues and many other issues. The homogenisation theory is to divide the microstructure of many different holes in the design area (see Figure 1) to perform topology optimisation on the continuum.

2.1

Mathematical formula of the homogenisation process

In the process of homogenisation, it is difficult to keep all important physical quantities conserved before and after homogenisation. To prove this, consider the multigroup approximation of the neutron transport equation and assume that we can get the solution of the equation: (1) $\nabla \cdot J_{g} (r) + Σ_{tg} (r) Φ_{g} (r) = \sum_{g = 1}^{G} [1 / k_{eff} M_{gg}^{'} (r) + Σ_{gg}^{'} (r)] Φ_{g}^{'} (r)$ \nabla \cdot {{\rm{J}}_{\rm{g}}}(r) + {\Sigma _{tg}}(r){\Phi _g}(r) = \sum\limits_{g = 1}^G \left[ {1/{k_{eff}}M_{gg}^\prime (r) + \Sigma _{gg}^\prime (r)} \right]\Phi _g^\prime (r) In the formula J_g(r) = ∫ dΩΩ · Φ_g(r,Ω), Φ_g(r) = ∫ dΩΦ_g(r,Ω), $M_{{gg}^{'}} (r) = χ_{g} υ Σ_{fg}^{'} (r)$ {{\rm{M}}_{{\rm{g}}{{\rm{g}}^\prime }}}(r) = {\chi _g}\upsilon \Sigma _{fg}^\prime (r) , $Σ_{{gg}^{'}} (r) = 1 / 2 \int d μ_{0} Σ_{gg}^{'} (r, μ_{0})$ {\Sigma _{{\rm{g}}{{\rm{g}}^\prime }}}(r) = 1/2\int d{\mu _0}\Sigma _{gg}^\prime (r,\mu_0) , μ₀ = Ω · Ω′, k_eff is the eigenvalue of the reactor, Φ_g(r,Ω) is the neutron angular flux of the g group, G is the number of energy groups, and the symbol of the reaction cross section is in the conventional form.

After solving equation (1), the physical quantities that need to be conserved during the cell homogenisation process can be easily calculated. The three most important physical quantities are the average response rate of the grid element, the average interfacial flow and the eigenvalue (or proliferation coefficient). Like equation (1), the homogenisation equation can be obtained [5]: (2) $\nabla \cdot {\hat{J}}_{g} (r) + {\hat{Σ}}_{tg} (r) {\hat{Φ}}_{g} (r) = \sum_{g = 1}^{G} [1 / k_{eff} {\hat{M}}_{gg}^{'} (r) + {\hat{Σ}}_{gg}^{'} (r)] {\hat{Φ}}_{g}^{'} (r)$ \nabla \cdot {{\rm{\hat J}}_{\rm{g}}}(r) + {\hat \Sigma _{tg}}(r){\hat \Phi _g}(r) = \sum\limits_{g = 1}^G \left[ {1/{k_{eff}}\hat M_{gg}^\prime (r) + \hat \Sigma _{gg}^\prime (r)} \right]\hat \Phi _g^\prime (r) In the formula ${\hat{J}}_{g} (r) = \int d Ω Ω \cdot {\hat{Φ}}_{g} (r, Ω)$ {{\rm{\hat J}}_{\rm{g}}}(r) = \int d\Omega \Omega \cdot {\hat \Phi _g}(r,\Omega ) , ${\hat{Φ}}_{g} (r) = \int d Ω {\hat{Φ}}_{g} (r, Ω)$ {\hat \Phi _{\rm{g}}}(r) = \int d\Omega {\hat \Phi _g}(r,\Omega ) , ${\hat{M}}_{gg}^{'} (r) = χ_{g} υ {\hat{Σ}}_{fg}^{'} (r)$ {\rm{\hat M}}_{{\rm{gg}}}^\prime (r) = {\chi _g}\upsilon \hat \Sigma _{fg}^\prime (r) , ${\hat{Σ}}_{gg}^{'} (r) = 1 / 2 \int d μ_{0} {\hat{Σ}}_{gg}^{'} (r, μ_{0})$ \hat \Sigma _{{\rm{gg}}}^\prime (r) = 1/2\int d{\mu _0}\hat \Sigma _{gg}^\prime (r,{\mu _0}) , ${\hat{Σ}}_{gg}^{'}$ \hat \Sigma _{gg}^\prime , ${\hat{Σ}}_{fg}^{'}$ \hat \Sigma _{fg}^\prime , and ${\hat{Σ}}_{tg}$ {\hat \Sigma _{tg}} are the reaction cross sections that are independent of r after the cell homogenisation. In general, these sections are non-uniform throughout the area, that is, they are related to r.

The relationship between homogenisation parameters and non-uniform parameters is derived below according to the conservation of important physical quantities before and after homogenisation: (3) $\begin{array}{l} {\hat{Σ}}_{α g} \int_{V} {\hat{Φ}}_{g} (r) dr = \int_{V} Σ_{α g} Φ_{g} (r) dr \\ g = 1, 2, ..., G \\ α = t, g g^{'}, a, etc . \end{array}$ \matrix{{{{\hat \Sigma }_{\alpha g}}\int_V {{\hat \Phi }_g}(r)dr = \int_V {\Sigma _{\alpha g}}{\Phi _g}(r)dr} \hfill \cr {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;g = 1,2,...,G} \hfill \cr {\;\;\;\;\,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\alpha = t,g{g^\prime },a,etc.} \hfill \cr } (4) $\int_{S} {\hat{J}}_{g} (r) \cdot dS = \int_{S} J_{g} (r) \cdot dS$ \int_S {\hat J_g}(r) \cdot dS = \int_S {J_g}(r) \cdot dS In the formula, S and V are the surface area and volume of the homogenised region, respectively.

2.2

The mathematical paradox in homogenisation

Due to the strong coupling between the homogenisation cross section and the homogenisation flux in the transport equation, it is not easy to calculate the homogenisation parameters to make equations (3) and (4) true. On the other hand, during the homogenisation process, a total of G × (N + 1) homogenisation parameters (G is the total number of energy groups, N + 1 is the number of reaction types plus the diffusion coefficient) are required, and G × (N + 2 × D) physical quantities are to be conserved (G × N reaction rates and G × 2 × D homogenisations). The interface flow on each surface of the region, where D represents the dimension of the problem in the Cartesian coordinate system, and the degree of freedom of the homogenisation problem is insufficient. This reality shows that it is impossible to accurately maintain the conservation of important physical quantities during cell homogenisation. Equivalent homogenisation theory and generalised equivalent homogenisation theory bypass this mathematical problem by introducing discontinuous factors. For example, the SPH method adds G degrees of freedom (SPH factors), which can be applied to one-dimensional or multidimensional and symmetrical problems but is not applicable to general two-dimensional or three-dimensional problems. For all important physical quantities to be accurately conserved, more degrees of freedom need to be introduced. The symmetry of the module structure can effectively reduce the complexity in the design of the reactor and make it possible to achieve accurate and uniform theory. However, modern and future nuclear reactor designs will be more non-uniform, and as a result, the external environment of different cells will be more asymmetric.

Another feature of homogenisation is that it is sensitive to the external environment of the cell. It can be assumed that the core is the same for all cells, and there is no need to uniformise each cell at this time. However, in most realistic reactor problems, due to the different positions in the component, even the same cells will cause a difference in its external environment. People want to reduce the number of homogenisations of the same type of cells (basically, each type is performed only once), so the homogenisation parameters of a certain type of cells in a certain environment may not guarantee important physical quantities in the same cell in different environments.

Conservation

This issue will be addressed in the proposed homogenisation technique.

2.3

Semi-uniform iteration method

The semi-homogenisation method guarantees the conservation of all-important physical quantities by introducing enough degrees of freedom. An iterative method is used to solve the homogenisation parameters.

Definition

Let R_αg = §_V Σ_αg(r)Φ_g(r)dr. In the first iteration, the homogenisation parameters were obtained using the FVW method from the following formula: (5) ${\hat{Σ}}_{α g}^{1} = \frac{R_{α g}}{\int_{V} Φ_{g} (r) dr}$ \hat \Sigma _{\alpha g}^1 = {{{R_{\alpha g}}} \over {\int_V {\Phi _g}(r)dr}} (6) ${\hat{D}}_{g}^{1} = \frac{\int_{V} D_{g} Φ_{g} (r) dr}{\int_{V} Φ_{g} (r) dr}$ \hat D_g^1 = {{\int_V {D_g}{\Phi _g}(r)dr} \over {\int_V {\Phi _g}(r)dr}} Solve equation (2) to get ${\hat{Φ}}_{g}^{1} (r)$ \hat \Phi _g^1(r) . The interface flow is denoted as J_g,_k (g is the energy group label and k is the surface label), the uniformised interface flow is denoted as ${\hat{J}}_{g, k}$ {\hat J_{g,k}} and the semi-uniform diffusion coefficient is denoted as ${\hat{D}}_{g, j}$ {\hat D_{g,j}} . Among them, the diffusion coefficient is related to spatial coordinates. For example, for a two-dimensional problem, the cell can be divided into four parts, each of which corresponds to a diffusion coefficient, and the homogenisation section is for the entire cell.

For one-dimensional problems with two surfaces, the homogenisation parameters can be determined by: (7) ${\hat{Σ}}_{α g}^{2} = \frac{R_{α g}}{\int_{V} {\hat{Φ}}_{g}^{1} (r) dr}$ \hat \Sigma _{\alpha g}^2 = {{{R_{\alpha g}}} \over {\int_V \hat \Phi _g^1(r)dr}} (8) ${\hat{D}}_{g, 1}^{2} = {\hat{D}}_{g}^{1} (1 - α (1 - \frac{{\hat{J}}_{g, 1}^{1}}{J_{g, 1}}))$ \hat D_{g,1}^2 = \hat D_g^1(1 - \alpha (1 - {{\hat J_{g,1}^1} \over {{J_{g,1}}}})) (9) ${\hat{D}}_{g, 2}^{2} = {\hat{D}}_{g}^{1} (1 - α (1 - \frac{{\hat{J}}_{g, 2}^{1}}{J_{g, 2}}))$ \hat D_{g,2}^2 = \hat D_g^1(1 - \alpha (1 - {{\hat J_{g,2}^1} \over {{J_{g,2}}}})) where α is an accelerated convergence coefficient. In this way, the second iteration calculation can be completed. In the i-th step iteration calculation, the homogenisation parameters are calculated as follows: (10) ${\hat{Σ}}_{α g}^{i} = \frac{R_{α g}}{\int_{V} {\hat{Φ}}_{g}^{i - 1} (r) dr}$ \hat \Sigma _{\alpha g}^i = {{{R_{\alpha g}}} \over {\int_V \hat \Phi _g^{i - 1}(r)dr}} (11) ${\hat{D}}_{g, 1}^{i} = {\hat{D}}_{g, 1}^{i - 1} (1 - α (1 - \frac{{\hat{J}}_{g, 1}^{i - 1}}{J_{g, 1}}))$ \hat D_{g,1}^i = \hat D_{g,1}^{i - 1}(1 - \alpha (1 - {{\hat J_{g,1}^{i - 1}} \over {{J_{g,1}}}})) (12) ${\hat{D}}_{g, 2}^{i} = {\hat{D}}_{g, 2}^{i - 1} (1 - α (1 - \frac{{\hat{J}}_{g, 2}^{i - 1}}{J_{g, 2}}))$ \hat D_{g,2}^i = \hat D_{g,2}^{i - 1}(1 - \alpha (1 - {{\hat J_{g,2}^{i - 1}} \over {{J_{g,2}}}})) It can be seen from the above that the homogenisation parameters obtained after the convergence of the iterative calculation can ensure the accurate conservation of all reaction rates and interfacial flows. The next section compares this method with other traditional homogenisation methods.

Numerical solution to structural instability optimisation based on ANSYS software and partial differential equations

3.1

ANSYS software

ANSYS program is a powerful and flexible large-scale general-purpose finite element commercial analysis software for design analysis and optimisation, integration of structure, heat, fluid, electromagnetic and acoustics. It is widely used in various industrial and scientific researches. The ANSYS program has basically realised seamless integration with CAD software: it has a strong and reliable meshing capability; the program has explicit integral calculations in the structure and multi-field analysis so that it can solve from linear problems to nonlinear problems and has the evolution from solving the structure field to solving the wobble field; at the same time, the ANSYS program has an open environment, and users can add their own tasks to the software to complete special tasks as needed [6].

The ANSYS program itself provides a direct use program for structural topology optimisation. At the same time, because of its powerful compatibility, if you want to use ANSYS to design a program for topology optimisation, you only need to combine it with the FORTRAN language for hybrid analysis. You can get the desired result. Because of the visual effects of ANSYS, the results of structural topology optimisation can be clearly seen, so the ANSYS program can be said to be a good implementation tool for structural topology optimisation.

3.2

Grid dependency

3.2.1

Cause

The topology optimisation problem of a structure is equivalent to an inverse problem of a partial differential equation. The finite element method is a weak solution to the original problem on a finite-dimensional space S_h: (13) $a (μ_{h}, φ) = (f, φ) \forall φ \in S_{h}$ a\left( {{\mu _h},\varphi } \right) = \left( {f,\varphi } \right)\quad \forall \varphi \in {S_h} As the size of the finite element mesh becomes more and more fine, its numerical solution gradually approaches the physical understanding of the original problem, that is, the solution of the finite element method should be convergent. When using a thin initial mesh for structural topology optimisation, the original intention is to obtain a structural topology optimisation form that is satisfactory to the designer to avoid using the shape optimisation method to perform tedious post-processing on the obtained topology structure form. However, a strange phenomenon appears in the actual topology optimisation design, that is, the grid dependence of the structural topology optimisation solution. The so-called grid dependency means that for the same initial design area, when using different initial grids for optimisation, the optimal topological form is different, that is, the topology optimisation results depend on the division of the initial grid (see Figure 2). In this way, the thinner initial mesh corresponds to the more complex topological structure, and some of the details in the optimisation result are already smaller [7], which is very fragile relative to the overall structure size, and it does not meet the actual engineering requirements.

In fact, the grid dependency problem of the topology optimisation solution can be given a clear explanation by the characteristics of the integral equation solution. As can be seen from the previous section, the solutions to structural topology optimisation problems are ill-posed, but there are certain connections between the solutions. However, to find the numerical solution of the first type of Fredholm integral equation, the integral equation is first discretised to obtain: (14) $\sum_{j = 1}^{n} K (x, y_{j}) ω_{j} f (y_{j}) = μ (x)$ \sum\limits_{j = 1}^n K\left( {x,{y_j}} \right){\omega _j}f\left( {{y_j}} \right) = \mu \left( x \right) In the formula, y_j (j = 1,2,3,...,n) represents the integration node, and ω_j represents the weight factor of the numerical integration. The value of the above formula on x = x_i (i = 1,2,3,...,n) is obtained as: (15) $\sum_{j = 1}^{n} K (x_{i}, y_{j}) ω_{j} f (y_{j}) = μ (x_{i})$ \sum\limits_{j = 1}^n K\left( {{x_i},{y_j}} \right){\omega _j}f\left( {{y_j}} \right) = \mu \left( {{x_i}} \right) If the discrete node (x_i,y_j) for x,y has been selected, the above equation is a system of linear algebraic equations. Given the matrix K = {k (x_i,y_j)}, the right term U = u{(x_i)} and the unknown quantity {ω_j f (y_j), obviously when two different numerical integration schemes are used, there are: (16) $ω_{j}^{(1)} f^{(1)} (y_{j}) = ω_{j}^{(2)} f^{(2)} (y_{j}) (j = 1, 2, 3, ..., n)$ \omega _j^{\left( 1 \right)}{f^{\left( 1 \right)}}\left( {{y_j}} \right) = \omega _j^{\left( 2 \right)}{f^{\left( 2 \right)}}\left( {{y_j}} \right)\quad \left( {j = 1,2,3,...,n} \right) Because ${ω_{j}^{(1)}} \neq {ω_{j}^{(2)}}$ \left\{ {\omega _j^{\left( 1 \right)}} \right\} \ne \left\{ {\omega _j^{\left( 2 \right)}} \right\} , {f⁽¹⁾ (y_j)} ≠ {f⁽²⁾ (y_j)}. So when the integral equation is discretised in different ways, the resulting numerical solution will not be stable, and often the finer the mesh is divided, resulting in matrix condition number of the linear algebraic equations. The bigger the solution, the more unstable the solution is and the more complicated the topology is.

3.2.2

Finite element analysis method

Because the finite element method is an approximate method and the word length of a computer is limited, the calculated solution is only an approximate solution of a numerical model, that is, there is a certain error between the numerical solution and the physical understanding of the original problem. The so-called high-pass and low-pass terms of the error are relative to the smoothness of the numerical solution of the original problem. When the scale parameter of the function is constant relative to the selected grid size, the high-pass error components that can be eliminated for different initial grid functions are different. The multigrid method is an iterative method for solving large-scale equations and has been widely used in various disciplines and engineering and technical problems. The multigrid method can eliminate the low-pass part of the numerical solution error and accelerate the convergence of the solution to effectively solve the grid dependency problem in topology optimisation. On the premise of considering only static loads, according to the principle of virtual work in elasticity, the finite element equilibrium equation for plane stress problems is: (17) $KU = F$ KU = F Its node displacement expression is: (18) $U = K^{- 1} F$ U = {K^{ - 1}}F In the formula, F represents the node load vector; U represents the node displacement vector; K represents the total stiffness matrix of the system, which is assembled from all element stiffness matrices and is a symmetric matrix, K = K^T.

Isoperimetric elements use two sets of coordinate systems, one is the actual coordinate system where the irregular quadrilateral elements are located, called the xy coordinate system, and the other is the local coordinate system where the square elements (parent elements) are located, called the ξ η coordinate system. To consider the actual unit as the image of the parent unit, it is necessary to establish the transformation relationship between the two sets of coordinate systems, and this is a problem describing the coordinate vector field. In unit analysis, many of them are related to the partial derivative of the row function with respect to the overall coordinates. According to the derivation rules of the composite function, there are (19) ${\begin{matrix} \frac{\partial N_{i}}{\partial x} = \frac{\partial N_{i}}{\partial ξ} \frac{\partial ξ}{\partial x} + \frac{\partial N_{i}}{\partial η} \frac{\partial η}{\partial x} \\ \frac{\partial N_{i}}{\partial y} = \frac{\partial N_{i}}{\partial ξ} \frac{\partial ξ}{\partial y} + \frac{\partial N_{i}}{\partial η} \frac{\partial η}{\partial y} \end{matrix}$ \left\{ {\matrix{{{{\partial {N_i}} \over {\partial x}} = {{\partial {N_i}} \over {\partial \xi }}{{\partial \xi } \over {\partial x}} + {{\partial {N_i}} \over {\partial \eta }}{{\partial \eta } \over {\partial x}}} \cr {{{\partial {N_i}} \over {\partial y}} = {{\partial {N_i}} \over {\partial \xi }}{{\partial \xi } \over {\partial y}} + {{\partial {N_i}} \over {\partial \eta }}{{\partial \eta } \over {\partial y}}} \cr } } \right. x_i,y_i represent unit node coordinates; N_i represents row functions and local coordinates represent interpolation functions. To facilitate numerical integration, transform it into a local coordinate system. According to the area where the parent element is located, the expression of the stiffness matrix of the 4-node isoperimetric element is: (20) $K^{ε} = t \int_{- 1}^{1} \int_{- 1}^{1} {[B]}^{J} [D] [B] Jd ξ d η$ {K^\varepsilon } = t\int_{ - 1}^1 \int_{ - 1}^1 {\left[ B \right]^J}\left[ D \right]\left[ B \right]Jd\xi d\eta Among them, [B] represents the geometric matrix, [J] represents the Jacobian matrix, [D] represents the elastic matrix of the stress–strain relationship and t represents the element thickness. The Gaussian numerical integration method with higher efficiency is used for the above equation integration. The two-dimensional Gaussian integral formula is: (21) $\int_{- 1}^{1} \int_{- 1}^{1} f (s, t) dsdt = \sum_{i = 1}^{n} \sum_{j = 1}^{n} H_{i} H_{j} f (s_{i}, t_{j})$ \int_{ - 1}^1 \int_{ - 1}^1 f\left( {s,t} \right)dsdt = \sum\limits_{i = 1}^n \sum\limits_{j = 1}^n {H_i}{H_j}f\left( {{s_i},{t_j}} \right) In the formula, f (s_i,t_j) represents the function value of the integrand at the integration point s_i,t_j, H_i, H_j represents the two-dimensional weighting coefficient of the integration point and n represents the number of Gaussian integration points.

3.2.3

Grid-independent filtering method

The grid-independent filtering method is to modify the sensitivity information value of a unit by weighting the sensitivity information of neighbouring units in the surrounding neighbourhood of a specific unit to obtain the redistribution of the filtered unit sensitivity information as the next iteration.

Calculated initial value. Compared with other methods that try to eliminate the high-pass part of the grid-dependent numerical solution error, the grid-independent filtering method has little effect on the calculation time, and the algorithm is easy to implement and has been widely used in structural topology optimisation [8].

As mentioned earlier, the filtering method is to use different filter functions as the convolution kernel to convolve with the original digital signal. Here are some explanations on the selection of the convolution kernel (filter function). There are many classic filter functions in digital signals. In general, the selection of the filter function should follow the two principles: (1) The main lobe width of the frequency response of the filter function should be small, and the energy it contains should be as large as possible in the total energy. (2) The energy contained in the side lobe of the frequency response of the filter function varies with the bandwidth, which tends to decrease rapidly.

The above two principles are selected to consider the phenomenon of spectral leakage when the original function and the convolution kernel are convolved. Since the purpose of suppressing spectrum leakage can be achieved by selecting an appropriate filter function, how to choose an optimal filter function becomes the key. According to the Hegenss principle, the product of the time width T and the bandwidth ω of any signal satisfies $T ω \geq \frac{1}{4} π$ T\omega \ge {1 \over 4}\pi . The filter function obtained under the optimisation criterion $T ω = \frac{1}{4} π$ T\omega = {1 \over 4}\pi of the minimum time width and bandwidth product is a Gaussian function: (22) $G (x) = {\begin{array}{l} \frac{1}{\sqrt{2 π τ^{2}}} exp (- \frac{x 2}{2 τ 2}) & | x | τ \\ 0 & | x | \geq τ \end{array}$ G\left( x \right) = \left( {\matrix{{{1 \over {\sqrt {2\pi {\tau ^2}} }}\exp \left( { - {{x2} \over {2\tau 2}}} \right)} \hfill & {\left| x \right|\tau } \hfill \cr 0 \hfill & {\left| x \right| \ge \tau } \hfill \cr } } \right. where τ represents the scale parameter of the Gaussian function.

Application of topology optimisation theory in optimal design of gravity dam

4.1

Basic situation of load and section

According to the formula of Granting Reservoir, the wave height is 2HL = 1.8 m, the wavelength is 2LL = 16.64 m and the average wave centreline is 0.85 m above the surface of the still water. The water depth in front of the dam is 95 m, which is greater than LL and is a deep-water wave. The calculated wave pressure per metre of single width is 109.9 kN. The load considers the basic load combination situation and considers upstream water pressure, downstream water pressure, lift pressure and gravity and wave pressure and does not consider sediment pressure and earth pressure. The topological optimisation of the structure is performed for the downstream water levels of 0 and 10 m. According to the basic requirements of the section design of the gravity dam, the height of the concrete gravity dam is initially set at 100 m. Because of the drainage pipe, the lifting pressure is effectively reduced, and so the width of the dam bottom is set to 75 m and the width of the dam top is 8 m [9, 10].

4.2

Structural topology optimisation design of a non-null section of concrete gravity dam when the downstream water level is 0

Set the basic design area as a 100 m × 75 m rectangular plane (see Figure 3), where W is the dead weight, P₁ is the upstream water pressure, U is the lift pressure, and P_L is the wave pressure.

Schematic diagrams of the basic structural optimisation of the non-overflow section of a concrete gravity dam when the downstream water level is 0

The grid size is divided into 1 m × 1 m, the concrete bulk density is 24 kN/m³ and the water bulk density is 10 kN/m³. The downstream water level is 0. The width of the dam top is 8 m. In order to make the dam top have enough width, topology optimisation is forbidden at 0–8 m upstream to obtain a practical profile. The number of deleted element proportions is 55% of the optimised area. After the optimised results are obtained, the results are checked for sliding resistance and stress. To illustrate the problem, we do not take any measures to eliminate numerical instability in topology optimisation such as checkerboard format and grid dependency [11]. The parts of the dam crest width range are not subjected to topology optimisation. The results obtained are shown in Figure 4(a) and (b). Figure 4(a) shows the pseudo-density of the entire structure after topology optimisation, and the dashed box represents the original design area. Figure 4(b) shows the structure of the topology optimised in the form of a unit.

Topology optimisation unit result display

Summary

The thesis is closely related to the actual situation of hydraulic engineering. Based on the fractional partial differential equations and the theory of finite element analysis, the topological optimisation design of the non-overflow section of the concrete gravity dam, which is one of the most common forms of hydraulic engineering building structures, is proposed, a new type of dam section. The influence of the downstream water level pressure in the form of topology optimisation was compared. The results of the existing examples are used to compare the results of topology optimisation with commonly used optimisation results. The results prove the superiority of the topology optimisation idea.

Ilustracje i tabele

Referencje

[1] Efremova L S. Small C1-smooth perturbations of skew products and the partial integrability property[J]. Applied Mathematics and Nonlinear Sciences, 2020, 5(2):317–328. EfremovaL S Small C1-smooth perturbations of skew products and the partial integrability property [J]. Applied Mathematics and Nonlinear Sciences 2020 5 2 317 328 10.2478/amns.2020.2.00057 Search in Google Scholar

[2] Itil H G. Investigation of A Fuzzy Problem by the Fuzzy Laplace Transform[J]. Applied Mathematics and Nonlinear Sciences, 2019, 4(2):407–416. ItilH G Investigation of A Fuzzy Problem by the Fuzzy Laplace Transform [J]. Applied Mathematics and Nonlinear Sciences 2019 4 2 407 416 10.2478/AMNS.2019.2.00039 Search in Google Scholar

[3] Banerjee A, Paul D K, Acharyya A. Optimization and safety evaluation of concrete gravity dam section[J]. KSCE Journal of Civil Engineering, 2015, 19(6):1612–1619. BanerjeeA PaulD K AcharyyaA Optimization and safety evaluation of concrete gravity dam section [J]. KSCE Journal of Civil Engineering 2015 19 6 1612 1619 10.1007/s12205-015-0139-0 Search in Google Scholar

[4] Meena C, Praveen T V, Parvathi I S. Coupled Seepage–Stress Analysis of Selected Roller-Compacted Concrete Gravity Dam with Upstream Impervious Layer Considering the Seepage Behavior of Interface Thickness[J]. Indian Geotechnical Journal, 2020, 50(4):634–645. MeenaC PraveenT V ParvathiI S Coupled Seepage–Stress Analysis of Selected Roller-Compacted Concrete Gravity Dam with Upstream Impervious Layer Considering the Seepage Behavior of Interface Thickness [J]. Indian Geotechnical Journal 2020 50 4 634 645 10.1007/s40098-019-00376-4 Search in Google Scholar

[5] Chen D H, Yang Z H, Wang M, et al. Seismic performance and failure modes of the Jin’anqiao concrete gravity dam based on incremental dynamic analysis[J]. Engineering Failure Analysis, 2019, 100:227–244. ChenD H YangZ H WangM Seismic performance and failure modes of the Jin’anqiao concrete gravity dam based on incremental dynamic analysis [J]. Engineering Failure Analysis 2019 100 227 244 10.1016/j.engfailanal.2019.02.018 Search in Google Scholar

[6] Su H, Li J, Guo Z, et al. Nonprobabilistic Reliability Evaluation for In-Service Gravity Dam Undergoing Structural Reinforcement[J]. IEEE Transactions on Reliability, 2018, 67(3):970–986. SuH LiJ GuoZ Nonprobabilistic Reliability Evaluation for In-Service Gravity Dam Undergoing Structural Reinforcement [J]. IEEE Transactions on Reliability 2018 67 3 970 986 10.1109/TR.2018.2827919 Search in Google Scholar

[7] Petaccia G, Lai C G, Milazzo C, et al. The collapse of the Sella Zerbino gravity dam[J]. Engineering Geology, 2016, 211:39–49. PetacciaG LaiC G MilazzoC The collapse of the Sella Zerbino gravity dam [J]. Engineering Geology 2016 211 39 49 10.1016/j.enggeo.2016.06.024499770827594709 Search in Google Scholar

[8] Wang Y Y, Liu J Y, Ning B K, et al. Seismic Response Analysis of Reservoir Water-Gravity Dam-Foundation System[J]. Key Engineering Materials, 2017, 737:494–499. WangY Y LiuJ Y NingB K Seismic Response Analysis of Reservoir Water-Gravity Dam-Foundation System [J]. Key Engineering Materials 2017 737 494 499 10.4028/www.scientific.net/KEM.737.494 Search in Google Scholar

[9] Pausz S, Nowotny H, Jung G. Rock mass classification and geotechnical model for the foundation of a RCC gravity dam[J]. Geomechanik Und Tunnelbau, 2015, 8(5):436–440. PauszS NowotnyH JungG Rock mass classification and geotechnical model for the foundation of a RCC gravity dam [J]. Geomechanik Und Tunnelbau 2015 8 5 436 440 10.1002/geot.201500020 Search in Google Scholar

[10] Hu J, Wu S. Statistical modeling for deformation analysis of concrete arch dams with influential horizontal cracks[J]. Structural Health Monitoring, 2019, 18(2):546–562. HuJ WuS Statistical modeling for deformation analysis of concrete arch dams with influential horizontal cracks [J]. Structural Health Monitoring 2019 18 2 546 562 10.1177/1475921718760309 Search in Google Scholar

[11] Wang G, Zhang S, Lu W, et al. Damage effects of concrete gravity dams subjected to underwater explosion[J]. Journal of Hydraulic Engineering, 2015, 46(6):723–731. WangG ZhangS LuW Damage effects of concrete gravity dams subjected to underwater explosion [J]. Journal of Hydraulic Engineering 2015 46 6 723 731 Search in Google Scholar