Experimental and Numerical Studies on a Glubam Spherical Dome

Ruijia WU, Ke MA, Pengyu LI, Yubing HOU, Peixiang WANG, Binbin LI, Yan XIAO, Experimental and numerical studies on a glubam spherical dome, Engineering Structures, Volume 303, 2024, 117462, ISSN 0141-0296, https://doi.org/10.1016/j.engstruct.2024.117462.


Abstract: This paper explores static experiments and dynamic tests, along with structural analysis of certain glubam (glued laminated bamboo) dome models. A self-balancing lever system was employed for the static experiment, revealing mechanical performance details such as load-displacement relationships and internal axial forces in glubam dome structures. Dynamic behaviors were assessed through triaxial accelerometers on specific joints, employing the Bayesian FFT modal identification method to calculate structural modal parameters, including natural frequency and mode shape. Experimental outcomes closely aligned with results from the analytical approach using the finite element method. Overall, the paper contributes valuable insights into the static and dynamic characteristics of glubam dome structures, bridging experimental findings with analytical predictions. Keywords: Bamboo-based dome; Static experiment; Dynamic test; Finite element analysis

Design a Space Elvator

TV series Foundation shows a space elevator 911 event. How should a space elevator be designed and what will its collapse look alike?

Consider the space elevator to be a link with linear mass \hat{m}(l) and a total length L.

To have a straight link, everything below the Geosynchronous orbit, or GSO would be rotating too slow to support its own mass. To balance the effect one can choose from 2 sources of support: base reaction from the ground, or tension from above GSO.

The tension of link at height x

T(x) = - F_{r} - \int_0^x{\cfrac{GM\hat{m}}{(R_{e}+l)^2}}dl + \int_0^L{\hat{m}\omega_{e}^2(R_{e}+l)}dl

where R_e is earth radius, \omega_e is earth rotating angular speed, F_r is base reaction force.

The boundry condition

T(L) = 0

Intuitively, a reasonable model will have the ground part in compression, but not that much of the entire gravity of things below synchronous orbit. The link at synchronous orbit will be in tension, so somewhere below there would be a zero force point.

The major payload should be deployed right at GSO, such that it doesn't pose any extra load into the link. Another benefit of the design is that the major payload can always survive any 911 attacks not targeted at it directly by simply cutting itself off the link. In such design, cut the link anywhere will have the lower part fall to Earth and upper part go into space. The TV series made a reasonable assumption.

Though seemingly unavoidable, it is by design inevitable for a (common civillian) building to fail a 911 attack due to economic concern, even a space elevator.

Axial force diagram should look like moment diagram of a cantilever beam under distributed load, and so should a reasonable section design be alike.

Project: Prefabricated Engineered Bamboo Mobile House – BIM Modeling

Glubam is short for glue laminated bamboo, which is a new material1.

This is the BIM modeling project for a new code: Technical specification for prefabricated engineered bamboo mobile house.

  1. Wu, Y., and Y. Xiao. "Steel and glubam hybrid space truss." Engineering Structures 171 (2018): 140-153. 

The 2nd Civil Engineering Computing and Simulation Technology Academic Conference

I would not put my half organized notes here.

I met some of the greatest Chinese scholars in the field. Watched their report, and found myself at a very strange position.

I love computing and simulation technology, but the deeper it goes the more clear its boundry is - technology is technology, civil engineering is engineering, not science.

Changsha Trip

True Fans Changsha Culture and Art Center, by Zaha Hadid.

A true fans sees him everywhere.

20210509 College of Civil Engineering, Hunan University

Odd enough, the metro exit reads 'College of Civil Engineering, Hunan University' instead of just 'Hunan University', despite the entire campus is just outside of the exit and College of Civil Engineering is the one most far from the exit.

Core Concepts in FEA


Before we start, some basic mathematic theorem12.


∇·v = v_{i,i}
∇·T = T_{ij,i} e_j

Gradient (dim+1)

∇\phi = \phi_{,i} e_i
∇v = v_{i,j} e_i⊗e_j
∇T = T_{ij,k} e_i⊗e_j⊗e_k

Curl (dim=)

∇ × v =ε_{ijk} v_{j,i} e_k
∇ × T =ε_{ijk} T_{mj,i} e_k⊗e_m

The product rule of differentiation

\because \bm{σ}·\bm{\nu} = \bm{σ}_{ij}\bm{\nu}_j e_i
\therefore ∇·(\bm{σ}·\bm{\nu}) = (\bm{σ}_{ij}\bm{\nu}_j)_{,i} = \bm{σ}_{ij,i}\bm{\nu}_j + \bm{σ}_{ij}\bm{\nu}_{j,i}


(∇·\bm{σ})·\bm{\nu} = \bm{σ}_{ij,i} \bm{\nu}_j
∇\bm{\nu} : \bm{σ} = \bm{\nu}_{i,j}\bm{σ}_{ji}

or as Prof. Zohdi suggested

∇\bm{\nu} : \bm{σ} = \bm{\nu}_{i,j}\bm{σ}_{ij}

but since \bm{σ} is symmetric

\bm{σ}_{ij} = \bm{σ}_{ji}


\tag{1} ∇·(\bm{σ}·\bm{\nu}) = (∇·\bm{σ})·\bm{\nu} + ∇\bm{\nu} : \bm{σ}

Divergence Theorem

∫_{∂R} \phi n dA = ∫_{R} ∇·\phi dV
\tag{2} ∫_{∂R} \bm{v}·n dA = ∫_{R} ∇·\bm{v} dV
∫_{∂R} \bm{T} n dA = ∫_{R} ∇·\bm{T} dV


