198 lines
5.8 KiB
Mathematica
198 lines
5.8 KiB
Mathematica
|
function [] = XCOR1Db()
|
||
|
|
% MATLAB based 1-D XFEM Solver
|
||
|
|
% J. Grogan (2012)
|
||
|
|
clear all
|
||
|
|
% Define Geometry
|
||
|
|
len=1.;
|
||
|
|
% Define Section Properties
|
||
|
|
rho=1.;
|
||
|
|
% Generate Mesh
|
||
|
|
numElem=4;
|
||
|
|
charlen=len/numElem;
|
||
|
|
ndCoords=linspace(0,len,numElem+1);
|
||
|
|
numNodes=size(ndCoords,2);
|
||
|
|
indx=1:numElem;
|
||
|
|
elemNodes(:,1)=indx;
|
||
|
|
elemNodes(:,2)=indx+1;
|
||
|
|
% dofs per node
|
||
|
|
ndof=2;
|
||
|
|
% initial interface position
|
||
|
|
dpos=0.6;
|
||
|
|
% Initial temperatures
|
||
|
|
Tnew=zeros(numNodes*2,1);
|
||
|
|
Bound=zeros(numNodes*2,1);
|
||
|
|
%storage
|
||
|
|
stored(1)=dpos;
|
||
|
|
for e=1:numElem
|
||
|
|
for n=1:2
|
||
|
|
crdn1=ndCoords(elemNodes(e,n));
|
||
|
|
if crdn1<=dpos
|
||
|
|
Tnew(2*elemNodes(e,n)-1)=1.;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
Bound(1)=1.;
|
||
|
|
% Define Time Step
|
||
|
|
dtime=0.01;
|
||
|
|
tsteps=10;
|
||
|
|
time=0.;
|
||
|
|
% penalty term
|
||
|
|
beta=100.;
|
||
|
|
% Loop through time steps
|
||
|
|
for ts=1:tsteps
|
||
|
|
eNodes=zeros(2*numNodes,1);
|
||
|
|
% Get interface velocity
|
||
|
|
d(1)=dpos+charlen;
|
||
|
|
d(2)=dpos+3*charlen/4;
|
||
|
|
d(3)=dpos+charlen/4;
|
||
|
|
d(4)=dpos;
|
||
|
|
for e=1:numElem
|
||
|
|
crdn1=ndCoords(elemNodes(e,1));
|
||
|
|
crdn2=ndCoords(elemNodes(e,2));
|
||
|
|
for j=1:4
|
||
|
|
if d(j)>=crdn1 && d(j)<crdn2
|
||
|
|
elen=abs(crdn2-crdn1);
|
||
|
|
ajacob=elen/2.;
|
||
|
|
point=(d(j)-crdn1)/ajacob-1.;
|
||
|
|
theta(1)=abs(crdn1-dpos)*sign(crdn1-dpos);
|
||
|
|
theta(2)=abs(crdn2-dpos)*sign(crdn2-dpos);
|
||
|
|
tmp1a=Tnew(elemNodes(e,1)*2-1);
|
||
|
|
tmp1b=Tnew(elemNodes(e,1)*2);
|
||
|
|
tmp2a=Tnew(elemNodes(e,2)*2-1);
|
||
|
|
tmp2b=Tnew(elemNodes(e,2)*2);
|
||
|
|
xi=point;
|
||
|
|
gm(1)=(1.-xi)/2.;
|
||
|
|
gm(3)=(1.+xi)/2.;
|
||
|
|
term=theta(1)*gm(1)+theta(2)*gm(3);
|
||
|
|
gm(2)=gm(1)*(abs(term)-abs(theta(1)));
|
||
|
|
gm(4)=gm(3)*(abs(term)-abs(theta(2)));
|
||
|
|
t(j)=gm(1)*tmp1a+gm(2)*tmp1b+gm(3)*tmp2a+gm(4)*tmp2b;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
vel=0.5*(0.5/charlen)*(2*t(1)+t(2)-t(3)-2*t(4));
|
||
|
|
vel=0.
|
||
|
|
% Update interface position
|
||
|
|
dpos=dpos+vel*dtime;
|
||
|
|
stored(ts+1)=dpos;
|
||
|
|
K=zeros(numNodes*ndof,numNodes*ndof);
|
||
|
|
M=zeros(numNodes*ndof,numNodes*ndof);
|
||
|
|
pforce=zeros(numNodes*ndof,1);
|
||
|
|
% Loop Through Elements
|
||
|
|
for e=1:numElem
|
||
|
|
Ke=zeros(2*ndof);
|
||
|
|
Me=zeros(2*ndof);
|
||
|
|
crdn1=ndCoords(elemNodes(e,1));
|
||
|
|
crdn2=ndCoords(elemNodes(e,2));
|
||
|
|
theta(1)=abs(crdn1-dpos)*sign(crdn1-dpos);
|
||
|
|
theta(2)=abs(crdn2-dpos)*sign(crdn2-dpos);
|
||
|
|
enr=2;
|
||
|
|
elen=abs(crdn2-crdn1);
|
||
|
|
ajacob=elen/2.;
|
||
|
|
if sign(theta(1))~=sign(theta(2))
|
||
|
|
% enriched element
|
||
|
|
eNodes(2*elemNodes(e,1))=1;
|
||
|
|
eNodes(2*elemNodes(e,2))=1;
|
||
|
|
enr=4;
|
||
|
|
% get interface position on element
|
||
|
|
point=(dpos-crdn1)/ajacob-1.;
|
||
|
|
% devide element for sub integration
|
||
|
|
len1=abs(-point-1.);
|
||
|
|
len2=abs(1.-point);
|
||
|
|
mid1=-1+len1/2.;
|
||
|
|
mid2=1-len2/2.;
|
||
|
|
gpx(1)=-(len1/2.)/sqrt(3.)+mid1;
|
||
|
|
gpx(2)=(len1/2.)/sqrt(3.)+mid1;
|
||
|
|
gpx(3)=-(len2/2.)/sqrt(3.)+mid2;
|
||
|
|
gpx(4)=(len2/2.)/sqrt(3.)+mid2;
|
||
|
|
w(1)=(len1/2.);
|
||
|
|
w(2)=(len1/2.);
|
||
|
|
w(3)=(len2/2.);
|
||
|
|
w(4)=(len2/2.);
|
||
|
|
else
|
||
|
|
% regular element - fix extra dofs
|
||
|
|
gpx(1)=-1/sqrt(3.);
|
||
|
|
gpx(2)=1/sqrt(3.);
|
||
|
|
w(1)=1.;
|
||
|
|
w(2)=1.;
|
||
|
|
end
|
||
|
|
% Loop Through Int Points
|
||
|
|
for i=1:enr;
|
||
|
|
c=gpx(i);
|
||
|
|
phi(1)=(1.-c)/2.;
|
||
|
|
phi(3)=(1.+c)/2.;
|
||
|
|
term=theta(1)*phi(1)+theta(2)*phi(3);
|
||
|
|
cond=1.;
|
||
|
|
spec=1.;
|
||
|
|
phi(2)=phi(1)*(abs(term)-abs(theta(1)));
|
||
|
|
phi(4)=phi(3)*(abs(term)-abs(theta(2)));
|
||
|
|
phic(1)=-0.5;
|
||
|
|
phic(3)=0.5;
|
||
|
|
dterm=sign(term)*(phic(1)*theta(1)+phic(3)*theta(2));
|
||
|
|
phic(2)=phic(1)*(abs(term)-abs(theta(1)))+phi(1)*dterm;
|
||
|
|
phic(4)=phic(3)*(abs(term)-abs(theta(2)))+phi(3)*dterm;
|
||
|
|
phix(1)=phic(1)/ajacob;
|
||
|
|
phix(2)=phic(2)/ajacob;
|
||
|
|
phix(3)=phic(3)/ajacob;
|
||
|
|
phix(4)=phic(4)/ajacob;
|
||
|
|
we=ajacob*w(i);
|
||
|
|
Ke=Ke+we*cond*phix'*phix;
|
||
|
|
Me=Me+(we*rho*spec*phi'*phi)/dtime;
|
||
|
|
end
|
||
|
|
% Add penalty term and get temp gradient on interface
|
||
|
|
if enr==4;
|
||
|
|
xi=point;
|
||
|
|
gm(1)=(1.-xi)/2.;
|
||
|
|
gm(3)=(1.+xi)/2.;
|
||
|
|
term=theta(1)*gm(1)+theta(2)*gm(3);
|
||
|
|
gm(2)=gm(1)*(abs(term)-abs(theta(1)));
|
||
|
|
gm(4)=gm(3)*(abs(term)-abs(theta(2)));
|
||
|
|
pen=beta*(gm'*gm);
|
||
|
|
pfL=beta*1*gm';
|
||
|
|
Ke=Ke+pen;
|
||
|
|
else
|
||
|
|
pen=zeros(4);
|
||
|
|
pfL=zeros(4,1);
|
||
|
|
end
|
||
|
|
% Assemble Global Matrices
|
||
|
|
gnum=2.*elemNodes(e,1)-1.;
|
||
|
|
for i=1:4;
|
||
|
|
for j=1:4;
|
||
|
|
K(gnum+j-1,gnum+i-1)=K(gnum+j-1,gnum+i-1)+Ke(j,i);
|
||
|
|
M(gnum+j-1,gnum+i-1)=M(gnum+j-1,gnum+i-1)+Me(j,i);
|
||
|
|
end
|
||
|
|
pforce(gnum+i-1)=pforce(gnum+i-1)+pfL(i);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
A=K+M;
|
||
|
|
Sub=A*Bound;
|
||
|
|
RHS=M*Tnew-Sub+pforce;
|
||
|
|
% Apply Boundary Conditions
|
||
|
|
Biindex=0.;
|
||
|
|
for i=1:ndof*numNodes;
|
||
|
|
if Bound(i)==0.;
|
||
|
|
Biindex=Biindex+1;
|
||
|
|
RHSR(Biindex)=RHS(i);
|
||
|
|
jindex=0;
|
||
|
|
for j=1:ndof*numNodes;
|
||
|
|
if Bound(j)==0.;
|
||
|
|
jindex=jindex+1;
|
||
|
|
AR(Biindex,jindex)=A(i,j);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
%Solve
|
||
|
|
Tnewr=(AR^-1)*RHSR';
|
||
|
|
% Restore Matrices
|
||
|
|
Biindex=0;
|
||
|
|
for i=1:ndof*numNodes;
|
||
|
|
if Bound(i)==0.;
|
||
|
|
Biindex=Biindex+1;
|
||
|
|
Tnew(i)=Tnewr(Biindex);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
Tnew
|
||
|
|
end
|
||
|
|
stored';
|