推导Eshelby夹杂解析解的maple程序

计算内点的弹性场

#计算内点:restart:
> with(linalg):
> assume(a>0):
> a1:=a:
> a2:=a:
> a3:=a:
> engstrain:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> engstrain[i,i]:=10^(-6);
> od: 
> deta:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> deta[i,i]:=1:
> end do:
> C:=array(1..3,1..3,1..3,1..3):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> C[i,j,k,l]:=2*mu*nu/(1-2*nu)*deta[i,j]*deta[k,l]+mu*deta[i,k]*deta[j,l]+mu*deta[i,l]*deta[j,k]:
> od: od: od: od:
> #下面计算I积分
> iI:=array(1..3,sparse):
> ii:=array(1..3,1..3,symmetric):
> y:=(a1^2+s)^(1/2)*(a2^2+s)^(1/2)*(a3^2+s)^(1/2):
> for i from 1 to 3 do:
> iI[i]:=int(2*Pi*a1*a2*a3/((a||i^2+s)*y),s=0..infinity):
> od:
> for i from 1 to 3 do:
> for j from i to 3 do:
> ii[i,j]:=int(2*Pi*a1*a2*a3/((a||i^2+s)*(a||j^2+s)*y),s=0..infinity):
> od: od:
> ####
> S:=array(1..3,1..3,1..3,1..3,sparse):
> for i from 1 to 3 do:
> S[i,i,i,i]:=3/(8*Pi*(1-nu))*a||i^2*ii[i,i]+(1-2*nu)/(8*Pi*(1-nu))*iI[i]:
> od:
> S[1,1,2,2]:=1/(8*Pi*(1-nu))*a||2^2*ii[1,2]-(1-2*nu)/(8*Pi*(1-nu))*iI[1]:
> S[1,1,3,3]:=1/(8*Pi*(1-nu))*a||3^2*ii[1,3]-(1-2*nu)/(8*Pi*(1-nu))*iI[1]:
> S[2,2,1,1]:=1/(8*Pi*(1-nu))*a||1^2*ii[2,1]-(1-2*nu)/(8*Pi*(1-nu))*iI[2]:
> S[2,2,3,3]:=1/(8*Pi*(1-nu))*a||3^2*ii[2,3]-(1-2*nu)/(8*Pi*(1-nu))*iI[2]:
> S[3,3,2,2]:=1/(8*Pi*(1-nu))*a||2^2*ii[3,2]-(1-2*nu)/(8*Pi*(1-nu))*iI[3]:
> S[3,3,1,1]:=1/(8*Pi*(1-nu))*a||2^2*ii[3,1]-(1-2*nu)/(8*Pi*(1-nu))*iI[3]:
> S[1,2,1,2]:=1/(16*Pi*(1-nu))*(a||1^2+a||2^2)*ii[1,2]+(1-2*nu)/(16*Pi*(1-nu))*(iI[1]+iI[2]):
> S[2,3,2,3]:=1/(16*Pi*(1-nu))*(a||2^2+a||3^2)*ii[2,3]+(1-2*nu)/(16*Pi*(1-nu))*(iI[2]+iI[3]):
> S[1,3,1,3]:=1/(16*Pi*(1-nu))*(a||1^2+a||3^2)*ii[1,2]+(1-2*nu)/(16*Pi*(1-nu))*(iI[1]+iI[3]):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for m from 1 to 3 do:
> for n from 1 to 3 do:
> if(S[i,j,m,n]<>0) then
> S[j,i,m,n]:=S[i,j,m,n]:
> S[i,j,n,m]:=S[i,j,m,n]:
> fi:
> if(S[j,i,m,n]<>0) then
> S[i,j,n,m]:=S[j,i,m,n]:
> fi:
> od: od: od: od:
> strain:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for m from 1 to 3 do:
> for n from 1 to 3 do:
> strain[i,j]:=S[i,j,m,n]*engstrain[m,n]+strain[i,j]:
> od: od: od: od:
> estrain:=evalm(strain-engstrain):
> stress:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for m from 1 to 3 do:
> for n from 1 to 3 do:
> stress[i,j]:=C[i,j,m,n]*estrain[m,n]+stress[i,j]:
> od: od: od: od:
> print(simplify(stress)):

计算外点的弹性场