∫_{∂R} \phi n_i dA = ∫_{R} \phi_{,i} dV
∫_{∂R} v_i n_i dA = ∫_{R} v_{i,i} dV
∫_{∂R} T_{ij}n_j dA = ∫_{R} T_{ij,j} dV

Stokes theorem

∫_{C} \phi d\bm{x} = ∫_{S} n × ∇\phi dA
∫_{C} v·d\bm{x} = ∫_{S} n·∇×v dA
∫_{C} T d\bm{x} = ∫_{S} (∇×T)^T n dA


∫_{C} \phi dx_i = ∫_{S} ε_{ijk} n_j \phi_{,k} dA
∫_{C} v_i dx_i = ∫_{S} n_i ε_{ijk} v_{k,j} dA
∫_{C} T_{ij} dx_j = ∫_{S} ε_{jpq} T_{iq,p} n_j dA


The SMALL deformation of the body is governed by (strong form):



\bm{σ} = \bm{IE}:∇\bm{u}


\bm{f} = ρ\bm{g}


∇·\bm{σ} + \bm{f} = 0

so we can write

∫_Ω (∇·\bm{σ} + \bm{f} ) · \bm{\nu} dΩ = 0

using (1)

∫_Ω (∇·(\bm{σ}·\bm{\nu})-∇\bm{\nu}:\bm{σ}) dΩ + ∫_Ω \bm{f}·\bm{\nu} dΩ = 0

using (2)

∫_Ω ∇\bm{\nu}:\bm{σ} dΩ = ∫_Ω \bm{f}·\bm{\nu} dΩ + ∫_{∂Ω} \bm{σ}·\bm{\nu}·n dA


\bm{t} = \bm{σ}·\bm{n}


∫_Ω ∇\bm{\nu}:\bm{σ} dΩ = ∫_Ω \bm{f}·\bm{\nu} dΩ + ∫_{Γ_t} \bm{t}·\bm{\nu} dA

so we have the weak form

\begin{array}{lcr} \text{Find }\bm{u}, \bm{u}|_{Γ_u} = \bm{u}^* \text{, such that }∀\bm{\nu}, \bm{\nu}|_{Γ_u}= 0 \\ \displaystyle∫_Ω ∇\bm{\nu}:\bm{IE}:∇\bm{u} dΩ = ∫_Ω \bm{f}·\bm{\nu} dΩ + ∫_{Γ_t} \bm{t}·\bm{\nu} dA \end{array}

where u^ is the applied boundary displacement on Γ_u, t = t^ on Γ_t

since the integrals must be finite, we improve the weak form as below

\begin{array}{lcr} \text{Find }\bm{u} ∈ \bm{H}^1(Ω), \bm{u}|_{Γ_u} = u^* \text{, such that }∀\bm{\nu} ∈ \bm{H}^1(Ω), \bm{\nu}|_{Γ_u}= 0 \\ \displaystyle∫_Ω ∇\bm{\nu}:\bm{IE}:∇\bm{u} dΩ = ∫_Ω \bm{f}·\bm{\nu} dΩ + ∫_{Γ_t} \bm{t}·\bm{\nu} dA \end{array}


\bm{H}^1(Ω) \vcentcolon= [H^1(Ω)]^3

explained as below

similar to the definition of u ∈ H^1(Ω) which states

u ∈ H^1(Ω) \text{ if } ||u||_{H^1(Ω)}^2 \vcentcolon= ∫_Ω u_{,j}u_{,j} dΩ + ∫_Ω uudΩ

the definition of \bm{u} ∈ \bm{H}^1(Ω) states as below

\bm{u} ∈ \bm{H}^1(Ω) \text{ if } ||\bm{u}||_{\bm{H}^1(Ω)}^2 \vcentcolon= ∫_Ω u_{i,j}u_{i,j} dΩ + ∫_Ω u_iu_idΩ


should \bm{u}^* be different from the applied boundary displacement on Γ_u, weak form can be stated as below

\tag{3} \begin{array}{lcr} \text{Find }\bm{u} ∈ \bm{H}^1(Ω) \text{, such that }∀\bm{\nu} ∈ \bm{H}^1(Ω) \\ \displaystyle∫_Ω ∇\bm{\nu}:\bm{IE}:∇\bm{u} dΩ = ∫_Ω \bm{f}·\bm{\nu} dΩ + ∫_{Γ_t} \bm{t}·\bm{\nu} dA + P^\star∫_{Γ_u} (\bm{u}^*-\bm{u})·\bm{\nu} dA \end{array}


\begin{array}{lcr} \text{Find }\bm{u} ∈ \bm{H}^1(Ω) \text{, such that }∀\bm{\nu} ∈ \bm{H}^1(Ω) \\ \displaystyle∫_Ω \nu_{i,j} IE_{ijkl} u_{k,l} dΩ = ∫_Ω f_i \nu_i dΩ + ∫_{Γ_t} t_i \nu_i dA + P^\star∫_{Γ_u} u^*_i \nu_i - u_i \nu_i dA \end{array}

where the (penalty) parameter P^\star is a large positive number.


To have

\begin{cases} \hat{\phi}_i(\zeta_1,\zeta_2,\zeta_3) = 1 &\text{ , } \zeta_1=\zeta_{i1},\zeta_2=\zeta_{i2},\zeta_3=\zeta_{i3} \\ \hat{\phi}_i(\zeta_1,\zeta_2,\zeta_3) = 0 &\text{ , } \text{others} \end{cases}

