On consistent nonlinear analysis of soil–structure interaction problems

The CPDM model can be extended to capture the viscoelastic creep and ageing based on the EC2 standard. To simplify the formulation, nonlinear creep effects that may appear for larger compressive stresses are neglected, and no distinction is made between the creep effects in compression and tension. In the EC2 standard, creep and ageing are represented by the creep coefficient ϕ(t, t_o), understood as a ratio between the creep and an instantaneous elastic strain computed for concrete loaded at t=28 days, and defined as follows:

ϕ(t,to)=ϕoβc(t,to)$$\begin{array}{} \displaystyle \phi(t,t_{o})=\phi_{o}\beta_{c}(t,t_{o}) \end{array}$$(2)

ϕo=ϕRHβ(fcm)β(to)$$\begin{array}{} \displaystyle \phi_{o}=\phi_{RH}\beta(f_{cm})\beta(t_{o}) \end{array}$$(3)

ϕRH=1+1−RH1000.1ho1/3forfcm≤35MPa1+1−RH1000.1ho1/3α1α2forfcm>35MPa$$\begin{array}{} \displaystyle \phi_{RH}= \left\{\begin{array}{}1+\frac{1-\frac{RH}{100}}{0.1\,h_{o}^{1/3}}\qquad\qquad\qquad~~\,\text{for}\quad f_{cm}\le35\text{MPa}\\ \displaystyle\left(1+\frac{1-\frac{RH}{100}}{0.1~h_{o}^{1/3}}\alpha_{1}\right)\alpha_{2}\qquad \text{for}\quad f_{cm}\gt 35\text{MPa} \end{array}\right\} \end{array}$$(4)

β(fcm)=16.8fcm$$\begin{array}{} \displaystyle \beta(f_{cm})=\frac{16.8}{\sqrt{f_{cm}}} \end{array}$$(5)

β(to)=10.1+to0.2$$\begin{array}{} \displaystyle \beta(t_{o})=\frac{1}{0.1+t_{o}^{0.2}} \end{array}$$(6)

βc(t,to)=[t−toβH+t−to]0.3$$\begin{array}{} \displaystyle \beta_{c}(t,t_{o})=\bigg[\frac{t-t_{o}}{\beta_{H}+t-t_{o}}\bigg]^{0.3} \end{array}$$(7)

