# A Form Finding Method for Arch Bridges Using Parametric Level Set Method

Authors: Yadong Shen; Jianhu Feng; Xiaohan Cheng; Xuntao Wang; Changhao Zhang
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Topology optimization has developed rapidly in the past three decades; as a creative and efficient optimization technique, it has been applied in engineering fields of aerospace and mechanical. However, there are a few attempts in bridge form design. In this paper, the parametric level set method is utilized to solve the form finding of arch bridges. The optimization model for minimizing the structural compliance under the volume constraint is built. Three numerical examples of form finding of arch bridges are studied. Results show that the optimal structures which have well-distributed stress and smooth force transmission are almost identical with the actual forms of arch bridges. The optimal forms can be treated as alternatives in the preliminary design stage, and topology optimization has a bright prospect in form finding of arch bridges.


---

## Body

## 1. Introduction

In conventional process of engineering design, there is blindness that the initial form of structure is achieved according to the knowledge and experience of engineers. Topology optimization refers to achieving optimal distribution of a material in a given design domain and also obtaining the optimal mechanical performance under specified constraints simultaneously. It can provide a novel and reasonable design in the preliminary design stage. Based on these advantages, topology optimization is widely applied in different fields of engineering, for instance, aerospace, mechanical, electromagnetic, and electrochemical [1–4].Topology optimization has been developed rapidly in the past three decades. Numerous topology optimization methods including homogenization method (HM) [5, 6], density method (DM) [7–9], evolutionary structural optimization (ESO) [10, 11], and level set method (LSM) [12, 13] have been proposed. Compared with the previous three methods, LSM has a smooth boundary which contains a lot of geometry information, and it is easier to convert the optimal configuration into the CAD model for manufacturing. In addition, LSM handles complex boundary change easily. LSM refers to using of the zero isosurfaces of the level set functions to represent the structural boundaries. In conventional LSM, the evolution of structural boundaries is tracked by solving a Hamilton–Jacobi partial differential equation (H-J PDE). An upwind scheme was employed by Mei et al. [14] to solve the H-J PDE using the discrete solution. In this method, the accuracy and rate of convergence are affected as they are often too steep or too flat near the boundaries during the evolution of structural boundaries. Periodically, reinitializing the signed distance function is needed to alleviate numerical instability. And the optimal result generally depends on the initial topology because the new holes cannot be created in the design domain in the most conventional LSM.Parametric level set method (PLSM) solves these problems effectively without directly solving the H-J PDE. In this method, the implicit level set function is interpolated by radial basis functions (RBFs), the H-J PDE is converted into a series of ordinary differential equations which are easier to handle, and the propagation of boundaries is driven by updating interpolated coefficients using mathematical programming methods. The PLSM not only inherits the favorable features of the conventional LSM but also avoids the difficulties in solving the HJ-PDE. Wang et al. [15] firstly introduced globally supported radial basis functions (GSRBFs) into topology optimization of minimizing the structural compliance under the volume constraints. Luo et al. [16] applied the compactly supported radial basis functions (CSRBFs) [17, 18] into topology optimization of compliant mechanisms and established the optimization model with a nonconvex objective function and specified constraints. The evolution of the boundary was the process of updating the interpolated coefficients using the method of moving asymptotes (MMA) [19]. Luo et al. [20, 21] also proposed an effective method for structural stiffness design which combines the CSRBFs with well-established optimization algorithms, such as MMA and OC.Up to date, topology optimization has been widely used in fields of aerospace, mechanical, etc. However, it has seen little use in form finding of bridge. Chen et al. [22] obtained the optimal form of arch bridge using BESO and discussed the optimization results with different rise-to-span ratios. Felicetti et al. [23] applied BESO to bridge design in consideration of different requirements, for instance, support types, constraint types, and selections of span. However, this method has the numerical problem of checkerboard. Furthermore, the process of converting the optimized result to the CAD model needs a lot of postprocessing efforts.In this paper, PLSM which has a lot of geometry information is applied into form finding of arch bridges. The optimization model of minimizing the structural compliance under the volume constraint is proposed. And three examples of arch bridges are studied. The paper is composed as follows. In Section2, level set representation and parameterization using CSRBF are presented. After that, the topology optimization model is formulated. The sensitivity analysis and the flowchart of optimization are given in Section 3. Section 4 discusses the results of three examples in form finding of arch bridges. Finally, some conclusions are presented in Section 5.