for each \hat{\phi}_i 8 functions (for each of the 8 points) with 8 unknowns, aka. \zeta_1^m \zeta_2^n \zeta_3^k where m,n,k=0,1 can be derived, leading to the only solution as following

\begin{cases} \hat{\phi}_1 = \frac{1}{8} (1-\zeta_1) (1-\zeta_2) (1-\zeta_3) \\ \hat{\phi}_2 = \frac{1}{8} (1+\zeta_1) (1-\zeta_2) (1-\zeta_3) \\ \hat{\phi}_3 = \frac{1}{8} (1+\zeta_1) (1+\zeta_2) (1-\zeta_3) \\ \hat{\phi}_4 = \frac{1}{8} (1-\zeta_1) (1+\zeta_2) (1-\zeta_3) \\ \hat{\phi}_5 = \frac{1}{8} (1-\zeta_1) (1-\zeta_2) (1+\zeta_3) \\ \hat{\phi}_6 = \frac{1}{8} (1+\zeta_1) (1-\zeta_2) (1+\zeta_3) \\ \hat{\phi}_7 = \frac{1}{8} (1+\zeta_1) (1+\zeta_2) (1+\zeta_3) \\ \hat{\phi}_8 = \frac{1}{8} (1-\zeta_1) (1+\zeta_2) (1+\zeta_3) \end{cases}


\bm{IE} has the following symmetries 23

IE_{ijkl} = IE_{klij} = IE_{jikl} = IE_{ijlk}

allowing us to rewrite (3) in a more compact matrix form

\tag{4} \begin{aligned} ∫_Ω ([\bm{D}]\{\bm{\nu}\})^T [\bm{IE}] ([\bm{D}]\{\bm{u}\}) dΩ = &∫_Ω \{\bm{\nu}\}^T\{\bm{f}\} dΩ \\ &+ ∫_{Γ_t} \{\bm{\nu}\}^T\{\bm{t}^*\} dA \\ &+ P^\star∫_{Γ_u} \{\bm{\nu}\}^T\{\bm{u}^*-\bm{u}\} dA \end{aligned}


\tag{5} [\bm{D}] \vcentcolon = \begin{bmatrix} \cfrac{∂}{∂x_1} & 0 & 0 \\ 0 & \cfrac{∂}{∂x_2} & 0 \\ 0 & 0 & \cfrac{∂}{∂x_3} \\ \cfrac{∂}{∂x_2} & \cfrac{∂}{∂x_1} & 0 \\ 0 & \cfrac{∂}{∂x_3} & \cfrac{∂}{∂x_2} \\ \cfrac{∂}{∂x_3} & 0 & \cfrac{∂}{∂x_1} \\ \end{bmatrix}, \{\bm{u}\} \vcentcolon = \begin{Bmatrix} u_1 \\ u_2 \\ u_3 \\ \end{Bmatrix}, \{\bm{f}\} \vcentcolon = \begin{Bmatrix} f_1 \\ f_2 \\ f_3 \\ \end{Bmatrix}, \{\bm{t}^*\} \vcentcolon = \begin{Bmatrix} t^*_1 \\ t^*_2 \\ t^*_3 \\ \end{Bmatrix},

and [\bm{IE}] is the elastic stiffness matrix of the material.

[\bm{IE}] \vcentcolon = \begin{bmatrix} IE_{1111} & IE_{1122} & IE_{1133} & IE_{1112} & IE_{1123} & IE_{1113} \\ IE_{2211} & IE_{2222} & IE_{2233} & IE_{2212} & IE_{2223} & IE_{2213} \\ IE_{3311} & IE_{3322} & IE_{3333} & IE_{3312} & IE_{3323} & IE_{3313} \\ IE_{1211} & IE_{1222} & IE_{1233} & IE_{1212} & IE_{1223} & IE_{1213} \\ IE_{2311} & IE_{2322} & IE_{2333} & IE_{2312} & IE_{2323} & IE_{2313} \\ IE_{1311} & IE_{1322} & IE_{1333} & IE_{1312} & IE_{1323} & IE_{1313} \\ \end{bmatrix}

Further rewrite {\bm{u}^h}

\tag{6} \{\bm{u}^h\} = [\phi]\{\bm{a}\}


\{\bm{a}\} \vcentcolon = \begin{Bmatrix} a_1 \\ a_2 \\ a_3 \\ .\\ .\\ .\\ a_{3N}\\ \end{Bmatrix}, [\phi] \vcentcolon = \begin{bmatrix} \phi_1 & \phi_2 & ... & \phi_N & & & & & & & & \\ & & & & \phi_1 & \phi_2 & ... & \phi_N & & & & \\ & & & & & & & & \phi_1 & \phi_2 & ... & \phi_N \\ \end{bmatrix},

choose \bm{\nu} with the same basis, but a different linear combination

\tag{7} \{\bm{\nu}^h\} = [\phi]\{\bm{b}\}

with (6),(7) we can rewrite (4)

\tag{8} \begin{aligned} ∫_Ω ([\bm{D}][\phi]\{\bm{b}\})^T [\bm{IE}] ([\bm{D}][\phi]\{\bm{a}\}) dΩ = &∫_Ω ([\phi]\{\bm{b}\})^T\{\bm{f}\} dΩ \\ & + ∫_{Γ_t} ([\phi]\{\bm{b}\})^T\{\bm{t}^*\} dA \\ & + P^\star∫_{Γ_u} ([\phi]\{\bm{b}\})^T\{\bm{u}^*-[\phi]\{\bm{a}\}\} dA \end{aligned}