βH=min(1.5[1+(0.012RH)18]ho+250,1500)forfcm≤35MPamin(1.5[1+(0.012RH)18]ho+250α3,1500α3)forfcm>35MPa$$\begin{array}{} \displaystyle \beta_{H}=\left\{\begin{array}{} \min(1.5[1+(0.012RH)^{18}]h_{o}+250,1500)\qquad~~\,\text{for}~~f_{cm}\le35\text{MPa}\\ \min(1.5[1+(0.012RH)^{18}]h_{o}+250\alpha_{3},1500\alpha_{3})~~\text{for}~~f_{cm}\gt35\text{MPa} \end{array}\right. \end{array}$$(8)

In these expressions, α₁, = (35/f_cm)^0.7, α₂ = (35/f_cm)^0.2, α₃ = (35/f_cm)^0.5; the relative humidity RH is expressed in percentage, the equivalent member size h_o is expressed in millimetre, an averaged compressive strength f_cm is expressed in megapascals and the loading time t_o is expressed in days. The implementation scheme, in which the creep curve is approximated by a set of nonageing Kelvin units, is partially based on the implicit algorithm proposed by Havlásek [14]. To eliminate the time dependency of Kelvin units, the amplitude A_μ (in the ageing material) of the current creep strain increment (see Eq. 9) is divided by a factor ν^cr, which amplifies, at early stages, or reduces, for old concrete, the current creep rate. Therefore, the increment of the current creep strains is expressed as follows:

Δεn+1cr=Do−11νn+1/2cr∑μ=1NAμ(1−βμ,n+1)σνμ,n$$\begin{array}{} \displaystyle \Delta\varepsilon_{n+1}^{cr}=\mathbf{D}_{o}^{-1}\frac{1}{\nu_{n+1/2}^{cr}}\sum_{\mu=1}^{N}A_{\mu}(1-\beta_{\mu,n+1})\sigma_{\nu\mu,n} \end{array}$$(9)

In the above expression, the visco-elastic projection matrix is denoted by Do−1$\begin{array}{} \displaystyle \mathbf{D}_{o}^{-1} \end{array}$ = E₂₈D^−l, νn+1/2cr$\begin{array}{} \displaystyle \nu_{n+1/2}^{cr} \end{array}$ , is an extra scaling factor amplifying the creep rate due to the ageing phenomenon, σ_νμ,n+1 represents the viscous effective stresses in the μ – th Kelvin unit (number of Kelvin units in the chain is denoted by N), the algorithmic factor β_μ,n+1 = exp(–Δt_n+1/τ_μ, the retardation time of μ -th Kelvin unit is denoted by τ_μ and the ultimate creep strain value for the μ -th Kelvin unit, obtained for unit stress, is denoted by A_μ.

The following updated procedure for the viscous stresses, in each Kelvin unit, is used:

σνμ,n+1=βμ,n+1σνμ,n+λμ,n+1Δσ¯n+1$$\begin{array}{} \displaystyle \boldsymbol{\sigma}_{\nu\mu,n+1}=\beta_{\mu,n+1}\boldsymbol{\sigma}_{\nu\mu,n}+\lambda_{\mu,n+1}\Delta\overline{\boldsymbol{\sigma}}_{n+1} \end{array}$$(10)

The algorithmic term λ_μ,n+1 is equivalent to λμ,n+1=(1−βμ,n+1)τμΔtn+1$\begin{array}{} \displaystyle \lambda_{\mu,n+1}=(1-\beta_{\mu,n+1})\frac{\tau_{\mu}}{\Delta t_{n+1}} \end{array}$ , while Δσ_n+1 represents the stress increment in effective (undamaged) configuration at any integration point in the structure. Assuming that the structure is loaded after 28 days, the algorithmic effective Young’s modulus is expressed as follows:

E¯=11E28+1νcr∑μ=1N(1−λμ,n+1)Aμ$$\begin{array}{} \displaystyle \overline{E}=\frac{1}{\frac{1}{E_{28}}+\frac{1}{\nu^{cr}}\sum\limits_{\mu=1}^{N}(1-\lambda_{\mu,n+1})A_{\mu}} \end{array}$$(11)

The ν^cr term appearing in Eq. 9 and Eq. 11 can be derived based on the following reasoning. Evolution of the creep strain in time, according to EC2, is expressed by the following equation:

εcr=A1(t−toβH+t−to)0.3β(to)$$\begin{array}{} \displaystyle \varepsilon^{cr}=A_{1}\bigg(\frac{t-t_{o}}{\beta_{H}+t-t_{o}}\bigg)^{0.3}\beta(t_{o}) \end{array}$$(12)

where A₁ϕ_RHβ(f_cm)/E₂₈,, while t – t_o is the time of creeping. Let us now introduce the reference creep strain curve corresponding to the loading time t_o = 2 8 days:

εrefcr=A1(t−toβH+t−to)0.3β(to=28)$$\begin{array}{} \displaystyle \varepsilon^{cr}_{ref}=A_{1}\bigg(\frac{t-t_{o}}{\beta_{H}+t-t_{o}}\bigg)^{0.3}\beta(t_{o}=28) \end{array}$$(13)

where t – t_o is also the time of creeping. This reference curve is approximated by the chain of Kelvin elements in which the retardation times τ_μ are adjusted with respect to the predicted time of analysis to be carried out. Once the retardation times τ_μ are set, the A_μ coefficients, in the chain of non-ageing Kelvin units, can be optimised using standard optimisation methods. Keeping the amplitudes of Kelvin elements, defined for the reference curve, constant while hiding the ageing effects in the ν^cr term is the main goal of this approach. To derive the ν^cr term, we assume the following compatibility condition for the creep strain rates at any creeping time t – t_o.

ε˙cr=1νcrε˙refcr$$\begin{array}{} \displaystyle \dot{\varepsilon}^{cr}=\frac{1}{\nu^{cr}}\dot{\varepsilon}_{ref}^{cr} \end{array}$$(14)

This yields the following expression for the ν^cr term:

νcr=β(to=28)β(to)$$\begin{array}{} \displaystyle \nu^{cr}=\frac{\beta(t_{o}=28)}{\beta(t_{o})} \end{array}$$(15)

where t_o is the age of concrete at the beginning of analysis. The proposed simplified approach yields an equivalent, relative to the EC2, visco-elastic implicit creep model.

In order to perform the preliminary wall dimensioning, two simulations (Cases A1 and A2), based on the NSO-LST computational strategy, were carried out. In Case A1, a nominal stiffness modulus E_cm =31000 MPa for the concrete class C25/30 was used, while in Case A2, a reduced (by 20%) stiffness modulus (E_cm = 25000 MPa) was used just to check the influence of this parameter.

The envelope of bending moments and the associated membrane forces in the wall, for both Cases A1 and A2, are shown in Fig. 5. As expected, larger moments are obtained for Case A1, but a reduction of the stiffness modulus by 20% yields <7% reduction in the maximum bending moment. A comparison of the envelopes of bending moments in the wall, for Case A1, based on all time instances registered until the last excavation step (labelled as A1^∗), and on the extended set of time instances registered until the foundation raft installation time, is shown in Fig. 6. One can notice that the maximum bending moment can be reduced at best by 24% if ground supports are used. The main role of these supports is to reduce wall movement until the time when the foundation raft is installed. This result is shown here just for a qualitative assessment of the factors influencing the distribution of bending moments in the wall and its deformations. Wall deformations for Cases A1, A1 ^∗, A2 and A2^∗ are shown in Fig. 7. As expected, larger deformations (by ≈ 7%) are obtained for Case A2, but exclusively at the time instance when the foundation raft is installed. At the time instance corresponding to the last excavation step, differences between the resulting deformations for the two cases are negligible. It is well visible that for the anchored diaphragm walls, protecting deep excavations, foundation rafts should be installed as soon as possible and ground supports should be used as well. Without ground supports, wall deformations increase >30%.

Figure 5

Figure 6

Figure 7

Based on the bending moment envelopes for Case A1 and the associated membrane forces, the following distribution of reinforcement is assumed (see Table 1) (A_s1 is the reinforcement placed at the internal wall face, while A_s2 is placed at the external one). This reinforcement will further be verified using non-linear subsoil–non-linear structure (NSO-NST) computational strategy. It has to be emphasised that due to the combined action of the bending moment and the membrane force, in the wall cross-section, but also the large eccentricity (M_xx/N_xx) case, the characteristic bending moments are amplified by a factor 1.35 (dead load partial factor), to design reinforcement, while the membrane forces are kept unchanged. This assumption leads to conservative dimensioning of reinforcement.

Table 1

Preliminary design of reinforcement in the wall based on the results achieved for Case A1.

General rules concerning the conduct of a fully non-linear analysis of subsoil and RC structures are not explicitly specified in EC2, nor in EC7. In the course of the designing process, the two limit states are always analysed, i.e. the ULS and the serviceability limit state (SLS). In the latter one, cracks in the RC structures and deformations, both in the structure and in subsoil, are checked. In the NSO-LST computational strategy, the characteristic stress resultants in the diaphragm wall are selectively increased (only bending moments) by the dead load partial safety factor, equal to 1.35. Then, standard dimensioning procedure is applied to set up the required reinforcement. It is worth mentioning that the soil parameters are not scaled by the partial safety factors to avoid non-physical situations (negative porosity, for instance). Therefore, the same strategy is used for RC structures in which the characteristic values of concrete strength (f_ck, f_ctk0,05 ) and stiffness E_cm are used to calibrate the modified Lee–Fenves model, while the steel reinforcement is characterised by the characteristic strength value f_yk. The SLS state, limited to the estimation of deformations, is relatively easy to analyse as all loadings are scaled with the partial safety factor equal to 1.0, the creep phenomenon is activated in the structure, while all parameters in the subsoil and concrete are taken as characteristic (this notion is rather problematic for the subsoil where HS model parameters are derived from a very limited set of experiments). However, assessment of the maximum crack opening is definitely more complicated.

In the computational FE model (FEM), the diaphragm wall is discretised using shell elements, but no extra interface between the concrete core and the steel reinforcement is introduced. Therefore, strains in concrete and steel are compatible, and such a model is unable to reproduce crack opening in a direct manner. To remedy this serious deficiency, one can assume that the difference between the averaged value of strain in steel and concrete (ε_sm – ε_cm) is approximately equal to the value of strain in the steel reinforcement, resulting from the conducted FEM analysis. If we accept this approximation, then standard EC2 procedure for computing maximum crack opening, w_k = S_r,max (ε_sm – ε_cm), can easily be adopted (see EC2 for more details).

The assessment of ULS state is definitely more problematic. Here, we assume that all parameters in the subsoil and structure are taken as characteristic, or let us say derived for the subsoil, while the computed stress resultants (bending moment and membrane force) at any point of the RC structure are projected on the domain bound by the bending moment–membrane force ultimate interaction curve (see Fig. 8) derived for concrete strength parameters f_cd and f_ctd, as well as the steel strength f_yd. Here, we assume that any stress resultant point (B) having coordinates {NxxB,MxxB⋅γ~}$\begin{array}{} \displaystyle \{N_{xx}^{B},M_{xx}^{B}\cdot\tilde{\gamma}\} \end{array}$ stays inside of the ultimate bending moment–membrane force interaction curve (see Fig. 8). The γ̃ value combines the two partial safety factors, i.e. the one corresponding to the dead load (1.35) and the one corresponding to the material one (1.4 for concrete [according to the Polish EC2 version] and 1.15 for steel). The upper bound is γ̃ = 1.35 1.4 ≈ 1.9, while the lower bound is γ̃ = 1.35 · 1.15 ≈ 1.55. Here, we will use the upper bound value γ̃ = 1.9.

Figure 8

As already mentioned, in the NSO-NST computational strategy, the modified Lee–Fenves CPDM is used, in which ageing and creep, compatible with the EC2 standard, is introduced [13, 5]. In the analysed case study, we assume that concrete age, at the beginning of the excavation stage, is approximately 28 days. The following set of material properties for concrete is used in this case study: E_cm = E₂₈ = 31000 MPa. ν = 0.2, γ = 25 kN/m³, f_c = 25 MPa, f_cbo /f_c = 0.4, f_cbo /f_c = 1.16, D̃_c = 0.435 at σ̃_c/f_c =1.0, G_c = 13.5^∗10^-3 MN/m, f_t = 1.8 MPa, D̃_t = 0.5 at σ̃_t/f_t = 0.5, G_t = 0.135^∗10^-3 MN/m, S_o = 0.2, α_p = 0.2 and α_d =1.0. Notion of all these parameters is explained in the original paper by Lee and Fenves [3] and in a previous paper by the author [5]. The above set of parameters yields sufficiently good match with the EC2 uniaxial stress–strain characteristics obtained in the uniaxial compression test. Additional parameters to characterise the creep (for RH =0.8) are as follows: ϕ_o β(f_cm)/E₂₈ = 1.14 ^∗ 10^–4 MPa ^–1, β_H = 2000 days and s = 0.38.

The last important aspect concerning modelling RC is related to the estimation of the so-called characteristic length l_RC, which should be used to scale the fracture energies in compression G_c and tension G_t in the FEM. For plain concrete, this length is equivalent to the FE size h^e, while in the case of RC structures, its value must carefully be estimated. Here, we assume that in the properly reinforced RC structures, total loss of concrete tensile strength appears at the strain at which steel reinforcement becomes plastic. This yields the following estimate for the l_RC value:

lRC=2GfEsfctk0,05fyk$$\begin{array}{} \displaystyle l_{RC}=\frac{2~G_{f}~E_{s}}{f_{ctk0,05~}f_{yk}} \end{array}$$(16)

The above expression is derived assuming a linear decrease in the concrete tensile strength with the axial tensile strain. In the considered case study, the characteristic length l_RC ≈ 0.06 m.

The two cases were analysed at first. In all of them, the NSO-NST computational strategy was used. Case B1 was dedicated to assess the ULS for preliminary design of reinforcement (creep is not activated), while Case B2 was designed to assess the SLS state (with creep).

The resulting bending moment envelopes and the associated membrane force profiles for Case B1 are shown in Fig. 9. In order to obtain the upper bound estimate for the bending moments in the case when ground supports are used and the foundation raft is installed as soon as possible, the two envelopes (B1 and B1 ∗) are shown in Fig. 10. The maximum bending moment in Case B1 is nearly 930 kNm/m. If we scale it by γ̃ = 1.9, then we will get nearly 1770 kNm/m. This value is larger than the one obtained from Case A1, for which the maximum registered characteristic moment was about 1180 kNm/m (moment used for dimensioning was 1.35∗1180 ≈ 1590 kNm/m). So, we see that non-linear computational strategy for γ̃ = 1.9 leads to a more conservative design.

Figure 9

Figure 10

When analysing results of the SLS state (see Fig. 11), we can see that in the best possible case (ground supports are used and foundation raft is installed fast), horizontal deformations of the wall are practically not influenced by creep, and small differences, in general, between the two computational strategies are observed. However, in the less-optimal construction technology, these differences become much bigger and creep effects become visible too. The maximum difference in terms of horizontal wall displacements between the Cases B2 and A1 is about 18%.

Figure 11

In order to check the maximum crack opening in the zone of the maximum bending moment, the maximum tensile strains have been traced in the steel reinforcement placed at the internal wall face at a depth of 14.5 m, for Case B2 (Fig. 12). The maximum registered value was εsmax$\begin{array}{} \displaystyle \varepsilon_{s}^{max} \end{array}$ = 6.2e – 4. For an assumed net steel cover c = 10 cm and equivalent steel bars with diameter ϕ_eq = 0.0225 (mixture of 20/25 mm bars), the h_c,ef = 0.275 m, ρ_p,eff = 0.0182 and s_r,max = 0.55 m (see the chapter on SLS state assessment in EC2). So, the estimated maximum crack opening was w_k = 0.55 · 6.2e – 4 · 1000 = 0.34 mm. This value should not be >0.3 mm; hence, the major longitudinal reinforcement was increased from 50 cm² to A_s1 = 57 cm². For this newly designed reinforcement, the maximum registered strain in steel was εsmax$\begin{array}{} \displaystyle \varepsilon_{s}^{max} \end{array}$ = 5.35e – 4, which yields the maximum crack opening of w_k = 0.028 mm.

Figure 12