## 2. Parametric Level Set Method

### 2.1. Level Set Representation

In the level set method, the structural boundaries are represented by the zero-level set of a higher dimensional level set function, which is usually defined as follows:(1)Φx>0,∀x∈Ω\Γ,Φx=0,∀x∈Γ,Φx<0,∀x∈D\Ω∪Γ,where Φx is the level set function, D is the design domain, Ω is the positive value of the level set function, and Γ denotes the zero-level set.A pseudotimet is introduced to manage the evolution of structural boundaries; differentiating both the sides of Φxt,t=0 with respect to t, we can obtain that(2)∂Φx,t∂t−vn∇Φx,t=0,where vn is the normal velocity of the level set function.


### 2.2. Parameterization Using CSRBF

The PLSM with compactly supported radial basis function (CSRBF) [20, 21] which has sufficiently smooth and accurate to approximate the level set function is applied into topology optimization of arch bridge. The level set function can be expressed as(3)Φx,t=∑i=1Nwiαit,where N is the total number of computational points, wi denotes the CSRBF at the ith point, and αi is the interpolated coefficient at the ith point.In this work, CSRBF with C2 continuity [17, 18] which has strictly positive definiteness and the sparseness of collection matrices is employed to approximate the level set function, and it can be expressed as(4)wix=max0,1−rix44rix+1,i=1,2,…,N,where r is the radius of support, and it is usually defined in a 2D Euclidean space as(5)rix=x−xiR,i=1,2,…,N,where xi is the coordinates of the ith computational point and R is the radius of support domain of the computational point x, which is defined as 2.0−4.0.By substituting Equation (3) into Equation (2), the HJ-PDE in Equation (2) can be rewritten as follows:(6)∑i=1Nwixα˙it−vn∑i=1N∇wixαit=0.The normal velocity of the level set function can be calculated by(7)vn=∑i=1Nwiα˙it∑i=1N∇wiαit,where a˙it=dai/dt.Equation (7) shows that vn is obtained by the variation of interpolated coefficients, and the shape change problem is transformed to a parameterization problem.


## 2.1. Level Set Representation

In the level set method, the structural boundaries are represented by the zero-level set of a higher dimensional level set function, which is usually defined as follows:(1)Φx>0,∀x∈Ω\Γ,Φx=0,∀x∈Γ,Φx<0,∀x∈D\Ω∪Γ,where Φx is the level set function, D is the design domain, Ω is the positive value of the level set function, and Γ denotes the zero-level set.A pseudotimet is introduced to manage the evolution of structural boundaries; differentiating both the sides of Φxt,t=0 with respect to t, we can obtain that(2)∂Φx,t∂t−vn∇Φx,t=0,where vn is the normal velocity of the level set function.


## 2.2. Parameterization Using CSRBF

The PLSM with compactly supported radial basis function (CSRBF) [20, 21] which has sufficiently smooth and accurate to approximate the level set function is applied into topology optimization of arch bridge. The level set function can be expressed as(3)Φx,t=∑i=1Nwiαit,where N is the total number of computational points, wi denotes the CSRBF at the ith point, and αi is the interpolated coefficient at the ith point.In this work, CSRBF with C2 continuity [17, 18] which has strictly positive definiteness and the sparseness of collection matrices is employed to approximate the level set function, and it can be expressed as(4)wix=max0,1−rix44rix+1,i=1,2,…,N,where r is the radius of support, and it is usually defined in a 2D Euclidean space as(5)rix=x−xiR,i=1,2,…,N,where xi is the coordinates of the ith computational point and R is the radius of support domain of the computational point x, which is defined as 2.0−4.0.By substituting Equation (3) into Equation (2), the HJ-PDE in Equation (2) can be rewritten as follows:(6)∑i=1Nwixα˙it−vn∑i=1N∇wixαit=0.The normal velocity of the level set function can be calculated by(7)vn=∑i=1Nwiα˙it∑i=1N∇wiαit,where a˙it=dai/dt.Equation (7) shows that vn is obtained by the variation of interpolated coefficients, and the shape change problem is transformed to a parameterization problem.