We can rewrite (8)

\{\bm{b}\}^T \{ [\bm{K}] \{\bm{a}\} - \{\bm{R}\} \} = 0


\tag{9} [\bm{K}] \vcentcolon = ∫_Ω ([\bm{D}][\phi])^T [\bm{IE}] ([\bm{D}][\phi]) dΩ + P^\star∫_{Γ_u} [\phi]^T[\phi] dA
\tag{10} \{\bm{R}\} \vcentcolon = ∫_Ω [\phi]^T\{\bm{f}\} dΩ + ∫_{Γ_t} [\phi]^T\{\bm{t}^*\} dA + P^\star∫_{Γ_u} [\phi]^T\{\bm{u}^*\} dA


[\bm{K}] = \sum_{e=1}^{N_e}{[\bm{K}]^e}

we can rewrite (9),(10) for e=1,2..,N_e

\tag{11} [\bm{K}]^e \vcentcolon = ∫_{Ω_e} ([\bm{D}][\phi])^T [\bm{IE}] ([\bm{D}][\phi]) dΩ_e + P^\star∫_{Γ_{u,e}} [\phi]^T[\phi] dA_e
\tag{12} \{\bm{R}\}^e \vcentcolon = ∫_{Ω_e} [\phi]^T\{\bm{f}\} dΩ_e + ∫_{Γ_{t,e}} [\phi]^T\{\bm{t}^*\} dA_e + P^\star∫_{Γ_{u,e}} [\phi]^T\{\bm{u}^*\} dA_e


Γ_{u,e} = Γ_u ∩ ∂Ω_e
Γ_{t,e} = Γ_t ∩ ∂Ω_e


In 3D,

[\hat\phi] \vcentcolon = \begin{bmatrix} \hat\phi_1 & \hat\phi_2 & ... & \hat\phi_8 & & & & & & & & \\ & & & & \hat\phi_1 & \hat\phi_2 & ... & \hat\phi_8 & & & & \\ & & & & & & & & \hat\phi_1 & \hat\phi_2 & ... & \hat\phi_8 \\ \end{bmatrix}

express [\bm{D}] in terms ζ_1, ζ_2, ζ_3

[D(\phi(x1, x2, x3))] = [\hat D \phi(M_{x_1}(ζ1, ζ2, ζ3),M_{x_2}(ζ1, ζ2, ζ3),M_{x_3}(ζ1, ζ2, ζ3))]

so that

[\hat \bm{D}] = \begin{bmatrix} \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_1} & 0 & 0 \\ 0 & \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_2} & 0 \\ 0 & 0 & \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_3} \\ \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_2} & \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_1} & 0 \\ 0 & \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_3} & \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_2} \\ \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_3} & 0 & \cfrac{∂}{ζ_i}\cfrac{ζ_i}{∂x_1} \\ \end{bmatrix}


[\hat \bm{D}][\hat\phi] = \begin{bmatrix} \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_1} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_1} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_1} & & & & & & & \\ & & & & \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_2} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_2} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_2} & & & \\ & & & & & & & & \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_3} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_3} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_3} \\ \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_2} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_2} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_2} & \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_1} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_1} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_1} & & & \\ & & & & \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_3} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_3} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_3} & \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_2} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_2} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_2} \\ \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_3} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_3} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_3} & & & & & \cfrac{∂\hat\phi_1}{ζ_i}\cfrac{ζ_i}{∂x_1} & \cfrac{∂\hat\phi_2}{ζ_i}\cfrac{ζ_i}{∂x_1} & ... & \cfrac{∂\hat\phi_8}{ζ_i}\cfrac{ζ_i}{∂x_1} \\ \end{bmatrix}

with Gaussian Quadrature

\int_0^L{F(x)dx} = \int_{-1}^1{F(x(\zeta))J(\zeta)d\zeta} = \sum_1^G{w_i}F(\zeta_i)J(\zeta_i)

we can rewrite (11),(12)

[\bm{K}]^e = \sum_1^G{\sum_1^G{\sum_1^G{w_q w_r w_s ([\hat \bm{D}][\hat\phi])^T[\hat{\bm{IE}}][\hat \bm{D}][\hat\phi]J}}}+ P^\star\sum_1^G{\sum_1^G{w_q w_r [\phi]^T[\phi] J_s}}
\{\bm{R}\}^e = \sum_1^G{\sum_1^G{\sum_1^G{w_q w_r w_s [\phi]^T\{\bm{f}\} J}}} + \sum_1^G{\sum_1^G{w_q w_r[\phi]^T\{\bm{t}^*\} J_s}} + P^\star\sum_1^G{\sum_1^G{w_q w_r [\phi]^T\{\bm{u}^*\} J_s}}


\bm{F} \vcentcolon= ∇x(ζ_1,ζ_2,ζ_3)
J \vcentcolon= |\bm{F}|


J_s = J\bm{F}^{-T}·\bm{N}·\bm{n}

derived as below

Nanson’s formula

\bm{n} dA_e = J\bm{F}^{-T}·\bm{N} d\hat{Ae}

where A_e is an area of a region in the current configuration, \hat{Ae} is the same area in the reference configuration, and \bm{n} is the outward normal to the area element in the current configuration while \bm{N} is the outward normal in the reference configuration.

multiply \bm{n} on both sides