restart:with(linalg):
> a1:=a:
> a2:=a:
> a3:=a:
> engstrain:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> engstrain[i,i]:=10^(-6):
> od: 
> deta:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> deta[i,i]:=1:
> end do:
> C:=array(1..3,1..3,1..3,1..3):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> C[i,j,k,l]:=2*mu*nu/(1-2*nu)*deta[i,j]*deta[k,l]+mu*deta[i,k]*deta[j,l]+mu*deta[i,l]*deta[j,k]:
> od: od: od: od:
> iI:=array(1..3,sparse):
> ii:=array(1..3,1..3):
> y:=(a1^2+s)^(1/2)*(a2^2+s)^(1/2)*(a3^2+s)^(1/2):
> for i from 1 to 3 do:
> iI[i]:=int(2*Pi*a1*a2*a3/((a||i^2+s)*y),s=lambda..infinity):
> od:
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> ii[i,j]:=int(2*Pi*a1*a2*a3/((a||i^2+s)*(a||j^2+s)*y),s=lambda..infinity):
> od: od:
> #这里先计算lambda的导数:
> xx:=array(1..3,[x1,x2,x3]):
> deI:=array(1..3,1..3):#iI的第一阶导数
> for i from 1 to 3 do:
> for l from 1 to 3 do:
> deI[i,l]:=diff(iI[i],lambda)*((2*xx[l]/(a||l^2+lambda))/(xx[1]^2/(a1^2+lambda)^2+xx[2]^2/(a2^2+lambda)^2+xx[3]^2/(a3^2+lambda)^2)):
> od: od:
> deI2:=array(1..3,1..3,1..3,sparse):#Ii的第二阶导数
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for l from 1 to 3 do:
> deI2[i,j,l]:=diff(deI[i,j],lambda)*((2*xx[l]/(a||l^2+lambda))/(xx[1]^2/(a1^2+lambda)^2+xx[2]^2/(a2^2+lambda)^2+xx[3]^2/(a3^2+lambda)^2)):
> od: od: od:
> deII1:=array(1..3,1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for l from 1 to 3 do:
> deII1[i,j,l]:=diff(ii[i,j],lambda)*((2*xx[l]/(a||l^2+lambda))/(xx[1]^2/(a1^2+lambda)^2+xx[2]^2/(a2^2+lambda)^2+xx[3]^2/(a3^2+lambda)^2)):
> od: od: od:
> deII2:=array(1..3,1..3,1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> deII2[i,j,k,l]:=diff(deII1[i,j,k],lambda)*((2*xx[l]/(a||l^2+lambda))/(xx[1]^2/(a1^2+lambda)^2+xx[2]^2/(a2^2+lambda)^2+xx[3]^2/(a3^2+lambda)^2)):
> od: od: od: od:
> S:=array(1..3,1..3,1..3,1..3):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> S[i,j,k,l]:=deta[i,j]*deta[k,l]*(2*nu*iI[i]-iI[k]+a||i^2*ii[k,i])+(deta[i,k]*deta[j,l]+deta[j,k]*deta[i,l])*(a||i^2*ii[i,j]-iI[j]+(1-nu)*(iI[k]+iI[l])):
> end do:
> end do:
> end do:
> end do:
> DD:=array(1..3,1..3,1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> DD[i,j,k,l]:=S[i,j,k,l]+2*nu*deta[k,l]*xx[i]*deI[i,j]+(1-nu)*(deta[i,l]*xx[k]*deI[k,j]+deta[j,l]*xx[k]*deI[k,i]+deta[i,k]*xx[l]*deI[l,j]+deta[j,k]*xx[l]*deI[l,i])-deta[i,j]*xx[k]*(deI[k,l]-a||i^2*deII1[k,i,l])-(deta[i,k]*xx[j]+deta[j,k]*xx[i])*(deI[j,l]-a||i^2*deII1[i,j,l])-(deta[i,l]*xx[j]+deta[j,l]*xx[i])*(deI[j,k]-a||i^2*deII1[i,j,k])-xx[i]*xx[j]*(deI2[j,l,k]-a||i^2*deII2[i,j,l,k]):
> od: od: od: od:
> strain:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> DD[i,j,k,l]:=DD[i,j,k,l]/(8*Pi*(1-nu)):
> od: od: od: od:
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> strain[i,j]:=DD[i,j,k,l]*engstrain[k,l]+strain[i,j]:
> od: od: od: od:
> stress:=array(1..3,1..3,sparse):
> for i from 1 to 3 do:
> for j from 1 to 3 do:
> for k from 1 to 3 do:
> for l from 1 to 3 do:
> stress[i,j]:=C[i,j,k,l]*(strain[k,l])+stress[i,j]:
> od: od: od: od:
> #print(C):
> print(simplify(strain)):
> #print(deI):
> #print(ii):
> #print("//////////////////////////////////////////////////////////////////////////////////////"):
> #print(iI):
> #print(deI):