## 3. Topology Optimization for Form Finding of Arch Bridges

### 3.1. Optimization Model Using PLSM

In this optimization model, the interpolated coefficients of CSRBFs are adopted as design variables, optimization model of minimizing the structural compliance under the volume constraint is proposed, and it can be mathematically expressed as follows:(8)findαi,i=1,2,…,N,min.J=∫DεuDΦεudΩ,s.t.au,v,Φ=lv,Φ,GΦ≤fV0,α¯≤αi≤α¯,where J denotes the objective function, u is the displacement field, and v is the virtual temperature field. εu is the strain, V0 is the volume of design domain, GΦ=∫DHΦdΩ is the volume of structure, f denotes the volume constraint, α¯ and α¯ are the lower and the upper bounds of the design variables, and DΦ is the elastic stiffness which can be calculated as follows:(9)DΦ=HΦE,where E is the material elastic stiffness and HΦ is the Heaviside function. In order to avoid regenerating the element mesh, we employ a smoothed approximation of HΦ, which is defined as(10)HΦ=χ,Φ≤−Δ,31−χ4ΦΔ−13ΦΔ3+1+χ2,−Δ≤Φ≤Δ,1,Φ≥Δ,(11)δΦ=31−χ4Δ1−ΦΔ2,Φ≤Δ,γ,Φ≥Δ,where χ and γ are small positive numbers to avoid singularity of stiffness matrix and Δ is the width of the numerical approximation.a
u
,
v
,
Φ
=
l
v
,
Φ is the weak form of the state equation, and the energy bilinear form au,v,Φ and the load linear form lv,Φ can be expressed as follows:(12)au,v,Φ=∫DεuDΦεvdΩ,lv,Φ=∫DpvHΦdΩ+∫DcvδΦ∇ΦdΩ,where p is the body force and c is the boundary traction.


### 3.2. Sensitivity Analysis

The optimization problem in Equation (8) can be solved by sensitivity-based optimization algorithms. Thus, the derivatives of the objective function and volume constraints with respect to the design variables are required. The Lagrangian method is employed to transform the constrained optimization problem into an unconstrained problem, which can be expressed as(13)L=Ju,Φ+au,v,Φ−lv,Φ+λGΦ−fV0,where λ is the Lagrangian multiplier.The material derivative is used to obtain the shape derivative ofL. For more detailed theoretical proofs, the readers can refer to [20, 21]. The shape derivative of L can be obtained by(14)∂Lu,Φ∂t=∑i=1N∫Dβu,v,Φwi∂αi∂tdΩ+λ∑i=1N∫DδΦ∂αi∂tdΩ,and(15)βu,v,Φ=fu,u−εuCεv+pv+∇cv•n+κcv,where κ=divn=∇•−∇Φ/∇Φ is the mean curvature.The derivative of LagrangianL with respect to t can also be obtained by the chain rule as follows:(16)∂Lu,Φ∂t=∑i=1N∂Ju,Φ∂αi+λ∑i=1N∂GΦ∂αi∂αi∂t.By comparing the corresponding terms in Equations (14) and (16), the derivatives of the objective function and the volume constraint with respect to the interpolated coefficients can be obtained by(17)∂Ju,Φ∂αi=∫Dβu,v,ΦwidΩ,(18)∂GΦ∂αi=∫DδΦdΩ.When Equations (17) and (18) are obtained, the optimization problem in Equation (8) can be solved by sensitivity-based optimization algorithms.


### 3.3. Optimization Procedure

In optimization problem, with the sensitivities of objective function and volume constraint with respect to interpolated coefficients obtained, the sensitivity-based optimization algorithms such as optimality criteria and mathematical programming can be applied straightforwardly. Then, the shape and topology optimization problem converts into the parameter optimization problem. In this work, the method of moving asymptotes (MMA) which is a widely used mathematical programming algorithm is employed to solve the optimization problem. The process of optimization is shown in Figure1, and the corresponding introduction is presented as follows:Step 1: initialize the design domain of arch bridge. Define the load, boundary conditions, and set the volume constraint.Step 2: FE analysis and calculate the real and virtual displacement fields.Step 3: calculate the derivatives of objective function and the volume constraint with respect to the interpolated coefficients in Equations (17) and (18).Step 4: calculate the values of objective function and the volume constraint. Update the interpolated coefficients with MMA optimization algorithm, and then, the update level set function.Step 5: check whether the convergence condition is shown in Equation (19). If the terminal condition is met, stop, and the optimal structure is found. Otherwise, go back to Step 2. The terminal condition is defined as(19)Js−Js−1Js≤ε,where Js denotes the value of the objective function at sth iteration, and ε is a small positive value.Figure 1
Flowchart of form finding of arch bridge.