dA_e = J\bm{F}^{-T}·\bm{N}·\bm{n} d\hat{Ae}


J_s = J\bm{F}^{-T}·\bm{N}·\bm{n}




let the heat flux density4 q

\bm{q} = \bm{IK}⋅∇θ

and let

f = ρs


∇⋅\bm{q} + f = 0

so we can write

\tag{13} ∫_Ω (∇·\bm{q} + f ) · \nu dΩ = 0

now consider

(∇⋅\bm{q})·\nu = q_{i,i}\nu
∇⋅(\bm{q}\nu) = (q_i \nu)_{,i} = q_{i,i}\nu + q_i \nu_{,i}
(∇\nu)·\bm{q} = q_i \nu_{,i}

therefore we can write (13) as

∫_Ω (∇⋅(\bm{q}\nu) - (∇\nu)·\bm{q}) dΩ + ∫_Ω f\nu dΩ= 0
∫_Ω (∇\nu)·\bm{q} dΩ = ∫_Ω f\nu dΩ + ∫_Ω ∇⋅(\bm{q}\nu)

using (2)

∫_Ω (∇\nu)·\bm{q} dΩ = ∫_Ω f\nu dΩ + ∫_{∂Ω} \bm{q}⋅\bm{n} \nu dA

so we have the weak form

\begin{array}{lcr} \text{Find }θ ∈ H^1(Ω), θ|_{Γ_θ} = θ^* \text{, such that }∀\nu ∈ H^1(Ω), \nu|_{Γ_θ}= 0 \\ \displaystyle∫_Ω (∇\nu)·\bm{IK}⋅∇θ dΩ = ∫_Ω f \nu dΩ + ∫_{∂Ω} \bm{q}⋅\bm{n} \nu dA \end{array}


should θ^* be different from the applied boundary displacement on Γ_θ, weak form can be stated as below

\tag{14} \begin{array}{lcr} \text{Find }θ ∈ H^1(Ω), θ|_{Γ_θ} = θ^* \text{, such that }∀\nu ∈ H^1(Ω), \nu|_{Γ_θ}= 0 \\ \displaystyle∫_Ω (∇\nu)·\bm{IK}⋅∇θ dΩ = ∫_Ω f \nu dΩ + ∫_{∂Ω} \bm{q}⋅\bm{n} \nu dA + P^\star∫_{Γ_θ} (θ^*-θ) \nu dA \end{array}

where the (penalty) parameter P^\star is a large positive number.


same as 3


rewrite θ^h

\tag{15} θ^h = [\phi]\{\bm{a}\} = \phi_i a_i


\{\bm{a}\} \vcentcolon = \begin{Bmatrix} a_1 \\ a_2 \\ a_3 \\ .\\ .\\ .\\ a_{N}\\ \end{Bmatrix}, [\phi] \vcentcolon = \begin{bmatrix} \phi_1 & \phi_2 & ... & \phi_N \end{bmatrix},

choose \bm{\nu} with the same basis, but a different linear combination

\tag{16} \nu^h = [\phi]\{\bm{b}\} = \phi_j b_j

with (15),(16) we can rewrite (14)

\tag{17} \begin{aligned} ∫_Ω (∇[\phi]\{\bm{b}\})·\bm{IK}⋅∇[\phi]\{\bm{a}\} dΩ = & ∫_Ω f [\phi]\{\bm{b}\} dΩ \\ & + ∫_{∂Ω} \bm{q}⋅\bm{n} [\phi]\{\bm{b}\} dA \\ & + P^\star∫_{Γ_θ} (θ^*-[\phi]\{\bm{a}\})·[\phi]\{\bm{b}\} dA \end{aligned}

We can rewrite (17)

\{\bm{b}\}^T \{ K \{\bm{a}\} - \{\bm{R}\} \} = 0


\tag{18} K \vcentcolon = ∫_Ω (∇[\phi])·\bm{IK}⋅∇[\phi] dΩ + P^\star∫_{Γ_θ} [\phi]^T[\phi] dA
\tag{19} \{\bm{R}\} \vcentcolon = ∫_Ω [\phi]^Tf dΩ + ∫_{∂Ω} [\phi]^T\bm{q}⋅\bm{n} dA + P^\star∫_{Γ_θ} [\phi]^Tθ^* dA


K= \sum_{e=1}^{N_e}K^e

we can rewrite (18),(19) for e=1,2..,N_e

K^e \vcentcolon = ∫_{Ω_e} (∇[\phi])·\bm{IK}⋅∇[\phi] dΩ_e + P^\star∫_{Γ_{u,e}} [\phi]^T[\phi] dA_e
\{\bm{R}\}^e \vcentcolon = ∫_{Ω_e} [\phi]^Tf dΩ_e + ∫_{∂Ω,e} [\phi]^T\bm{q}⋅\bm{n} dA_e + P^\star∫_{Γ_{u,e}} [\phi]^Tθ^* dA_e


Γ_{θ,e} = Γ_θ ∩ ∂Ω_e


Key difference is that the diemsion of the equations derived for thermodynamics are generally lower than those for linear elasticity.

This is because θ and s are scalers rather than vectors as \bm{u} and \bm{g}, and \bm{IK} is a second order tensor rather than a forth order tensor as \bm{IE}.


  1. Tensor derivative (continuum mechanics) Cartesian coordinates, https://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Cartesian_coordinates_2 

  2. Lallit Anand and Sanjay Govindjee, 2019, Introduction to Continuum Mechanics of Solid 

  3. Allan F. Bower, 2012, Applied Mechanics of Solids, Chapter 3 

  4. Thermal conduction differential form, https://en.wikipedia.org/wiki/Thermal_conduction#Differential_form 

