Examples • Example Applications
Wharf Subjected to Earthquake Load
Problem Statement
Note
View this project in FLAC3D:
Pick “ExampleApplications”
Pick “Wharf”
Select “wharf.prj” to load
The project’s main data file is shown at the end of this example.
This is an illustrative example of a wharf subjected to earthquake load. The soil profile with three layers (4.5, 6.5 and 11.0 m) is shown in Figure 1. The water table is 2 m below that slope ground. The wharf and its piles are shown in Figure 2 in threedimension. Sand is at the slope which may have later spreading during an earthquake. The slope later spreading will push the piles and may affect the safety of the wharf. On the contrary, the piles have effect of spinning which will reduce the slope lateral spreading. This is an example of soilstructure interaction.
The liquefiable part of the middle layer will be modeled by the P2PSand model (colored by green in Figure 2 and Figure 3) and other soils will be modeled by the MohrCoulomb model plus hysteretic damping.
Several analysis stages are required:
Setup initial stress and porepressure;
Install structures;
Assign P2PSand model;
Dynamic analysis.
Setup initial stress and porepressure
This model requires configures of fluid and dynamic (model configure
fluid dynamic
. At this stage, both fluid and dynamic are off. This model follows the standard boundary conditions: fixed degrees of freedom at the bottom and roller boundaries on the side surfaces. Take advantage of the command zone face skin
, the boundary conditions can be
zone face apply velocityx 0 range union group 'west' group 'east'
zone face apply velocityy 0 range union group 'north' group 'south'
zone face apply velocity (0,0,0) range group 'bottom'
At this stage, all soils are assumed MohrCoulomb. The soil matrix material properties are list in Table 1.
Soil 
Density 
Bulk 
Shear 
Friction 
Cohesion 
\(Mg/m^3\) 
kPa 
kPa 
deg. 
kPa 

#1 
1.7 
5.10e5 
2.35e5 
40 
4 
#2 
1.5 
1.36e5 
6.30e4 
35 
2 
#3 
1.4 
1.36e5 
6.30e4 
30 
2 
The fluid related material properties are water density \(\rho_f\) = 1e3 \(kg/m^3\), porosity = 0.3. Gridpoint above the water table is set to zero saturation while others are 1 by default. Note that the input density should be dry density since the fluid is configured.
The pore pressures are initialized as the static using a gradient. The static water pressure should be applied to the top surfaces below the water table. These can be realized by two simple commands:
zone gridpoint initialize porepressure 200 ...
gradient 0 0 10 range positionz 0 20
zone face apply stressnormal 200 gradient 0 0 10 ...
range group 'top' positionz 10.99 20
The model is cycled to equilibrium by the command:
model solve elastic convergence 1
Install Structures
This stage is to set up the wharf structures.
Before installation structures, the model should be reinitialized into zero displacements and velocities and be cleared all yield states because any previously calculated displacement/velocity and yield states are just to obtain an initial field of stress and porepressure for subsequent stages. These can be realized by three simple commands:
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
zone initialize state 0
The platform is modeled by shell elements and the foundations are modeled by pile elements. The thickness of the shell is 0.3 m. The diameter of the piles is 0.6 meter. The structure geometry and properties are summarized in Table 2. The friction between the pile and the soil is assumed to of 30 degrees.
Structure 
Density 
Young 
Poisson 
Thickness 
Diameter 
\(Mg/m^3\) 
kPa 
m 
m 

Shell 
2 
2e8 
0.2 
0.3 

Pile 
2 
2e7 
0.2 
0.6 
To include the endbearing effect of the piles, anther set of links are created at the pile ends because by default the pile links are frictional. If a node must have multiple links, the links can be created by a different ‘side’. Here ‘side’ is a generalized concept similar to ‘slot’. In this model ‘side=2’ is used for endbearing links (with a group name ‘End’) between the piles and soil. The endbearing links are assumed normyield type. Note that the link displacement is defined as target displacement minus node displacement, the (downward) driving is considered ‘tension’ of link and the (upward) pulling is considered ‘compression’ of link so the link’s yieldcompression is set zero and yieldtension is set a large value to be consistent to the effect of endbearing. The commands on endbearing links are:
struct link create side=2 target zone group 'End' range positionz 5.9 6.1
struct link attach x=normalyield y=free z=free range group 'End'
struct link attach rotx=free roty=free rotz=free range group 'End'
struct link property x area=[math.pi*d^4/64] stiffness=1e9 ...
yieldcompression=0 yieldtension=1e20 range group 'End'
To make sure the pile heads are rigidly linked to the shell, the following command is used
structure node join range group 'wharf'
Again, the model is cycled to equilibrium for the installation of the structures.
Assign P2PSand model
P2PSand model is reassigned to some zones in the middle soil layer (colored by green in Figure 2 and Figure 3). In this analysis for simplicity and demonstration, we just adopt the internally calibrated material parameters for the P2PSand model. The primary parameter relativedensityinitial is input as 0.50 and 0.75 to analyze the effect of sand densities on the dynamic performance of the slope and the wharf. The reference pressure is assigned as pressurereference = 100 since the model’s unit on stress/pressure is Pa.
Since the P2PSand model is a stressdependent model, the initial stress will be used to determine the initial moduli, a small FISH function is adopted here:
fish operator ini_stress(z, modelname)
if zone.model(z) == modelname then
local pp = zone.pp(z)
zone.prop(z,'stressxxinitial') = zone.stress.xx(z) + pp
zone.prop(z,'stressyyinitial') = zone.stress.yy(z) + pp
zone.prop(z,'stresszzinitial') = zone.stress.zz(z) + pp
zone.prop(z,'stressxyinitial') = zone.stress.xy(z)
zone.prop(z,'stressyzinitial') = zone.stress.yz(z)
zone.prop(z,'stressxzinitial') = zone.stress.xz(z)
endif
end
[ini_stress(::zone.list, 'p2psand')]
Note that herein a keyword ‘operator’ (instead of ‘fish’) is used for multithreading. The argument ‘z’ is a zone pointer and ‘modelname’ is a constitutive model name. The command [ini_stress(::zone.list, ‘p2psand’)] is to run the FISH function “ini_stress” in a multithreaded mode
After assigning the P2PSand model, the model needs another reinitialization of displacement/velocity and yield state to zero. This time, structure nodes needs the similar reinitialization as well:
structure node initialize displacement (0,0,0)
structure node initialize displacementrotational (0,0,0)
So far, we have set up all necessary constitutive models and structures with allowable stress, porepressure and structure forces. It is ready for dynamic analysis.
Dynamic Analysis
The dynamic configure needs to be turned on for dynamic analysis. For quicker calculation, the ‘multistep’ is turned on as well. Because the earthquake can be considered a quick loading, the response of the soil is approximately considered undrained process. The fluid modulus is assigned a nonzero value 2e5 Pa so that the deformation from the soil matrix can cause the change of porepressure. Note that a smaller water bulk modulus 2e5 Pa rather than the real water bulk modulus 2e6 Pa is assigned to compromise the effect the free surface of water table (water is not completely confined). There commands are:
model dynamic active on
zone dynamic multistep on
zone gridpoint initialize fluidmodulus 2.0e5
The base of the model is considered compliant, so the input motion at the base should be the upward motion, or (1/2) of the outcrop motion. The input outcrop motion (Nahanni) before scaling is shown in Figure 4
In practice, the adopted input motion is usually be scaled to abide by some design codes or guidelines. In this example, a scaling factor of 5 is used so that the maximum acceleration is around 0.75g.
At the base, the quiet boundaries is used to accommodate the compliant base. The input acceleration motion is integrated into a velocity motion and then into a shear stress history. The commands are
zone gridpoint free velocity
table 'acc' import "Nahanni.acc"
program call "..\..\Fish\Library\int.fis"
[integrate('acc','vel')]
; since g=10, scaling factor=5
[global mf = 1800*math.sqrt(2.35e8/1800)*10*5.0]
zone face apply stressdip [mf] table 'vel' range group 'bottom'
zone face apply stressstrike 0.0 range group 'bottom'
zone face apply stressnorm 0.0 range group 'bottom'
zone face apply quiet range group 'bottom'
For zones with the MohrCoulomb model, the hysteretic damping is used to include the \(G/G_{max}\) effect. For zones with the P2PSand model, the constitutive relation has intrinsic energy dissipation so that theoretically extra energy dissipation is unnecessary. However, in order to reduce the high frequency noises, a small Rayleigh damping is added to zones. The dynamic time step is rarely decreased because the Rayleigh damping is very small. For structures, things are different because Rayleigh damping will significantly reduce the time step. Instead, Maxwell damping with a target 5% damping ratio is used for structures in this model:
zone dynamic damp hysteretic hardin 0.2 range cmodel 'mohrcoulomb'
zone dynamic damp rayleigh 0.0015 1.0
structure dynamic damp maxwell 0.0385 0.5 0.0335 3.5 0.052 25.0 ;target 5%
Freefield boundary is then applied following the basic the guideline dynamic analysis:
zone dynamic freefield on
The results for sand \(D_r\) = 0.50 (50%) are plotted in Figure 5 to Figure 7. Figure 5 plots porepressure histories at selected points (35,25,13), (46,25,14) and (58,25,15). Apparent excess porepressure buildsup are observed. Figure Figure #wharfssdr50
plots the contour of the maximum shear strain. Figure 7 plots the Xdisplacement histories of platform, slope top and the base. It is interesting to observe that the lateral displacement of the platform is less than the lateral displacement of the slope, which is generally the effect of soilstructure interaction: lateral spreading of the soil pushes the pile foundation and the pile foundation somewhat reinforces the stability of the slope.
Similar results for sand \(D_r\) = 0.75 (75%) are plotted in Figure 8 to Figure 10. Excess porepressure buildsup are observed in Figure 8. Lateral displacements of platform and slope shown in Figure 10 are less than those in Figure 7 which is apparently due to the sand in Figure 10 is denser.
break
Data Files
wharfini.dat
model new
model configure fluid dynamic
model largestrain off
model fluid active off
model dynamic active off
model gravity 0 0 10
;
zone import 'grid.flac3d'
zone face skin
;
zone face apply velocityx 0 range union group 'west' group 'east'
zone face apply velocityy 0 range union group 'north' group 'south'
zone face apply velocity (0,0,0) range group 'bottom'
;
zone cmodel assign mohrcoulomb
zone prop density=1.8 bulk=5.10e5 shear=2.35e5 ...
cohesion=4 friction=40.0 range group 'wharf:soil 1'
zone prop density=1.6 bulk=1.36e5 shear=6.30e4 ...
cohesion=2 friction=35.0 range group 'wharf:soil 2'
zone prop density=1.5 bulk=1.36e5 shear=6.30e4 ...
cohesion=2 friction=30.0 range group 'wharf:soil 3'
;
zone fluid cmodel assign isotropic
zone initialize fluiddensity 1.0
zone fluid property porosity 0.3
zone gridpoint initialize saturation 0 range positionz 20 100
zone gridpoint initialize porepressure 200 ...
gradient 0 0 10 range positionz 0 20
zone face apply stressnormal 200 gradient 0 0 10 ...
range group 'top' positionz 10.99 20
;
model solve elastic convergence 1
;
zone gridpoint initialize velocity (0,0,0)
zone gridpoint initialize displacement (0,0,0)
zone initialize state 0
;
struct pile create byline 42 20.5 6 42 20.5 22.1 segments 8 id 1
struct pile create byline 45 20.5 6 45 20.5 22.1 segments 8 id 2
struct pile create byline 42 23.5 6 42 23.5 22.1 segments 8 id 3
struct pile create byline 45 23.5 6 45 23.5 22.1 segments 8 id 4
struct pile create byline 42 26.5 6 42 26.5 22.1 segments 8 id 5
struct pile create byline 45 26.5 6 45 26.5 22.1 segments 8 id 6
struct pile create byline 42 29.5 6 42 29.5 22.1 segments 8 id 7
struct pile create byline 45 29.5 6 45 29.5 22.1 segments 8 id 8
;
[global d = 0.6]
structure pile property density 2 young 2e7 poisson 0.2 ...
crosssectionalarea [math.pi*d^2/4] ...
torsionconstant [math.pi*d^4/32] ...
moiy [math.pi*d^4/64] moiz [math.pi*d^4/64] ...
couplingstiffshear 1e5 couplingcohesionshear 0 ...
couplingfrictionshear 30 couplingstiffnormal 1e5 ...
couplingcohesionnormal 0 couplingfrictionnormal 30 ...
perimeter [math.pi*d] couplinggapnormal off
; Include endbearing effect
struct link create side=2 target zone group 'End' range positionz 5.9 6.1
struct link attach x=normalyield y=free z=free range group 'End'
struct link attach rotx=free roty=free rotz=free range group 'End'
struct link property x area=[math.pi*d^4/64] stiffness=1e6 ...
yieldcompression=0 yieldtension=1e20 range group 'End'
;
zone create brick size 3 5 1 point 0 39 17.5 21.1 point 1 48 17.5 21.1 ...
point 2 39 32.5 21.1 point 3 39 17.5 22.1 group 'wharfdeck'
structure shell create byzoneface id 9 ...
range positionx 39 48 positiony 17.5 32.5 positionz 22.0 22.2
structure shell group 'wharf'
zone delete range group 'wharfdeck'
;
structure shell cmodel assign elastic
structure shell property density 2 young 2e8 poisson 0.2 ...
thick 0.3 range id 9
structure node join range group 'wharf'
;
model solve convergence 1
model save 'static'
wharfp2pdr75.dat
model restore 'static'
;
zone cmodel assign p2psand range group 'wharf:soil 2' ...
positiony 5 45 positionx 0 70
zone property relativedensityinitial 0.75 pressurereference 100 ...
frictioncritical 35.0 range cmodel 'p2psand'
;
fish operator ini_stress(z, modelname)
if zone.model(z) == modelname then
local pp = zone.pp(z)
zone.prop(z,'stressxxinitial') = zone.stress.xx(z) + pp
zone.prop(z,'stressyyinitial') = zone.stress.yy(z) + pp
zone.prop(z,'stresszzinitial') = zone.stress.zz(z) + pp
zone.prop(z,'stressxyinitial') = zone.stress.xy(z)
zone.prop(z,'stressyzinitial') = zone.stress.yz(z)
zone.prop(z,'stressxzinitial') = zone.stress.xz(z)
endif
end
[ini_stress(::zone.list, 'p2psand')]
;
model solve convergence 1
;
zone gridpoint initialize displacement (0,0,0)
zone gridpoint initialize velocity (0,0,0)
zone initialize state 0
structure node initialize displacement (0,0,0)
structure node initialize displacementrotational (0,0,0)
;
model dynamic active on
zone dynamic multistep on
zone gridpoint initialize fluidmodulus 2.0e5
;
zone gridpoint free velocity
table 'acc' import "Nahanni.acc"
[table.as.list('vel') = table.integrate('acc')]
; since g=10, scaling factor=5
[global mf = 1.8*math.sqrt(2.35e5/1.8)*10*5.0]
zone face apply stressdip [mf] table 'vel' range group 'bottom'
zone face apply stressstrike 0.0 range group 'bottom'
zone face apply stressnorm 0.0 range group 'bottom'
zone face apply quiet range group 'bottom'
;
zone dynamic damp hysteretic hardin 0.2 range cmodel 'mohrcoulomb'
zone dynamic damp rayleigh 0.0015 1.0
; target 5% for structure elements (& transferred to deformabale links)
structure dynamic damp maxwell 0.0385 0.5 0.0335 3.5 0.052 25.0
;
zone dynamic freefield on
;
history delete
history interval 10
model history name 'time' dynamic timetotal
program call 'excpphist.f3fis'
zone history name 'pp1' porepressure source gridpoint position 35 25 13
zone history name 'pp2' porepressure source gridpoint position 46 25 14
zone history name 'pp3' porepressure source gridpoint position 58 25 15
struct node history name 'dispx1' displacementx componentid 86
zone history name 'dispx2' displacementx position 50 24 0
zone history name 'dispx3' displacementx position 48 24 22
;
model dynamic timestep maximum 1.00e4
model solve timetotal 21
;
history export 'dispx1' vs 'time' table 'disp1'
table 'disp1' export 'disp1dr75' truncate
;
model save 'wharf_dr75'
⇐ Undrained Cylindrical Cavity Expansion in a CamClay Medium  Simulation of PullTests for Fully Bonded Rock Reinforcement ⇒
Was this helpful? ...  Itasca Software © 2022, Itasca  Updated: Mar 09, 2023 