% x: a 2*1 vector denoting the position of field point;
% xi: a 2*1 vector denoting the position of point force (source point);
% xc: cos_theta of point x on kth element;
% xs: sin_theta of point x on kth element;
global A B mu
z=x(1)+x(2)*mu; z_hat=xi(1)+xi(2)*mu;
dz=z-z_hat;
switch Ytype
case 0 % calculate [U_star,T_star];
logz=log(dz); f=(xc+xs*mu)./dz;
U0=imag(A*diag(logz)*A.')/pi; T0=imag(A*diag(f)*B.')/pi; % U0=U_star; T0=T_star;
case 1 % calculate [U1_star;U2_star,T1_star;T2_star];
zi=1./dz; fp=(xc+xs*mu)./(dz.^2);
U1_star=-imag(A*diag(zi)*A.')/pi; T1_star=imag(A*diag(fp)*B.')/pi;
U2_star=-imag(A*diag(mu.*zi)*A.')/pi; T2_star=imag(A*diag(mu.*fp)*B.')/pi;
U0=[U1_star;U2_star]; T0=[T1_star;T2_star];
end
15.2.2 A Half Plane – G2half2D (略) 15.2.3 Interfaces – G3interface2D (略) 15.2.4 Holes – G4hole2D (略)
15.2.5 Cracks (略)
15.2.6 Rigid Inclusions – G6Rinclusion2D (略) 15.2.7 Elastic Inclusions – G7Einclusion2D (略)
15.3 Fundamental Solutions for Coupled Stretching-Bending Analysis (略) 15.3.1 An Infinite Laminate – G1infCouple (略)
Function G1infCouple:
function [U0,T0]=G1infCouple(x,xi,xc,xs,jn,Ytype)
% U0,T0,Ytype: see the explanation stated in the function: GreenCouple;
% x: a 2*1 vector denoting the position of field point;
% xi: a 2*1 vector denoting the position of point force (source point);
% xc: cos_theta of point x on kth element, xc=nc=nt2-nt1 for Ytype=2 or 3 (corners);
% xs: sin_theta of point x on kth element, xs=sc=st2-st1 for Ytype=2 or 3 (corners);
% jn: jth node point; jn=0 stands for internal point;
global A B mu % employed;
global F Fd1 Fd2 Fint Fint_1 f3 f3d f3d2 f3int % generated;
z=x(1)+x(2)*mu; z_hat=xi(1)+xi(2)*mu;
dz=z-z_hat; dzi=1./dz; logz=log(dz);
logz_int=dz.*(logz-1); logz2_int=(dz.^2).*(2*logz-3)/4;
% calculation of function matrices F,f3,...,etc, eqn.(15.38b),(15.39),(15.41b,c,d);
switch Ytype
case 0 % calculate F,Fint,Fd1,Fd2,f3,f3int,f3d,f3d2, eqn.(15.38b),(15.39);
15‐12
F=diag(logz)*A.'/(2*pi*i); Fint=diag(logz_int)*A.'/(2*pi*i);
Fd1=diag(dzi)*A.'/(2*pi*i); Fd2=-diag(dzi.^2)*A.'/(2*pi*i);
f3=Fint(:,3); f3int=diag(logz2_int)*A(3,:).'/(2*pi*i); f3d=F(:,3); f3d2=Fd1(:,3);
case 1 % calculate F_1,Fint_1,Fd1_1,Fd2_1,f3_1,f3int_1,f3d_1,f3d2_1, eqn.(15.41b,c,d);
% F:F_1, Fint:Fint_1, Fd1:Fd1_1, Fd2:Fd2_1, f3:f3_1, f3int:f3int_1, f3d:f3d_1, f3d2:f3d2_1;
F=-diag(dzi)*A.'/(2*pi*i); Fint=-diag(logz)*A.'/(2*pi*i);
Fd1=diag(dzi.^2)*A.'/(2*pi*i); Fd2=-diag(dzi.^3)*A.'/(pi*i);
f3=Fint(:,3); f3int=-diag(logz_int)*A(3,:).'/(2*pi*i);
f3d=F(:,3); f3d2=Fd1(:,3);
case 2 % calculate Fint,Fd1,f3int,f3d, eqn.(15.38b),(15.39);
Fint=diag(logz_int)*A.'/(2*pi*i); Fd1=diag(dzi)*A.'/(2*pi*i);
f3int=diag(logz2_int)*A(3,:).'/(2*pi*i); f3d=diag(logz)*A(3,:).'/(2*pi*i);
case 3 % calculate Fint_1,Fd1_1f3int_1,f3d_1, eqn.(15.41b,c,d);
% Fint:Fint_1, Fd1:Fd1_1, f3int:f3int_1, f3d:f3d_1;
Fint=-diag(logz)*A.'/(2*pi*i); Fd1=diag(dzi.^2)*A.'/(2*pi*i);
f3int=-diag(logz_int)*A(3,:).'/(2*pi*i); f3d=-diag(dzi)*A(3,:).'/(2*pi*i);
case 10 % calculate F,Fint,f3,f3int, eqn.(15.38b),(15.39);
F=diag(logz)*A.'/(2*pi*i); Fint=diag(logz_int)*A.'/(2*pi*i);
f3=Fint(:,3); f3int=diag(logz2_int)*A(3,:).'/(2*pi*i);
case 11 % calculate F_1,Fint_1,f3_1,f3int_1, eqn.(15.41b,c,d), F:F_1, Fint:Fint_1, f3:f3_1, f3int:f3int_1;
F=-diag(dzi)*A.'/(2*pi*i); Fint=-diag(logz)*A.'/(2*pi*i);
f3=Fint_1(:,3); f3int=-diag(logz_int)*A(3,:).'/(2*pi*i);
end
[U0,T0]=GreenCouple(A,B,mu,xc,xs,jn,Ytype);
15.3.2 Holes – G4holeCouple (略) 15.3.3 Cracks (略)
15.3.4 Inclusions – G7inclusionCouple (略)
15.4 Two-Dimensional Anisotropic Elastic Analysis – Basic Version (略) 15.4.1 Mesh Generation of Boundary Element – BEMmesh (略)
Function BEMmesh is designed to organize the mesh information through three input textfiles: input_xn, input_node12 and input_bc for both 2D problems and coupled stretching-bending problems. The first one is for the coordinates of boundary nodes, the second one is for the element-node connection, and the third one is for the boundary conditions. Additional input files input_corner12.txt is required for coupled stretching-bending problems, which is the information about the corner-node connection. The computer code of this function is as follows.
function BEMmesh(Ltype1)
% BEMmesh generate boundary element mesh for 2D_BEM and couple_BEM through the input of three
% text files "input_xn.txt", "input_node12.txt" and "input_bc.txt" for both 2D_BEM & couple_BEM.
% Additional input files "input_corner12.txt" is required for coupled_BEM.
%
% Followings are the global variables employed in this function:
% Btype: Btyep=1: linear element; Btype=2: quadratic element; Btype=3:linear element with cubic deflection;
% dim: dimension of material eigenvector matrix A;
% xn(knode,i): knode: kth-node; i=1,x1-coordinate; i=2,x2-coordinate, input through 'input_xn.txt';
% ... coordinate system should be constructed based upon the coordinate of the Green's function;
% node12(kelement,i): kelement: kth-element;
% i=1,node 1 of kth-element; i=2,node 2 of kth-element, input through 'input_node12.txt';
% i=3,node 3 of kth-element for quadratic element;
% corner12(k,i): the node behind of the kth corner; corner12(k,2): the node ahead of the kth corner;
%
% bc(knode,i): knode: kth-node;
% plane problems (Dtype=1 or 2 or ...):
% i=1~3: displacement-prescribed or traction-prescribed, 1:displacement-prescribed, 0:traction-prescribed;
% i=4~6: prescribed values of displacements or tractions, input through 'input_bc.txt';
% i=7~9: prescribed values of acceleration for dynamic problems (if displacement is prescribed);
% i=7~9: initial values of displacement for dynamic problems (if displacement is unknown in time history, traction-prescribed);
% i=10~12: initial values of velocity for dynamic problems (if displacement is unknown in time history, traction-prescribed);
% coupled stretching-bending problems (Dtype=3 or ...):
% i=1~4: displacement-prescribed or traction-prescribed, 1:displacement-prescribed, 0:traction-prescribed;
% i=5~8: prescribed values of displacements or tractions, input through 'input_bc.txt';
% ...if meshes conatin corner points, additional variables - corner force should be considered,
%...which is added after nnode. Thus, when knode>nnode
% ...i=1: deflection-prescribed or corner force-prescribed, 1:deflection-prescribed, 0:corner force-prescribed;
%
% Followings are the global variables generated in this function:
% xm1(i,k): position of the 1st node of kth element; k: kth-element; i=1,x1-coordinate; i=2,x2-coordinate;
% xm2(i,k): position of the 2nd node of kth element; k: kth-element; i=1,x1-coordinate; i=2,x2-coordinate;
% xm3(i,k): position of the 3rd node of kth element; k: kth-element; i=1,x1-coordinate; i=2,x2-coordinate;
% ...xm3: is for quadratic element;
% element12(jr,i): i=1: the element behind node jr; i=2: the element ahead of node jr;
% ... i=3: the element containing node jr as the mid-point;
% theta_x(k): direction angle of the kth element;
% theta_xi(j): direction angle of the jth node, which may be different from theta(k) for the kth element;
% ds(k): length of the kth element;
%
% t_hat(i,j): prescribed traction value in the ith direction of the jth node;
% u_hat(i,j): prescribed displacement value in the ith direction of the jth node;
% t_vc(i,j): variable condition of traction in the ith direction of the jth node;
% ...t_vc(i,j)=1: is an unknown variable; t_vc(i,j)=0: is known and not a variable;
% u_vc(i,j): variable condition of displacement in the ith direction of the jth node;
% ...u_vc(i,j)=1: is an unknown variable; u_vc(i,j)=0: is known and not a variable;
%
% Followings are the global variables generated in this function for couple problems:
% tc_hat(j): prescribed corner force of the jth corner;
% uc_hat(j): prescribed corner deflection of the jth croner;
% tc_vc(j): variable condition of corner force in the jth corner node;
% ...tc_vc(j)=1: is an unknown variable; tc_vc(j)=0: is known and not a variable;
15‐14
% uc_vc(j): variable condition of deflection in the jth corner node;
% ...uc_vc(j)=1: is an unknown variable; uc_vc(j)=0: is known and not a variable;
% ci: ci(cnode)=0: cnode is a regular node, ci(cnode)=1: cnode is a node behind of the kth-corner;
% ....ci(cnode)=2: cnode is a node ahead of the kth-corner;
global dim E0 xnL node12 bc h0 corner12 ncorner Btype nelement nnode % generated in Main & BEMbankB;
global xm1 xm2 xm3 element12 theta_x theta_xi ds t_hat u_hat t_vc u_vc ua_hat u0_hat v0_hat % generated;
global tc_hat uc_hat tc_vc uc_vc ci % generated: for coupled_BEM;
if Btype==1 | Btype==3, element12=zeros(nnode,2); % linear element;
elseif Btype==2, element12=zeros(nnode,3); end % quadratic element;
for in=1:nnode, theta_xi(in)=100; end % 100 is an initial value which will not occur in a normal theta;
for k=1:nelement
node1=node12(k,1); node2=node12(k,2);
xm1(:,k)=xnL(node1,:); % xm1(i,k): position of the 1st node of kth element;
xm2(:,k)=xnL(node2,:); % xm2(i,k): position of the 2nd node of kth element;
switch Btype
case {1,3} % linear element or linear element with cubic deflection
element12(node2,1)=k; element12(node1,2)=k; % construct element12(jr,i):
dx=xm2(1,k)-xm1(1,k); dy=xm2(2,k)-xm1(2,k);
ds(k)=sqrt(dx^2+dy^2); % ds(k): length of the kth element;
cost(k)=dx/ds(k); sint(k)=dy/ds(k); % cost(k),sint(k): cos_theta and sin_theta of the kth element;
if sint(k)>=0, theta_x(k)=acos(cost(k));
else theta_x(k)=2*pi-acos(cost(k)); end
% construct theta_xi for all nodes which may be different from theta for elements
if theta_xi(node1)==100, theta_xi(node1)=theta_x(k); % if node1 is considered in the first time else theta_xi(node1)=(theta_x(k)+theta_xi(node1))/2; end
if theta_xi(node2)==100, theta_xi(node2)=theta_x(k);
else theta_xi(node2)=(theta_x(k)+theta_xi(node2))/2; end case 2 % quadratic element
node3=node12(k,3); xm3(:,k)=xnL(node3,:); % xm3(i,k): position of the 3rd node of kth element;
element12(node3,1)=k; element12(node1,2)=k; element12(node2,3)=k; % construct element12(jr,i);
end end
%
for j=1:nnode for i=1:dim
t_hat(i,j)=0; u_hat(i,j)=0; t_vc(i,j)=1; u_vc(i,j)=1;
if Ltype1==4, u0_hat(i,j)=0; v0_hat(i,j)=0; ua_hat(i,j)=0; end switch bc(j,i)
case 0
t_vc(i,j)=0; % traction boundary t_hat(i,j)=bc(j,i+dim); % prescribed traction value
if Ltype1==4, % initial displacement and velocity for dynamic problems;
u0_hat(i,j)=bc(j,i+2*dim); v0_hat(i,j)=bc(j,i+3*dim); end case 1
u_vc(i,j)=0; % displacement boundary u_hat(i,j)=bc(j,i+dim); % prescribed displacement value
if Ltype1==4, ua_hat(i,j)=bc(j,i+2*dim); end % prescribed acceleration for dynamic problems;
otherwise
disp('wrong input of boundary condition ') end
end end
t_hat=t_hat/E0; % 2D_BEM, nondimensionalization;
%
if Ltype1==5 % coupled stretching-bending analysis;
ci=zeros(1,nnode); % set default value zero to corner index ci, ci=0 for regular node;
for kc=1:ncorner
cnode=corner12(kc,1); ci(cnode)=1; % ci=1; the node behind of the kth-corner;
cnode=corner12(kc,2); ci(cnode)=2; % ci=2; the node ahead of the kth-corner;
end
t_hat(1:2,:)=t_hat(1:2,:)/h0; t_hat(3:4,:)=t_hat(3:4,:)/(h0^2); % nondimension (been divided by E0 for 2D_BEM);
for j=1:ncorner
jc=nnode+j; tc_hat(j)=0; uc_hat(j)=0; tc_vc(j)=1; uc_vc(j)=1;
switch bc(jc,1) case 0
tc_vc(j)=0; % traction boundary
tc_hat(j)=bc(jc,2)/(E0*h0*h0); % nondimensionalized prescribed traction value case 1
uc_vc(j)=0; % displacement boundary uc_hat(j)=bc(jc,2); % prescribed displacement value otherwise
disp('wrong input of boundary condition ') end
end end