Curve Generation Algorithm Development

UCB ARCH 259 Robotic Fabrication

Mark Ma, Yasaman Yavari


The final object is to make something fantastic, while creating a workflow including techniques of 3D scanning, parametric modeling, robotic arm operations, 3D printing and bio-mechanical knowledge.

While breaking down the final object, Yasaman and I found our interest in the parametric modeling part. With instructions from Professor Simon, we began our quest to develop a robot-printing friendly curve generation algorithm. To be robot-printing friendly, the curve generated should have the following features:

  • Artistic (subjective)
  • Reasonably structured (objective)
  • As continuous as possible so that saves the trouble of reconnecting during printing (objective)


Base shape

With the conclusion that we’re to use elastic canvas for basic form-finding, the first step of our modeling would be simulating this process. Several opinions are available:

  • Use finite element method to find the exact shape of the canvas under loading
  • Use Kangaroo to simulate the physical process
  • Use minimal surface

The minimal surface method is dropped due to the fact that we are to use form-finding in reality. The FEA method, despite exact and promising, does not go well with the rest of the project with its model format. Considering that the modeling precision would not be a big issue in this step, Yasaman and I decided to use Kangaroo for base shape simulation.

The drawback of Kangaroo is that we can hardly define an isotropic elastic material. The Quadrilateral Mesh, despite simple and straightforward, promises the same modulus only in the orthogonal direction of the mesh UV basis.

Through careful modeling, I add another restriction at 45° and 135° of the UV basis, making it closer to an isotropic material. Considering the exact model are to be found out using real form-finding and 3D-scanner, the remaining error is neglected.

20200501 Base Shape
20200501 Base Shape

Structural reasonableness

Considering the problem to be print on the basic form, several ideas of curve generation have been proposed. One plan is to adopt optimized principle stress lines. The advantage is obvious, as such a structure would be very reasonable under loading. However, principal stress lines turned out hard to be adopted. The continuity of the principle stress line does not guarantee the well-distribution of these lines, which often leaves large holes in structure generated. Due to the high sensitiveness to shape and load, the principal stress lines are hard to manipulate too. Consequently, we decided to start from the curve generation itself.

Curve generation

Step 1. Curve generated according to force distribution

One way of curve generation is to divide and lengthen parts of a continuous curve, while keeping the curve away from intersecting with itself, which defines all the goals we need to implement this idea with Kangaroo:

  • Lengthen curve sections
  • Prevent curve from intersecting with itself

To adapt the generated curves responsive to structural behavior, it's natural to densify the area with large stress, which gives the third goal:

  • Curve distribution responsive to stress distribution

With these goals the curves are successfully generated, meeting all our objective goals. Further development should focus on the subjective part, aka the artistic effects.

20200501 Step 1
20200501 Step 1

Step 2. Curve generated freely

To gain more control of the curve generation, some restrictions have to be dropped. Now take a look back at the initial three goals, if we neglect the curve distribution requirement, we can achieve something looks very alike Zaha’s project.

20200501 Step 2
20200501 Step 2

Recall that Zaha’s project uses a manually defined base curve to approximate initial stress distribution, we now certainly have a better solution: to use the curve generated in 1st step as the base shape for 2nd step, which gives us the following result.

20200501 Step 12
20200501 Step 12

Step 3. Combination

There is an infinite number of curves we can generate by tweaking with all the parameters, and an even larger infinite number of combinations with different layer position, tube thickness, color, etc. For a quick example, this is what it will look like if we combine Step 1+2 as a thicker base layer and Step 2 as a thinner decoration layer.

20200501 Step 3
20200501 Step 3


The curve generation is a complex workflow not fully automated, as lots of attention needs to be paid to calibrate all the parameters in order to get the best shape. All generated curves have the potential to become the final project, with a combination of each other and further improvement based on printing practice.

Due to the virus situation nowadays, it’s becoming harder for our project to be actually constructed. I hope all these quests into curve generation can be preserved, developed, and applied in future practice.


  1. Zaha Hadid Design - Thallus for for White in the City Animation - YouTube

  2. Curve growth - differential growth (Where to Start?) - Grasshopper - McNeel Forum

General parametric design of a steel-glubam hybrid space truss

Ma, K., and Xiao, Y. “General Parametric Design of a Steel-Glubam Hybrid Space Truss.” In Modern Engineered Bamboo Structures: Proceedings of the Third International Conference on Modern Bamboo Structures (ICBS 2018), June 25-27, 2018, Beijing, China, 1st ed., 223–29. CRC Press, 2019.

K. Ma

Zhejiang Univ.-Univ. of Illinois at Urbana Champaign Institute, Jiaxing, Zhejiang, China

Y. Xiao

Zhejiang Univ.-Univ. of Illinois at Urbana Champaign Institute, Jiaxing, Zhejiang, China

Nanjing Tech University, Nanjing, Jiangsu, China

Department of Civil Engineering, University of Southern California, Los Angeles, CA, USA

ABSTRACT: This paper introduces a parametric design method for a hybrid truss system composed of glued laminated bamboo (glubam) and steel. Experiments on determining material’s physical and mechanical parameters were carried out first, on basis of which design stages from modeling, analysis, optimization to manufacturing are all rendered possible through parametric ways by defining corresponding parameters within one single platform - Grasshopper. By maximizing automation during the process, efficiency and extensibility are taken into consideration for possibly further, larger, and more complex design.

