Computational Mechanics and Design Group

Department of Civil & Structural Engineering

MATLAB script for DLO analysis of geotechnical problems

Description

Discontinuity Layout Optimization (DLO) is a recently developed numerical limit analysis procedure, which can be programmed relatively easily. Here the demonstration MATLAB script for geotechnical problems using DLO presented at NUMGE conference, Norway, 2010 is provided.

Sample output

Example 1 - undrained retaining wall

dlo(9, 6, 0, 0, 2, 5, 1, 0, 0)

Solution 1.2% above exact value of (2 + pi) / 2

Example 2 - Prandtl problem

dlo(13, 7, 0, 0, [5; 2; 9], 1, 1, 0, 0)

Solution 1.2% above exact value of (2 + pi)

Example 3 - drained anchor

dlo(4, 4, [5; 0; 3], 0, 2, 1, 0, 20, 1)

Solution 0.6% above exact value of 6.912

Example 4 - surcharge near a vertical cut

dlo(3, 2, 0, 2, [4; 2; 1], 0, 1, 0, 0) 

Code

%===================================================================================%
%       Discontinuity Layout Optimization (DLO) Demonstration MATLAB Script         %
% (C) 2009 Computational Mechanics & Design Research Group, University of Sheffield %
%    v2.0 (Dec 2009) - added self-weight and up to 2 boundary conditions per edge   %
%===================================================================================%
%  edge nos: o C o  types: 0 = rigid,        1 = symmetry plane                     %
%            D   B         2 = free          3 = 'flexible' load                    %
%            o A o         4 = 'rigid' load  5 = 'rigid' load constrained laterally %
%                                                                                   %
% Usage: dlo(xmax, ymax, edgeA, edgeB, edgeC, edgeD, cohesion, phiDegs, unitWeight) %
% N.B. if two conditions along edgeX, replace with [edgeX1; edgeX2; lengthOfedgeX2] %
%                                                                                   %
%  Example 1: dlo(9, 6, 0, 0, 2, 5, 1, 0, 0) (Retaining wall; exact = (2 + pi) / 2) %
%                                                                                   %
%  Example 2: dlo(13, 7, 0, 0, [5; 2; 9], 1, 1, 0, 0) (Prandtl; exact = 2 + pi)     %
%                                                                                   %
%  Example 3: dlo(4, 4, [5; 0; 3], 0, 2, 1, 0, 20, 1) (Anchor; exact = 6.912)       %
%                                                                                   %
function dlo(xmax, ymax, edgeA, edgeB, edgeC, edgeD, cohesion, phiDegs, unitWeight)
        nodes = createNodes(xmax, ymax);
        discs = createDiscontinuities(nodes, xmax, ymax, edgeA, edgeB, edgeC, edgeD);
            B = compatibilityMatrix(nodes, discs);
[objP padN N] = plasticMultiplierTerms(nodes, discs, cohesion, phiDegs);
           fD = selfWeight(nodes, discs, ymax, unitWeight);
           fL = unitLoad(discs);
         tied = tieDiscs(discs);
 [vars, soln] = solve(B, N, padN, fD, fL, tied, objP);
                plotMechanism(nodes, discs, vars, xmax, ymax);
                soln
%===================================================================================%
% SETUP NODES %
function nodes = createNodes(nx, ny)
nodes = [];
for y = 0 : ny
    for x = 0 : nx
        nodes = [nodes; size(nodes, 1) + 1 x y];
    end