## 3.1. Optimization Model Using PLSM

In this optimization model, the interpolated coefficients of CSRBFs are adopted as design variables, optimization model of minimizing the structural compliance under the volume constraint is proposed, and it can be mathematically expressed as follows:(8)findαi,i=1,2,…,N,min.J=∫DεuDΦεudΩ,s.t.au,v,Φ=lv,Φ,GΦ≤fV0,α¯≤αi≤α¯,where J denotes the objective function, u is the displacement field, and v is the virtual temperature field. εu is the strain, V0 is the volume of design domain, GΦ=∫DHΦdΩ is the volume of structure, f denotes the volume constraint, α¯ and α¯ are the lower and the upper bounds of the design variables, and DΦ is the elastic stiffness which can be calculated as follows:(9)DΦ=HΦE,where E is the material elastic stiffness and HΦ is the Heaviside function. In order to avoid regenerating the element mesh, we employ a smoothed approximation of HΦ, which is defined as(10)HΦ=χ,Φ≤−Δ,31−χ4ΦΔ−13ΦΔ3+1+χ2,−Δ≤Φ≤Δ,1,Φ≥Δ,(11)δΦ=31−χ4Δ1−ΦΔ2,Φ≤Δ,γ,Φ≥Δ,where χ and γ are small positive numbers to avoid singularity of stiffness matrix and Δ is the width of the numerical approximation.a
u
,
v
,
Φ
=
l
v
,
Φ is the weak form of the state equation, and the energy bilinear form au,v,Φ and the load linear form lv,Φ can be expressed as follows:(12)au,v,Φ=∫DεuDΦεvdΩ,lv,Φ=∫DpvHΦdΩ+∫DcvδΦ∇ΦdΩ,where p is the body force and c is the boundary traction.


## 3.2. Sensitivity Analysis

The optimization problem in Equation (8) can be solved by sensitivity-based optimization algorithms. Thus, the derivatives of the objective function and volume constraints with respect to the design variables are required. The Lagrangian method is employed to transform the constrained optimization problem into an unconstrained problem, which can be expressed as(13)L=Ju,Φ+au,v,Φ−lv,Φ+λGΦ−fV0,where λ is the Lagrangian multiplier.The material derivative is used to obtain the shape derivative ofL. For more detailed theoretical proofs, the readers can refer to [20, 21]. The shape derivative of L can be obtained by(14)∂Lu,Φ∂t=∑i=1N∫Dβu,v,Φwi∂αi∂tdΩ+λ∑i=1N∫DδΦ∂αi∂tdΩ,and(15)βu,v,Φ=fu,u−εuCεv+pv+∇cv•n+κcv,where κ=divn=∇•−∇Φ/∇Φ is the mean curvature.The derivative of LagrangianL with respect to t can also be obtained by the chain rule as follows:(16)∂Lu,Φ∂t=∑i=1N∂Ju,Φ∂αi+λ∑i=1N∂GΦ∂αi∂αi∂t.By comparing the corresponding terms in Equations (14) and (16), the derivatives of the objective function and the volume constraint with respect to the interpolated coefficients can be obtained by(17)∂Ju,Φ∂αi=∫Dβu,v,ΦwidΩ,(18)∂GΦ∂αi=∫DδΦdΩ.When Equations (17) and (18) are obtained, the optimization problem in Equation (8) can be solved by sensitivity-based optimization algorithms.


## 3.3. Optimization Procedure

