% Matlab script twodd_half
% Two-dimensional half-space boundary element program
% with stress boundary conditions for Matlab 4.1
% This program is modified from Appendix B of the book
% "Boundary element methods in Solid Mechanics"
% by S.L. Crouch and A.M. Starfield, 1983
% Last revised on 2/8/94
clear %zero all variables
global PR1 PR2 CON CONS
global Sxxs Syys Sxys Sxxn Syyn Sxyn Uxs Uxn Uys Uyn
% Read elastic constants, remote stress field data, and boundary data from
% data files elastic.dat, remote.dat, and boundary.dat, respectively
load elastic.dat; % PR and E should be in one row;
PR = elastic(1); % PR = Poisson's ratio;
E = elastic(2); % E = Young's modulus;
clear elastic
load remote.dat % PXX, PYY, PXY should be in one row;
Pxx = remote(1); % Remote stress sigma xx;
Pyy = remote(2); % Remote stress sigma yy;
Pxy = remote(3); % Remote stress sigma xy;
clear remote
load boundary.dat
% Data in each row of this file pertains to one element
% Now separate out data types by column to form column vectors
XBEG = boundary(:,1); % element endpoint coordinate
YBEG = boundary(:,2); % element endpoint coordinate
XEND = boundary(:,3); % element endpoint coordinate
YEND = boundary(:,4); % element endpoint coordinate
BVs = boundary(:,5); % Shear traction boundary condition
BVn = boundary(:,6); % Normal traction boundary condition
NUM = length(BVs); % NUM = number of elements
clear boundary
% Define new elastic constants to be used internally
CON = 1/(4*pi*(1-PR));
CONS = E/(1+PR);
PR1 = 1-2*PR;
PR2 = 2*(1-PR);
% Define locations, lengths, orientations, and boundary conditions for boundary elements
% This information is stored here as column vectors of size (NUMx1)
x = (XBEG + XEND)/2; % element midpoints
y = (YBEG + YEND)/2; % element midpoints
XE = x;
YE = y;
XD = XEND-XBEG;
YD = YEND-YBEG;
A = sqrt(XD.*XD +YD.*YD)/2; % half-length of element
% In the lines below, B = beta = orientation of element relative to global x-axis
SINB = YD./(2*A); % sine beta
COSB = XD./(2*A); % cosine beta
SIN2B = 2*SINB.*COSB; % sine 2*beta
COS2B = COSB.*COSB-SINB.*SINB; % cosine 2*beta
SINB2 = SINB.*SINB; % sine squared of beta
COSB2 = COSB.*COSB; % cosine squared of beta
% Adjust stress boundary conditions to account for remote stresses.
% First resolve components of the remote stress onto the plane of each element...
SIGs = (Pyy-Pxx).*SIN2B/2 + Pxy.*COS2B;
SIGn = Pxx.*SINB2 - Pxy.*SIN2B + Pyy.*COSB2;
% and then subtract the resolved remote stress from the boundary stresses.
BVs = BVs - SIGs;
BVn = BVn - SIGn;
% The remote stress will be added back at the end of the solution.
% Compute influence coefficients between boundary elements.
% C(i,j) = effect at obs pt i due to a load at element j.
% The organization of the coefficients here is slightly different
% from Crouch & Starfield
% First dimension the vectors storing the influence coefficients
nobs=NUM*NUM;
set1;
for i=1:NUM
first = (i-1)*NUM+1;
last = i*NUM;
xe = XE(i);
ye = YE(i);
COSBe = COSB(i);
SINBe = SINB(i);
a = A(i);
Sd = 1;
Nd = 1;
coeffhs;
% Append sequentially produced vectors to master column vectors
append1;
end
% Calculate trig functions
% and direction cosine matrices for stress transformations
dircos1
% Reshape column vectors into square matrices, then clear column vectors
dimx = NUM;
dimy = NUM;
square1;
% Convert xy stresses returned from coeff to normal and shear tractions on elements
Cisjs = (Syys-Sxxs).*(SINN2B)/2 + Sxys.*COSS2B;
Cisjn = (Syyn-Sxxn).*(SINN2B)/2 + Sxyn.*COSS2B;
Cinjs = Sxxs.*SINNB2 - Sxys.*SINN2B + Syys.*COSSB2;
Cinjn = Sxxn.*SINNB2 - Sxyn.*SINN2B + Syyn.*COSSB2;
% Solve system of algebraic equations to get the displacement discontinuities
% First regroup the influence coefficient and boundary condition submatrices as shown:
% | [Cisjs] [Cisjn] | | [Ds] | = | [Bs] |
% | [Cinjs] [Cinjn] | | [Dn] | | [Bn] |
C = [Cisjs, Cisjn; Cinjs, Cinjn]; % C = influence coefficients;
% clear Cisjs Cisjn Cinjs Cinjn;
B = [BVs; BVn]; % B= Stress boundary conditions;
% Solve the matrix system [C][D]=[B] to get the displacement discontinuity vector D
D= C\B;
% Separate the D vector into subvectors Ds and Dn
Ds = D(1:NUM); % Ds elements are in upper half of D
Dn = D(NUM+1:2*NUM); % Dn elements are in lower half of D
% clear B C D
% Compute displacements and stresses on boundary elements
% The stresses obtained should match the boundary conditions
% First find the stresses and displacements, in an x-y reference frame, by superposition
UxNEG = Uxs*Ds + Uxn*Dn; % x-displacement on negative side of i;
UyNEG = Uys*Ds + Uyn*Dn; % y-displacement on negative side of i;
SIGxx = Pxx + Sxxs*Ds + Sxxn*Dn; % Note that remote stress is added back in;
SIGyy = Pyy + Syys*Ds + Syyn*Dn; % Note that remote stress is added back in;
SIGxy = Pxy + Sxys*Ds + Sxyn*Dn; % Note that remote stress is added back in;
% Now resolve DISPLACEMENTS in the xy frame into shear and normal components
% Start by solving for displacements on the NEGATIVE sides of the elements...;
UsNEG = UxNEG.*COSB + UyNEG.*SINB;
UnNEG = -UxNEG.*SINB + UyNEG.*COSB;
% and then add the displacement discontinuity at each element to get
% the displacement components on the positive sides of the elements
UsPOS = UsNEG - Ds;
UnPOS = UnNEG - Dn;
% The last step is to resolve the xy STRESSES to normal and shear stresses on the elements.
SIGs = (SIGyy-SIGxx).*(SIN2B/2) + SIGxy.*COS2B;
SIGn = SIGxx.*SINB2 - SIGxy.*SIN2B + SIGyy.*COSB2;
% Any writing to output files should occur now....
% profile % To plot profiles of displacements along boundaries
% Clear arrays in preparation for next section, saving vectors Ds and Dn
% clear xi yi xm ym
% clear Uxs Uxn Uys Uyn Sxxs Sxxn Syys Syyn Sxys Sxyn;
% clear SIGxx SIGyy SIGxy
% clear SIGs SIGn
% clear UxPOS UxNEG UyPOS UyNEG UsPOS UsNEG UnPOS UnNEG
% clear SINNB COSSB SINN2B COSS2B SINNB2 COSSB2
% COMPUTATION OF DISPLACEMENTS AND STRESSES AT SPECIFIED POINTS I IN BODY
check = 0;
if check ==1
% Begin by defining the coordinates of the gridpoints
twoddgrid
% Calculate influence coefficients at gridpoints
x = X(:); % Grid nodes in column vector form
y = Y(:); % Grid nodes in column vector form
% First dimension the vectors storing the influence coefficients
nobs=NUM2*NUM;
set1;
% Now loop through grid element by element
% The end result will be ten column vectors each having NUM2*NUM rows
for i=1:NUM
first = (i-1)*NUM2+1;
last = i*NUM2;
xe = XE(i);
ye = YE(i);
COSBe = COSB(i);
SINBe = SINB(i);
a = A(i);
Sd = Ds(i);
Nd = Dn(i);
coeffhs;
% Append sequentially produced vectors to master column vectors
append1;
end
% Reshape column vectors into square matrices the size of the
% observation grid, and in the process sum the contributions
% from each element
square2;
% Add back the remote field
Sxx = Pxx + Sxx;
Syy = Pyy + Syy;
Sxy = Pxy + Sxy;
S1 = SIG1(Sxx,Syy,Sxy);
S2 = SIG2(Sxx,Syy,Sxy);
tau = taumax(Sxx,Syy,Sxy);
mean = ave(Sxx,Syy,Sxy);
theta = angp(Sxx,Syy,Sxy);
end