end
% SETUP DISCONTINUITY ELEMENTS %
function discs = createDiscontinuities(nodes, xmax, ymax, edgeA, edgeB, edgeC, edgeD)
discs = []; 
edgeA = sparse(3, 1) + edgeA; edgeB = sparse(3, 1) + edgeB; edgeC = sparse(3, 1) + edgeC; edgeD = sparse(3, 1) + edgeD;
for i = 1 : size(nodes, 1) - 1
    for j = i + 1 : size(nodes, 1)              
        edgeType = edgeA(1) * ((nodes(i, 3) + nodes(j, 3) == 0)        && (nodes(i, 2)  <=  (xmax - edgeA(3))) && (nodes(j, 2)  <=  (xmax - edgeA(3))))...
                 + edgeA(2) * ((nodes(i, 3) + nodes(j, 3) == 0)        && (nodes(i, 2)  >=  (xmax - edgeA(3))) && (nodes(j, 2)  >=  (xmax - edgeA(3))))...
                 + edgeB(1) * (nodes(i, 2) + nodes(j, 2) == 2 * xmax) ...
                 + edgeC(1) * ((nodes(i, 3) + nodes(j, 3) == 2 * ymax) && (nodes(i, 2)  <=  (xmax - edgeC(3))) && (nodes(j, 2)  <=  (xmax - edgeC(3))))...
                 + edgeC(2) * ((nodes(i, 3) + nodes(j, 3) == 2 * ymax) && (nodes(i, 2)  >=  (xmax - edgeC(3))) && (nodes(j, 2)  >=  (xmax - edgeC(3))))...
                 + edgeD(1) * (nodes(i, 2) + nodes(j, 2) == 0);
        if (gcd(abs(nodes(j, 2) - nodes(i, 2)), abs(nodes(j, 3) - nodes(i, 3))) < 1.0001)
            discs = [discs; i j edgeType sqrt((nodes(i, 2) - nodes(j, 2))^2 + (nodes(i, 3) - nodes(j, 3))^2)];
        end
    end
end
% SETUP COMPATIBILIY MATRIX %
function B = compatibilityMatrix(nodes, discs);
[nCons nVars] = deal(2 * size(nodes, 1), 2 * size(discs, 1));
B = sparse(nCons, nVars);
for i = 1 : size(discs, 1)
    [n1 n2 len] = deal(discs(i, 1), discs(i, 2), discs(i, 4));
    [conN1 conN2 var] = deal(2 * nodes(n1, 1) - 1, 2 * nodes(n2, 1) - 1, 2 * i - 1);
    [cosine sine] = deal((nodes(n1, 2) - nodes(n2, 2)) / len, (nodes(n1, 3) - nodes(n2, 3)) / len);
    B_local = [cosine -sine; sine cosine];
    B([conN1 conN1 + 1], [var var + 1]) =  B_local;
    B([conN2 conN2 + 1], [var var + 1]) = -B_local;
end
% SETUP PLASTIC MULTIPLIER MATRIX %
function [objP, padN, N] = plasticMultiplierTerms(nodes, discs, c, phiDegrees);
objP = sparse(0, 0); padN = sparse(0, 0); N = sparse(0, 0); count = 1; tan_phi = tan(pi * phiDegrees / 180);
for i = 1 : size(discs, 1)
    if (discs(i, 3) < 2) % i.e. if rigid or symmetry
        [eff_c eff_tanPhi] = deal(c * (discs(i, 3) ~= 1), tan_phi * (discs(i, 3) ~= 1));
        [con var] = deal(2 * count - 1, 2 * count - 1);
        N_local = [1 -1; eff_tanPhi eff_tanPhi];
        N([con con + 1], [var var + 1]) = N_local;
        padN_local = sparse(2, 2 * size(discs, 1));
        padN_local(1, 2 * i - 1) = -1;
        padN_local(2, 2 * i) = -1;
        padN = [padN; padN_local];
        objP = [objP eff_c * discs(i, 4) eff_c * discs(i, 4)];
        count = count + 1;
    end
end
% SELF WEIGHT LOADS %
function fD = selfWeight(nodes, discs, ymax, unitWeight);
fD = [];
for i = 1 : size(discs, 1)
    [n1 n2 len] = deal(discs(i, 1), discs(i, 2), discs(i, 4));
    [cosine sine] = deal((nodes(n1, 2) - nodes(n2, 2)) / len, (nodes(n1, 3) - nodes(n2, 3)) / len);
    stripWidth = nodes(n2, 2) - nodes(n1, 2);
    stripWeight = unitWeight * stripWidth * (ymax - 0.5 * (nodes(n1, 3) + nodes(n2, 3)));
    fD = [fD; -sine * stripWeight; -cosine * stripWeight];