In optimization problem, with the sensitivities of objective function and volume constraint with respect to interpolated coefficients obtained, the sensitivity-based optimization algorithms such as optimality criteria and mathematical programming can be applied straightforwardly. Then, the shape and topology optimization problem converts into the parameter optimization problem. In this work, the method of moving asymptotes (MMA) which is a widely used mathematical programming algorithm is employed to solve the optimization problem. The process of optimization is shown in Figure1, and the corresponding introduction is presented as follows:Step 1: initialize the design domain of arch bridge. Define the load, boundary conditions, and set the volume constraint.Step 2: FE analysis and calculate the real and virtual displacement fields.Step 3: calculate the derivatives of objective function and the volume constraint with respect to the interpolated coefficients in Equations (17) and (18).Step 4: calculate the values of objective function and the volume constraint. Update the interpolated coefficients with MMA optimization algorithm, and then, the update level set function.Step 5: check whether the convergence condition is shown in Equation (19). If the terminal condition is met, stop, and the optimal structure is found. Otherwise, go back to Step 2. The terminal condition is defined as(19)Js−Js−1Js≤ε,where Js denotes the value of the objective function at sth iteration, and ε is a small positive value.Figure 1
Flowchart of form finding of arch bridge.


## 4. Numerical Examples

In this section, form finding of three kinds of arch bridges is presented to illustrate the effectiveness of the proposed method. The conventional level set method (CLSM), BESO, and SIMP are applied for comparison. In these examples, the computational points in PLSM are corresponding to the finite element nodes and four-node plane stress elements are used for finite element analysis (FEA). The parameters in Equation (8) are α¯k=6×maxαik0 and α¯k=6×minαik0, where αik0 denotes the initial values of interpolation coefficients, and the support radius for computational points is set to be R=2.5. Δ=1,χ=0.001, and γ=0.0005 are the parameters in Equations (10) and (11). The parameters in BESO are ER=2%,ARmax=50%, and rmin=2.5, while in SIMP, rmin=2.5. The Lenovo personal computer equipped with 1 Intel@ CoreTM i5-3230M processor with a clock speed of 2.6 GHz and 4.00 GB RAM is employed to calculate central processing unit (CPU) times.Example 1.
In this deck arch bridge design example, the admissible design domain, boundary conditions, and load conditions are shown in Figure2. The design domain is modeled into 180×60 elements, a uniform load F=0.01 is applied in the top of the design domain, and the prescribed volume fraction is f=0.3.
The optimal layouts obtained by PLSM, CLSM, BESO, and SIMP are shown in Figures3(a)–3(d), respectively, and it can be seen that the optimal layouts of different methods are almost identical and they are similar to the actual layout of deck arch bridge, which shows the effectiveness of PLSM. However, they also have some differences. The optimal layouts of PLSM and CLSM in Figures 3(a) and 3(b) have smooth boundaries, which are easier to transform into CAD model for bridge design. Figure 3(c) shows the optimal layout of BESO with the zigzag boundary. It is hard for follow-up design. Figure 3(d) shows the optimal layout of SIMP. There are a lot of gray elements, and the bars have split forks in the junction point of the bars and the bridge deck, which shows the optimal transmission path of force. In addition, the relatively smooth boundary can be produced by BESO and SIMP when the design space is modeled into a large number of finite elements. In this way, the computational efficiency is greatly reduced.
Figure4 gives the iteration history of the objective function and several intermediate results. The value of the objective function is 25.2 after 63 iterations. Corresponding to the PLSM topological designs in Figure 4, the level set surfaces at different iterations are displayed in Figure 5. Table 1 provides the results of different methods in form finding of deck arch bridge for comparison. It can be found that the compliance of PLSM is almost the same with that of other methods. And BESO produces the lowest compliance of 24.6 as it utilizes discrete variables. The PLSM costs the shortest computation time of 55.3 s which shows the high computational efficiency.Figure 2
Design space of the deck arch bridge.Figure 3
The optimal layout of deck arch bridge by different methods: (a) PLSM, (b) CLSM, (c) BESO, and (d) SIMP.
(a)
(b)
(c)
(d)Figure 4
Iteration history of structural compliance for form finding of deck arch bridge.Figure 5
Level set surfaces at different iterations.
(a)
(b)
(c)
(d)Table 1
Comparison of results of different methods in form finding of deck arch bridge.
PLSM
CLSM
BESO
SIMP
Compliance
25.2
25.7
24.6
25.9
Computation time(s)
85.3
284.5
253.1
193.2Example 2.
In this half-through arch bridge design example, the admissible design domain, boundary conditions, and load conditions are shown in Figure6. The design domain is modeled into 180×120 elements, a uniform load F=0.01 is applied in the height of h=30, and the prescribed volume fraction is f=0.2.
Figure7 shows the optimal layouts obtained by different optimization methods. It can be seen that the optimal layouts of different methods are similar to the actual layout of half-through arch bridge. The optimal configuration in PLSM is almost identical with the optimal configuration in CLSM as they have the same number of bars, which shows the validity of PLSM. Both of them have smooth boundaries. The optimal configurations of BESO and SIMP have more bars than that of PLSM and CLSM. The optimal configuration obtained by BESO in Figure 7(c) has inclined supports under the deck and also has a zigzag boundary. The configuration of SIMP in Figure 7(d) has a lot of gray elements.
The iteration history of the objective function and several intermediate results are shown in Figure8. The objective function is convergent to 19.9 after 37 iterations. For comparison, the results of different methods in form finding of deck arch bridge are presented in Table 2. It can be found that the compliance of PLSM is less than that of CLSM and SIMP, and BESO produces the lowest compliance of 18.6. The PLSM costs the shortest computation time of 100.9 s which is much less than other methods in calculation.Figure 6
Design space of the half-through arch bridge.Figure 7
The optimal layout of half-through arch bridge by different methods: (a) PLSM, (b) CLSM, (c) BESO, and (d) SIMP.
(a)
(b)
(c)
(d)Figure 8
Iteration history of structural compliance for form finding of half-through arch bridge.Table 2
Comparison of results of different methods in form finding of half-through arch bridge.
PLSM
CLSM
BESO
SIMP
Compliance
19.9
20.7
18.6
20.1
Computation time(s)
100.9
451.7
603.6
306.5Example 3.
In this through arch bridge design example, the admissible design domain, boundary conditions, and load conditions are shown in Figure9. The design domain is modeled into 180×120 elements, a uniform load F=0.01 is applied in the bottom of the design domain, and the prescribed volume fraction is f=0.2.
The optimal results of various topology optimization methods are illustrated in Figures10(a)–10(d), respectively. It can be seen that the optimal layouts of different methods are similar to the actual layout of through arch bridge. The optimal configurations in PLSM and CLSM have smooth boundaries, and the optimal configuration in PLSM is almost identical with the optimal configuration in CLSM. So, the effectiveness of the PLSM can be verified. The optimal configurations of BESO and SIMP have more bars than that of PLSM and CLSM. The configuration of BESO in Figure 10(c) has a zigzag boundary. The configuration of SIMP in Figure 10(c) has a lot of gray elements. And, it should be noticed that the bars have split forks in the junction point of the bars and the bridge deck.
Figure11 gives the iteration history of the objective function and several intermediate results. The value of objective function is 24.3 after 80 iterations. Table 3 summarizes the structural compliances and computation time of different methods in form finding of through arch bridge. It can be found that the compliance of PLSM is almost the same with that of other methods, and BESO produces the lowest compliance of 24.6. The PLSM costs the shortest computation time of 186.2 s. So, the accuracy and efficiency of PLSM are proved.Figure 9
Design space of the through arch bridge.Figure 10
The optimal layout of through arch bridge by different methods: (a) PLSM, (b) CLSM, (c) BESO, and (d) SIMP.
(a)
(b)
(c)
(d)Figure 11
Iteration history of structural compliance for form finding of through arch bridge.Table 3
Comparison of results of different methods in form finding of through arch bridge.
PLSM
CLSM
BESO
SIMP
Compliance
24.3
23.2
22.4
25.2
Computation time(s)
186.2
432.7
503.7
389.9


## 5. Conclusions

In this study, a form finding method for arch bridge is developed using the PLSM-based topology optimization, the solution procedure of the topology optimization method is described, and three numerical examples are employed to illustrate the feasibility and validity of PLSM in form finding of arch bridges. Numerical examples indicate that the configuration obtained by PLSM has smooth boundary. It is easier to convert the configuration to CAD model for manufacturing. In addition, the computational efficiency has been greatly improved.This method also can be extended to form finding of arch bridges with different rise-to-span, main tower of cable stayed bridge, and main cable of suspension bridge. Due to the effectiveness in form finding of arch bridges, it is expected that this method will be applied in bridge design in the future.



---

*Source: 2198696-2018-11-11.xml*