Good Lab

A structural lab that is clean is not a good lab. -- Khalid M. Mosalam on CE 249, Sept. 6, 2019

edit 20220801:

The ZJU Haining campus recruited a new Civ Lab engineer recently, who is dedicated to make the Civ lab 'as clean and beautiful as a museum'. He once proposed that wood cutting should be done outside of the Lab to minimize the noise and ash, and I had to call the Dean in to revert this policy.

With his help everyting is now put in order far away from its supposed working location, just we can not find things when we need, and the lab engineer has no idea too.

How ironically!

Project: Large Waste Transfer Station

20190628 Large Waste Transfer Station 20190628 Large Waste Transfer Station

20190811 Large Waste Transfer Station 20200811 Large Waste Transfer Station

20190815 Large Waste Transfer Station 20200815 Large Waste Transfer Station


Directed by Prof. Yan Xiao.

Designed by Mark Ma.

Non-commercial project for research and demo purpose.

Project: Concrete Dragon

About the project

Concrete is something not often linked to boat, let alone dragon boat. Traditional concept of concrete based building material is embedded deeply not only in common folks but also in Civil Engineering students.

Concrete, usually Portland cement concrete (for its visual resemblance to Portland stone), is a composite material composed of fine and coarse aggregate bonded together with a fluid cement (cement paste) that hardens over time

Concrete - Wikipedia

Still, concrete is just a material, and material evolves as scientific progress advances. New material such as FRP rebars show even better performance than traditional steel rebars. Moreover, as its Chinese name indicated, concrete is just 'some sort of ash mixed and hardens', expanding its concept to a even higher dimension. Comparing to traditional boat building material, concrete has its own pros and cons.


  • Stain proof (v.s. steel)
  • Flexible shape (v.s. timber)
  • Economical (v.s. steel and timber)


  • Heavier than water and timber
  • Brittle
  • Low tensile strength

The whole project is aimed at exploiting the advantages of concrete while fixing its problems, expanding our knowledge of concrete through engineering practice.

Prototype 1.x - The First FRP Concrete Dragon Boat


A new type of boat construction method is proposed firstly by Prof. Yan Xiao. The basic idea is to use a self-hardening FRP cloth material as both concrete model and structural part that provides tensile strength.

20190615 Prototype 1.x Cross Section
20190615 Prototype 1.x Cross Section


FRP Clothes

This is some material recently become popular for its ability of free shaping and fast hardening. Most common scenario of usage is replacement of medical bondage and plaster, or electrical wire joint protection. The mechanism is really simple as it's just piece of FRP soaked into glue. Exposition to air solidify the glue in several minutes, providing the structure with some strength.

Things become interesting when several layers of such kind of material are glued together. As FRP clothes is extremely strong in tensile direction, even glue can sustain the relatively less strong shearing force, guaranteeing high bending and shearing strength.

Plastic Hull

A plastic hull is needed for this model, primarily for the shape control of the self-hardening FRP clothes. To achieve better dynamic performance, hull is modeled to be smooth and continuous - NURBS interposition of key points. Therefore, 3D-printing technology is adopted for precise realization of the model.

20190406 3D Printing
20190406 3D Printing
20190223 Modal Boat Hull
20190223 Modal Boat Hull
20190317 Modal Boat Covered by FRP with Kaihang Zhang
20190317 Modal Boat Covered by FRP with Kaihang Zhang


ECC concrete is a new kind of concrete which includes PVA fiber to increase its tensile strength and crack resistance. Ingredients are listed as below:

Cement 1
Quartz Sand 20-40 1.1
Silica Fume 0.3
Flyash 0.15
Mineral Powder 0.1
Quartz Sand 325 0.1
PVA Fiber
Water 0.18
Water Reduce Agents 0.10%
Early Strength Agents 0.10%

Next Stage

Prototype 1.x focuses itself on boat construction. A boat floating on the water marked success of the first stage. On that occasion, however, path splits for two different prototypes though.

Prototype 2.x stays small and is optimized for wireless remote control. A concrete dragon boat (model) competition is later held based on this. Prototype 3.x is develped with optimization for man-powered sailing in a much larger scale - a real concrete dragon boat.

Prototype 2.x - ICDBC 2019


Fast, stable, swift steering

To make it fast, the boat has to be speed boat alike, narrow, and light. The concrete layer needs to be as thin as possible. For sure, boat surface still needs to be smooth and continuous.

To make it stable and swift steering, the boat has to be not that narrow in width, motors parallel to each other and as far away as each other as possible.

Power system

Duo RH-380 motors, 1500mAh Li-ion battery

Should gave adopted stronger motors.




CentOS7, php7.2, MySQL5.3, Nginx, WordPress Vultr Dallas

See Concrete Dragon Online.


Should have set limitations on motor power.


Should have provided ai version.

Should have made it clear that theme color can be switched.

Should have designed certificate templates.

Next stage

Should have adopted Arduino.

20190606 ICDBC Boat of My Team
20190606 ICDBC Boat of My Team
20190606 ICDBC Jury Committee with Boats
20190606 ICDBC Jury Committee with Boats
20190606 ICDBC Prof. Yan Xiao
20190606 ICDBC Prof. Yan Xiao
20190606 ICDBC Everyone
20190606 ICDBC Everyone

Prototype 3.x - The Great Dragon Boat


Strong, stable, light