end
% EXTERNAL LIVE LOADS %
function fL = unitLoad(discs);
fL = sparse(1, 2 * size(discs, 1));
for i = 1 : size(discs, 1)
    if (discs(i, 3) >= 3)
        fL(1, 2 * i) = 1;  %normal stress
    end
end
% LINK ELEMENT DISPLACEMENTS (ONLY IF RIGID LOADS) %
function tied = tieDiscs(discs); % for simplicity implement by adding constraints
tied = sparse(0, 2 * size(discs, 1)); ilast = 0;
for i = 1 : size(discs, 1)
    if (discs(i, 3) >= 4) % i.e. if rigid load
        if ((ilast ~= 0) || (discs(i, 3) == 5))
            extra = sparse(2, 2 * size(discs, 1));
            extra(1, 2 * i - 1) = 1;
            if ilast ~= 0
               extra(1, 2 * ilast - 1) = -1;
               extra(2, 2 * i) = 1; extra(2, 2 * ilast) = -1;
            end
            tied = [tied; extra];
        end
        ilast = i;
    end
end
% FORM AND SOLVE LINEAR PROGRAMMING PROBLEM %
function [variables, solution] = solve(B, N, padN, fD, fL, tied, objP);
bigNumber = 1000; warning off MATLAB:divideByZero;
[padB padFL padTied] = deal(sparse(size(B, 1), size(N, 2)), sparse(1, size(N, 2)), sparse(size(tied, 1), size(N, 2)));
equalityMatrix = [B padB; padN  N; fL padFL; tied padTied];
equalityRHS = sparse(size(equalityMatrix, 1), 1);
equalityRHS(size(B, 1) + size(N, 1) + 1, 1) = 1;
obj = [fD; objP']; % objP' contains plastic multiplier terms
lowerB = [-bigNumber * ones(size(B, 2), 1); sparse(size(N, 2), 1)];
upperB = [bigNumber * ones(size(B, 2), 1); bigNumber * ones(size(N, 2), 1)];
[variables, solution] = linprog(obj, [], [], equalityMatrix, equalityRHS, lowerB, upperB);
% PLOT MECHANISM %
function plotMechanism(nodes, discs, variables, xmax, ymax);
clf; hold on; axis equal; 
set(gca,'XLim',[0 xmax], 'XTick', [0:1:xmax], 'YLim',[0 ymax], 'YTick', [0:1:ymax], 'FontSize', 12)
xlabel('x', 'FontSize', 14); ylabel('y', 'FontSize', 12);
for pass = 1 : 2 % plot active elements on top of inactive elements
    for i = 1 : size(discs, 1)
       [xp yp] = deal([nodes(discs(i, 1), 2); nodes(discs(i, 2), 2)], [nodes(discs(i, 1), 3); nodes(discs(i, 2), 3)]);
       if ((pass == 2) && ((abs(variables(2 * i)) > 0.001) || (abs(variables((2 * i) - 1)) > 0.001)) && (discs(i, 3) ~= 1))
           plot(xp, yp, 'r-', 'LineWidth', 3);        % active = red
       elseif (pass == 1)
          plot(xp, yp, '-', 'Color', [0.9 0.9 0.9]); % inactive = light grey
       end
    end
end

Publication(s)

Gilbert M, Smith CC, Haslam IW, Pritchard TJ (2010). Application of discontinuity layout optimization to geotechnical limit analysis problems. In Numerical Methods in Geotechnical Engineering - Proceedings of the 7th European Conference on Numerical Methods in Geotechnical Engineering , Abstract: Limit analysis provides a long-established and powerful means of assessing the stability of...