Layered Damped Soil on Elastic Rock - Site Amplification in Kramer Book

Hi,

I would like to calculate the amplification factor in multilayer damped soil on elastic rock, which is also included in Kramer's book geotechnical earthquake engineering (page 268). After making the inputs of the rock and the soil as contained in the code, I need to extract the amplification that I expected to see on the surface, but I can't find where I made a mistake. While getting the correct result in another licensed software, matlab gives very high values, especially the first two waves. Is there anyone who is working on this issue who can help me ? Thank you

Kramer Book

Code:

% Define soil layer properties
H = [20, 30, 40];          % Thickness of the layers (m)
Vs = [200, 300, 400];      % Shear wave velocities of the layers (m/s)
rho = [1800, 1900, 2000];  % Densities of the layers (kg/m^3)
xi = [0.05, 0.04, 0.03];   % Damping ratios of the layers
% Rock properties
Vs_rock = 800;             % Shear wave velocity of the rock (m/s)
rho_rock = 2200;           % Density of the rock (kg/m^3)
xi_rock = 0.02;            % Damping ratio of the rock
% Define frequency range
f = linspace(0.1, 20, 1000);  % Frequency (Hz)
omega = 2 * pi * f;          % Angular frequency (rad/s)
% Create an empty array to store amplification values
A = zeros(size(f));
% Calculate amplification for each frequency
for j = 1:length(f)
    % Complex wave number for rock
    k_rock = (omega(j) / Vs_rock) * (1 - 1i * xi_rock);
    
    % Reflection at the rock base (initial transfer matrix is identity)
    T_total = eye(2);  % Identity matrix as the starting point
    
    % Calculate wave number and transfer matrix for each soil layer
    for i = 1:length(H)
        % Complex wave number for each soil layer (using updated formula)
        k = (omega(j) / Vs(i)) * (1 - 1i * xi(i));
        
        % Transfer matrix for the current layer
        Ti = [cos(k * H(i)), sin(k * H(i)) / (rho(i) * Vs(i)); 
            -rho(i) * Vs(i) * sin(k * H(i)), cos(k * H(i))];
        
        % Update the total transfer matrix
        T_total = Ti * T_total;
    end
    
    % Apply free surface boundary condition at the surface
    A(j) = 1 / abs(T_total(1,1));
end
% Plot the frequency-amplification graph
figure;
plot(f, A, 'LineWidth', 2);
xlabel('Frequency (Hz)');
ylabel('Amplification Factor');
title('Layered Soil Amplification over Elastic Rock (with Damping)');
grid on;

Kramer Book Screenshots:

[William Rose]

@ibrahim can beziklioglu, the book, pp. 268-269, does not indicate how k* and complex shear modulus G* are related to shear velocity vs, density rho, and damping xi.

Your code includes

k = (omega(j) / Vs(i)) * (1 - 1i * xi(i));

in other words,

Can you confirm that that is correct?

I do not see G or G* in your code. I am asking becausue I want to verify that your code implements equations 7.34a,b correctly (i.e. as the book says). Therefore I want to reproduce equation 7.35. Equation 7.35 says is complex, and it uses . How is (which is complex, I assume) related to (which is real, and is given in the table of constants)?

[ibrahim can]

Hi,

Actually this is,

[William Rose]

@ibrahim can beziklioglu,

Thank you for the additional information.

The new page of equations which you sent includes the following equation, which is not correct:

The approximation

is reasonable if . This same assumption about xi is also necessary to justify the approximation

which also appears on the page you posted above. You used xi=0.05, 0.04, 0.03, in your code, so the approximations are justified.

I will re-look at your code and equations. It will take me some time.

[ibrahim can]

I set up the equations again with the formulas specified in the book, this time I wanted to make two layers on the rock and see if the equations worked, but I got a similar error. The transfer function comes out as multiple and I can't find my error.

% MATLAB code to calculate F_{13} based on the given parameters

% Parameters

h1 = 15; % Thickness of layer 1 (m)

h2 = 30; % Thickness of layer 2 (m)

h3 = Inf; % Thickness of layer 3 (rock, infinite thickness)

Vs1 = 150; % Shear wave velocity in layer 1 (m/s)

Vs2 = 350; % Shear wave velocity in layer 2 (m/s)

Vs3 = 1050; % Shear wave velocity in layer 3 (m/s)

xi1 = 0.05; % Damping ratio in layer 1

xi2 = 0.03; % Damping ratio in layer 2

xi3 = 0.02; % Damping ratio in layer 3

rho1 = 2000; % Density of layer 1 (kg/m^3)

rho2 = 2100; % Density of layer 2 (kg/m^3)

rho3 = 2500; % Density of layer 3 (kg/m^3)

% Frequency range

f = linspace(0.01, 25, 500); % Frequency range from 0.01 to 25 Hz (avoid f=0)

omega = 2 * pi * f; % Angular frequency (rad/s)

% Complex shear wave velocities

Vs1_star = Vs1 * (1 + 1i * xi1);

Vs2_star = Vs2 * (1 + 1i * xi2);

Vs3_star = Vs3 * (1 + 1i * xi3);

% Wave numbers (complex)

k1_star = omega ./ Vs1_star;

k2_star = omega ./ Vs2_star;

k3_star = omega ./ Vs3_star;

% Impedance ratios (complex)

alpha1_star = (rho1 * Vs1_star) ./ (rho2 * Vs2_star);

alpha2_star = (rho2 * Vs2_star) ./ (rho3 * Vs3_star);

% Initialize A1 and B1

A1 = B1; % Shear stress must be zero

% Calculate A2 and B2

A2 = 0.5 * (A1 * (1 + alpha1_star) .* exp(1i * k1_star * h1) + B1 * (1 - alpha1_star) .* exp(-1i * k1_star * h1));

B2 = 0.5 * (A1 * (1 - alpha1_star) .* exp(1i * k1_star * h1) + B1 * (1 + alpha1_star) .* exp(-1i * k1_star * h1));

% Calculate A3 and B3

A3 = 0.5 * (A2 .* (1 + alpha2_star) .* exp(1i * k2_star * h2) + B2 .* (1 - alpha2_star) .* exp(-1i * k2_star * h2));

B3 = 0.5 * (A2 .* (1 - alpha2_star) .* exp(1i * k2_star * h2) + B2 .* (1 + alpha2_star) .* exp(-1i * k2_star * h2));

% Normalize A3 and B3 with respect to A1

a1 = A1 / A1; % Normalized to 1

b1 = B1 / B1; % Normalized to 1

a2 = A2 / A1; % Normalize to A1

b2 = B2 / B1; % Normalize to B1

a3 = A3 / A1; % Normalize to A1

b3 = B3 / B1; % Normalize to B1

% Calculate F13

F13 = abs((a1 + b1) ./ (a3 + b3)); % Take the modulus of the ratio

% Plot the result

figure;

plot(f, F13, 'LineWidth', 1.5);

grid on;

xlabel('Frequency (Hz)');

ylabel('|F_{13}|');

title('Frequency Response Function |F_{13}|');

legend('|F_{13}|');

[ibrahim can]

Thank you for your help. I solved it.

Answer 1

[William Rose]

@ibrahim can beziklioglu,

You included two plots: a matlab-generated plot of Amplification factor as a funciton of frequency, and a plot of transfer function H versus frequency, from another source. I think the transfer function plot (red trace) is from top of layer 2 to top of layer 3, due to the red rectangle around layer 2 at the left of the plot. Do you agree? If so, then your Matlab plot should also be for the amplificationm from top of layer 2 to top of layer 3.

[ibrahim can]

Thank you for your help. I solved it.