The large boat is very different from the small boat as construction method can not be inherited directly. It is too slow for the hull to be 3d-printed out. Instead, EPS foam slices cut by CO2 laser are glued (or tied) up to form the basic model of the boat. A FRP clothes layer is then put on the foam, then concrete with FRP bars inside.

20190526 Laser Cut Plastic Foam
20190526 Laser Cut Plastic Foam

To make the boat strong, the boat is modeled even more smooth, to the extent that FRP rebars can be put inside the concrete from head to tail without cut or connect. Such smooth shape without any edge can be regarded as a 2-d arch, making the most of the compressive strength of the concrete. To make the boat stable, the boat should have a lower gravity center.

20190615 Prototype 2.x Cross Section
20190615 Prototype 2.x Cross Section
20190531 Great Dragon Boat Foam and Dr. Dade Lai
20190531 Great Dragon Boat Foam and Dr. Dade Lai
20190601 Great Dragon Boat Covered by Concrete
20190601 Great Dragon Boat Covered by Concrete
20190604 Great Dragon Boat under Painting
20190604 Great Dragon Boat under Painting
20190606 Great Dragon Boat
20190606 Great Dragon Boat

Next stage

  • Should have the whole boat hull made by solid EPS form.
  • Should have adopted sprayed concrete technique.
  • Should have used more PVA.
  • Should have used more FRP fiber clothes.



Yan Xiao 肖岩

Boat Designer, Engineer, Manufacturer, Website Programmer and Designer

Ke Ma 马克


Ke Ma 马克, Dade Lai 赖大德, Zhiwei He 贺智伟

Also thanks for help by

Anqi Tan 谭安琪, Yiqi Feng 冯亦奇, Kaihang Zhang 张开航, Zhekai Li 李哲楷, Mengjun Wang 王梦君, Guoli Wang 汪郭立, Yang Zhou 周洋, Shangchun Jiang 江上春, Cristoforo Demartino, Zicheng Bao 包梓成

Decoration by

Jiahui Liang 梁嘉惠, Ke Ma 马克, Sicheng Zhou 周思成, Anqi Tan 谭安琪, Qingyun Liu 柳青云

Concrete Recipe

Bo Shan 单波, Dade Lai 赖大德

Logistic Support by

Sicheng Zhou 周思成

Competition Rule Maker

Yan Xiao 肖岩, Ke Ma 马克, Yiqi Feng 冯亦奇, Anqi Tan 谭安琪, Kaihang Zhang 张开航

Activity Organizer

Hang Wu 吴行, Ke Ma 马克, Binbin Li 李宾宾, Yanlong Xie 谢焱龙, Sichen Zhou 周思成, Tao Li 李涛, Jiahui Liang 梁嘉惠, Zhiwei He 贺智伟, Dade Lai 赖大德, Haixiang Zhu 朱海翔, Yi Zhang 张旖, Chenchen Ye 叶晨晨, Qian Yu 余倩, Jinyan Yu 俞静琰, Fengqing Jiang 江凤清 and other volunteers from ZJUI

Special Thanks for

Logistic Support

Jiyao Guo 郭霁瑶, Zhaijin Jia 贾翟菁, Qingbing Xie 谢庆兵

Technical Support

Tianyi Han 韩天屹, Si Li 李斯


Yue Feng 冯越, Hanyin Shao 邵寒吟



Project: Waste Transfer Station

20190324 Waste Transfer Station 1.0.1 Day 1 20190324 Waste Transfer Station, by Y. Xiao, K. Ma

202001030925 Waste Transfer Station 202001030925 Waste Transfer Station, Photo by S.C. Zhou


Directed by Prof. Yan Xiao.

Designed by Mark Ma.

Non-commercial project for research and demo purpose.

Employed as Solid Mechanics TA

20190114 Employed as TA of TAM 251
20190114 Employed as TA of TAM 251

http://tam251.intl.zju.edu.cn/people.html(LAN Access Only) It's been like a decade since I was designated as a class president or such sort of things. Still, I feel honoured. It's a fun fact that most of the TAM 251 website is edited by me (revised from UIUC), as years of playing with html, servers, etc. finally get me on the course 🙂 Last time I was soaked in with Solid Mechanics, I scored 94 out of 100. To me, it's almost the most interesting Mechanics subject I learned (comparing with structual mechanics and others). The point is, it's a theory-based subject for scientific research rather than an experience-based one for practical engineering. That means fun for me.

Project: RE Frame

Non-commercial project for research and demo purpose.

Project Introduction

This project is based on the 1st-prize-winning project 'RE Frame' of the Bamboo Pavilion Competition in Mediterranean University of Reggio Calabria by Stefano Vitale and his team. It's an honor to join him in such a great project. All following are based on my own understanding of the project.


马克. 基于Grasshopper的曲面木结构网壳参数化设计研究[D]:[学士学位论文].南京工业大学土木工程学院, 2018.









关键词:木网壳 参数化设计 Grasshopper

Voronoi 4.2
Voronoi 4.2



感谢Grasshopper的开发者David Rutten以及整个社区,特别是Paneling Tools, Kangaroo(2), LunchBox, NGon, Karamba等插件和他们的开发者。他们的工作为本课题的开展提供了巨大的便利,并在实际应用中带给了我一些思路上的启发。





[1] 关于公布南京工业大学2018届本科生优秀毕业设计(论文)、优秀指导教师的通知


Designed by Mark Ma.

Non-commercial project for research and demo purpose.

Produced with

  • Rhino 6, Educantional License
  • Lumion 6, Trial License
  • Adobe Creative Cloud, Educational License