This is the python routine I used to calculate the linear viscoelasticity (storage and loss shear moduli) of worm-like microstructures define in the paper on linear viscoelasticity of IPNs. This routine may be used for any two-phase materials with the microstructure of your choice.
You need two input.txt (phase1.txt and phase2.txt) files for the behavior constitutive phase. An example of these data and the python code may be dowloaded here.
Note that I am a poor code writer and if you are decent at it you can probably do better but this provide you with a ready to use code
"""
Created June 22, 2019 by Julie DIANI (julie.diani@polytechnique.edu)
FFT Homogeneization program for linear elastic and viscoelastic materials
Method introduced by H. P. Suquet (1998), A numerical method for computing
the overall response of nonlinear composites with complex microstructure.
Comp. Meth. Appl. Mech. Eng., 157, 69--94
Accelerated scheme from D. Eyre et G. Milton (1999), A fast numerical scheme
for computing the response of composites using grid refinement, Eur. Phys. J.,
6, 41-47.
Before writing the program below, I have used read programm from T.W.J. Geus
https://github.com/tdegeus/GooseFFT/blob/master/small-strain/linear-elasticity.py
T.W.J. de Geus, J. Vondˇrejc, J. Zeman, Peerlings, R.H.J. Peerlings, M.G.D. Geers
Comput. Methods Appl. Mech. Eng., 2017, 318:412430
doi: 10.1016/j.cma.2016.12.032.
I have used:
-Voigt notation to reduce the indices
-pyfftw package to be able to use FFT for other Ng than 2**n
-resize to flap the (Ng,Ng,Ng) matrices into vectors(Ng**3)
to reduce memory cost
If you use part or the entire program, please cite reference:
J. Diani, E. Strauch-Hausser, 2021. Linear viscoelasticity of an acrylate IPN, analysis and micromechanics modeling. Soft matter, 17, 7341-7349. https://doi.org/10.1039/d1sm00808k.
"""
###############################################################################
import pyfftw # Works with version 1.3 of scipy pip install scipy==1.3
import numpy as np
#import scipy.sparse.linalg as sp
import itertools
import time
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#import matplotlib.animation as animation
###############################################################################
## Mechanical behavior two viscoelastic material
###############################################################################
# Experimental data
# linear viscoelastic of both materials
# Material 1
MA1=np.loadtxt('Phase1.txt') # Temp/G'/G" for material 1
# Material 2
MA2=np.loadtxt('Phase2.txt') # Temp/G'/G" for material 2
# Assumption: Constant elastic bulk modulus, equal for both phase
Km=3000;
###############################################################################
# Possibly putting results in an output_file
file_out=open('FFT_output.txt', 'a')
file_out.write('Temp Storage shear modulus Loss shear'+'\n')
file_out.write('ligne 1 : G23 / ligne 2 : G13 / ligne 3 : G12'+'\n')
###############################################################################
# nbpoint:number of data point in between two calculations
# mo: number of point that will be calculated
###############################################################################
nbpoint=50
mo = int(min(len(MA2),len(MA1))/nbpoint)
###############################################################################
# Space dimension
ndim=3
# Lattice dimension Ngxndim
Ng=32
Nmg2=(Ng-1)**2
Nmg3=(Ng-1)**3
Ng3=Ng**3
###############################################################################
# Some useful functions
delta = lambda k,l: np.float(k==l) # \delta_{ij}
eij = np.array([[0,0],[1,1],[2,2],[3,3],[2,3],[1,3],[1,2]])
# Defines some fft functions
fft = lambda x: np.fft.fftshift(pyfftw.interfaces.numpy_fft.fftn(np.fft.ifftshift(x),[Ng,Ng,Ng]))
ifft = lambda x: np.fft.fftshift(pyfftw.interfaces.numpy_fft.ifftn(np.fft.ifftshift(x),[Ng,Ng,Ng]))
# If not using pyfftw but fft needs to have a grid size Ng=2^n
#fft = lambda x : np.fft.fftshift(np.fft.fftn (np.fft.ifftshift(x),[Ng,Ng,Ng]))
#ifft = lambda x : np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(x),[Ng,Ng,Ng]))
###############################################################################
## microstructure definition
## Exemple with heterogeneous 2 phases
# Target volume fraction vf of phase 2
###############################################################################
vf=0.3
nbpart=int(vf*Ng3) # number of pixel for phase 2
###############################################################################
# Example of a random microstructure
###############################################################################
# initialization microstructure
micro = int(np.zeros(1))*np.mgrid[:Ng, :Ng, :Ng][0]
#
while (np.sum(micro)<nbpart):
p=random.randint(0, Ng-2)
q=random.randint(0, Ng-2)
l=random.randint(0, Ng-2)
micro[p,q,l]=1
# for the symmetry of the microstructure
if p==0:
micro[Ng-1,q,l]=1
if p==0 and q==0:
micro[Ng-1,Ng-1,l]=1
if p==0 and l==0:
micro[Ng-1,q,Ng-1]=1
if p==0 and q==0 and l==0:
micro[Ng-1,Ng-1,Ng-1]=1
if q==0:
micro[p,Ng-1,l]=1
if q==0 and l==0:
micro[p,Ng-1,Ng-1]=1
if l==0:
micro[p,q,Ng-1]=1
# Restructuration of the microstructure
micro2=micro.reshape(-1,1)
# Actual volume fraction
volfrac=np.sum(micro2)/Ng3
print('volfrac=',volfrac)
## 3D Display of a microstructure
x, y, z = np.indices((Ng,Ng,Ng))
# combine the objects into a single boolean array
cube1 = micro[x,y,z]<1
cube2 = micro[x,y,z]>0.5
# combine the objects into a single boolean array
voxels = cube1 | cube2
# set the colors of each object
colors = np.empty(voxels.shape, dtype=object)
colors[cube1] = 'black'
colors[cube2] = 'white'
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.voxels(voxels, facecolors=colors, edgecolor='black', linewidth=0.1)
ax.axis('off')
#ax.set_aspect('equal')
plt.savefig('Micro_random.png', dpi=300)
#plt.show()
###############################################################################
# Ghs Storate and Ghl loss shear moduli of the heterogeneous material at
# temperatures Temph
###############################################################################
Ghs=np.zeros(mo)
Ghl=np.zeros(mo)
Temph=np.zeros(mo)
###############################################################################
# Clock to print the duration of calculations
# Duration parameter if needed
#t0= time.clock()
# loop lp on the number mo of homogeneized behaviors
# Initial temperature of interest
Temp_init=-80.0
# initial temperature
index_init = [i for i, j in enumerate(MA2[:,0]) if Temp_init-0.2 <j ][0]
for lp in range(mo):
print('loop #',lp,'/',mo)
lk=(nbpoint-1)*lp
Temph[lp]=MA2[lk+index_init,0]
ima1 = [i for i, j in enumerate(MA1[:,0]) if Temph[lp]-0.2 <j ][0]
ima2 = [i for i, j in enumerate(MA2[:,0]) if Temph[lp]-0.2 <j ][0]
###############################################################################
## viscoelastic behavior Shear modulus (complex) = mu
## Material 1 microstructure at 0
## Material 2 microstructure at 1
mu1 = MA1[ima1,1]+MA1[ima1,2]*1.j
mu2 = MA2[ima2,1]+MA2[ima2,2]*1.j
# Lamé constant = complex
lambda1=Km-(2*mu1/3)
lambda2=Km-(2*mu2/3)
############################################
# Elastic isotropic reference material behavior C0
# The reference material does not have to be viscoelastic
############################################
mu0=np.sqrt(mu1.real*mu2.real)
K0=Km # = np.sqrt(K1*K2)
lambda0=K0-(2*mu0/3)
lm0=(lambda0+mu0)/(mu0*(lambda0+2*mu0))
############################################
# Local behavior of the microstructure C(x)
mulocal=micro2*mu2+(1-micro2)*mu1
llocal=micro2*lambda2+(1-micro2)*lambda1
############################################
## Define usuful vector
freq = np.arange(-Ng/2.,+Ng/2.)
############################################
# 4-order tensor Sijkl=(C(x)+C0)^(-1)ijkl
# S11=S_1111 S12=S_1122 S14=S_1212
############################################
S11=2*(4*(mulocal+mu0))**(-1)-\
0.5*(llocal+lambda0)/((mulocal+mu0)*\
(3*(llocal+lambda0)+2*(mulocal+mu0)))
S12=-0.5*(llocal+lambda0)/((mulocal+mu0)*(3*(llocal+lambda0)+\
2*(mulocal+mu0)))
S44=(4*(mulocal+mu0))**(-1)
#####################################################
# Green tensor G0_ijkl function of reference behavior
# This calculations is necessary at each step
#####################################################
# Initialization of useful tensors
for l in range(1,7):
for k in range(l,7):
exec("G{}{}4f = np.zeros([Ng,Ng,Ng])".format(l,k))
exec("G{}{}4 = G114f.reshape(-1,1)".format(l,k))
#
for l in range(1,7):
exec("Epsn{}2 = G114f.reshape(-1,1)".format(l))
exec("Taun{}2 = G114f.reshape(-1,1)".format(l))
exec("Epsc{}2f = G114f.reshape(-1,1)".format(l))
exec("Epsi{}2f = G114f.reshape(-1,1)".format(l))
exec("Ec{}2f = np.zeros([Ng,Ng,Ng])".format(l))
exec("Tn{}2 = np.zeros([Ng,Ng,Ng])".format(l))
#
for x,y,z in itertools.product(range(Ng), repeat=3):
q = np.array([freq[x], freq[y], freq[z]])
if not q.dot(q) == 0:
for n1 in range(1,7):
for n2 in range(n1,7):
k=eij[n1,0]-1
kj=eij[n1,1]-1
l=eij[n2,0]-1
m=eij[n2,1]-1
exec("G{}{}4f[x,y,z] = (4*mu0*np.vdot(q,q))**(-1)*(\
delta(k,l)*q[kj]*q[m]+delta(k,m)*q[kj]*q[l]+\
delta(kj,l)*q[k]*q[m]+delta(kj,m)*q[k]*q[l] )-\
(lambda0+mu0)/(mu0*(lambda0+2*mu0)*np.vdot(q,q)**2)*\
q[k]*q[kj]*q[l]*q[m]".format(n1,n2))
#
for n1 in range(1,7):
for n2 in range(n1,7):
exec("G{}{}4 = G{}{}4f.reshape(-1,1)".format(n1,n2,n1,n2))
exec("del G{}{}4f ".format(n1,n2))
###############################################################################
# Loading tests 3 shear tests \eps_ij=dep
# iterative scheme until err<=errmax
# errmax may be decreased if necessary
###############################################################################
errmax=2E-6
dep=0.01
###############################################################################
# Three shearing tests according to 23, 13 and 12
###############################################################################
# Shearing test 23
nk=4
# initialization step \epsilon=Epsn2
exec("Epsn{}2=dep+Epsn{}2".format(nk,nk))
#
# Error initialization
err=1
#
# Itteration scheme while err>=errmax :
# initalize a counter for those who are interested in the number of steps
# needed to reach convergence
kcompt=0
while err>errmax:
Taun12=2*(mulocal-mu0)*Epsn12+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun22=2*(mulocal-mu0)*Epsn22+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun32=2*(mulocal-mu0)*Epsn32+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun42=2*(mulocal-mu0)*Epsn42
Taun52=2*(mulocal-mu0)*Epsn52
Taun62=2*(mulocal-mu0)*Epsn62
#
# Calculate the FFT of the polarization field
for l in range(1,7):
exec("Tn{}2=Taun{}2.reshape((Ng,Ng,Ng))".format(l,l))
exec("Taun{}2f=fft(Tn{}2).reshape(-1,1)".format(l,l))
exec("Tn{}2 = np.zeros([Ng,Ng,Ng])".format(l))
exec("Taun{}2 = Tn{}2.reshape(-1,1)".format(l,l))
#
# Calculate \epsilon_comp
Epsc12f=-(G114*Taun12f+G124*Taun22f+G134*Taun32f+\
2*(G144*Taun42f+G154*Taun52f+G164*Taun62f))
Epsc22f=-(G124*Taun12f+G224*Taun22f+G234*Taun32f+\
2*(G244*Taun42f+G254*Taun52f+G264*Taun62f))
Epsc32f=-(G134*Taun12f+G234*Taun22f+G334*Taun32f+\
2*(G344*Taun42f+G354*Taun52f+G364*Taun62f))
Epsc42f=-(G144*Taun12f+G244*Taun22f+G344*Taun32f+\
2*(G444*Taun42f+G454*Taun52f+G464*Taun62f))
Epsc52f=-(G154*Taun12f+G254*Taun22f+G354*Taun32f+\
2*(G454*Taun42f+G554*Taun52f+G564*Taun62f))
Epsc62f=-(G164*Taun12f+G264*Taun22f+G364*Taun32f+\
2*(G464*Taun42f+G564*Taun52f+G664*Taun62f))
# Calculate the strain at step i+1 in the Fourier space for [0,0,0]
# which is the average value of the applied strain field
for l in range(1,7):
exec("Ec{}2f=Epsc{}2f.reshape((Ng,Ng,Ng))".format(l,l))
exec("Ec{}2f[Ng//2,Ng//2,Ng//2]=0.0".format(l))
#
exec("Ec{}2f[Ng//2,Ng//2,Ng//2]=dep*Ng**ndim".format(nk))
#
# Inverse FFT then calculate \Epsn2
for l in range(1,7):
exec("Epsc{}2=ifft(Ec{}2f).reshape(-1,1)".format(l,l))
exec("Ec{}2f = np.zeros([Ng,Ng,Ng])".format(l))
exec("Epsc{}2f = np.zeros([Ng,Ng,Ng]).reshape(-1,1)".format(l))
#
#
Epsi12=Epsn12+2*(S11*((2*mu0+lambda0)*(Epsc12-Epsn12)+lambda0*\
(Epsc22-Epsn22+Epsc32-Epsn32))+S12*(\
(2*mu0+lambda0)*(Epsc22-Epsn22)+\
lambda0*(Epsc12-Epsn12+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc32-Epsn32)+\
lambda0*(Epsc22-Epsn22+Epsc12-Epsn12) ) )
Epsi22=Epsn22+2*(S11*((2*mu0+lambda0)*(Epsc22-Epsn22)+lambda0*\
(Epsc12-Epsn12+Epsc32-Epsn32))+S12*(\
(2*mu0+lambda0)*(Epsc12-Epsn12)+\
lambda0*(Epsc22-Epsn22+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc32-Epsn32)+\
lambda0*(Epsc22-Epsn22+Epsc12-Epsn12) ) )
Epsi32=Epsn32+2*(S11*((2*mu0+lambda0)*(Epsc32-Epsn32)+lambda0*\
(Epsc12-Epsn12+Epsc22-Epsn22))+S12*(\
(2*mu0+lambda0)*(Epsc12-Epsn12)+\
lambda0*(Epsc22-Epsn22+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc22-Epsn22)+\
lambda0*(Epsc32-Epsn32+Epsc12-Epsn12) ) )
Epsi42=Epsn42+8*S44*mu0*(Epsc42-Epsn42)
Epsi52=Epsn52+8*S44*mu0*(Epsc52-Epsn52)
Epsi62=Epsn62+8*S44*mu0*(Epsc62-Epsn62)
#
for l in range(1,7):
exec("Epsn{}2=Epsi{}2".format(l,l))
exec("Epsi{}2 = Ec{}2f.reshape(-1,1)".format(l,l))
# Evaluation of th error err
err=(1/dep)*np.sqrt(np.abs( np.mean( ( (Epsn12-Epsc12)**2+ \
(Epsn22-Epsc22)**2+(Epsn32-Epsc32)**2+\
2*(Epsn42-Epsc42)**2+2*(Epsn52-Epsc52)**2+\
2*(Epsn62-Epsc62)**2 ) ) ) )
#
# increment the counter
# kcompt+=1
# print('err',err)
#
# Error < errmax
exec("sig=np.mean(2*mulocal*Epsn{}2)".format(nk))
G12=sig/(2*dep)
Ghs[lp]=G12.real
Ghl[lp]=G12.imag
res=[Temph[lp],Ghs[lp],Ghl[lp]]
str_q = str(res)[1 : -1]
file_out.write(str_q+'\n')
#t1 = time.clock() - t0
#print("Time elapsed: ", t1)
##############################################################################
# Shearing test along 31
nk=5
#
# initialization step \epsilon=Epsn2
exec("Epsn{}2=dep+Epsn{}2".format(nk,nk))
#
# Error initialization
err=1
#
# Itteration scheme while err>=errmax :
# initalize a counter
kcompt=0
while err>errmax:
Taun12=2*(mulocal-mu0)*Epsn12+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun22=2*(mulocal-mu0)*Epsn22+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun32=2*(mulocal-mu0)*Epsn32+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun42=2*(mulocal-mu0)*Epsn42
Taun52=2*(mulocal-mu0)*Epsn52
Taun62=2*(mulocal-mu0)*Epsn62
#
# Calculate the FFT of the polarization field
for l in range(1,7):
exec("Tn{}2=Taun{}2.reshape((Ng,Ng,Ng))".format(l,l))
exec("Taun{}2f=fft(Tn{}2).reshape(-1,1)".format(l,l))
exec("Tn{}2 = np.zeros([Ng,Ng,Ng])".format(l))
exec("Taun{}2 = Tn{}2.reshape(-1,1)".format(l,l))
#
# Calculate \epsilon_comp
Epsc12f=-(G114*Taun12f+G124*Taun22f+G134*Taun32f+\
2*(G144*Taun42f+G154*Taun52f+G164*Taun62f))
Epsc22f=-(G124*Taun12f+G224*Taun22f+G234*Taun32f+\
2*(G244*Taun42f+G254*Taun52f+G264*Taun62f))
Epsc32f=-(G134*Taun12f+G234*Taun22f+G334*Taun32f+\
2*(G344*Taun42f+G354*Taun52f+G364*Taun62f))
Epsc42f=-(G144*Taun12f+G244*Taun22f+G344*Taun32f+\
2*(G444*Taun42f+G454*Taun52f+G464*Taun62f))
Epsc52f=-(G154*Taun12f+G254*Taun22f+G354*Taun32f+\
2*(G454*Taun42f+G554*Taun52f+G564*Taun62f))
Epsc62f=-(G164*Taun12f+G264*Taun22f+G364*Taun32f+\
2*(G464*Taun42f+G564*Taun52f+G664*Taun62f))
# Calculate the strain at step i+1 in the Fourier space for [0,0,0]
# which is the average value of the applied strain field
for l in range(1,7):
exec("Ec{}2f=Epsc{}2f.reshape((Ng,Ng,Ng))".format(l,l))
exec("Ec{}2f[Ng//2,Ng//2,Ng//2]=0.0".format(l))
#
exec("Ec{}2f[Ng//2,Ng//2,Ng//2]=dep*Ng**ndim".format(nk))
#
# Inverse FFT then calculate \Epsn2
for l in range(1,7):
exec("Epsc{}2=ifft(Ec{}2f).reshape(-1,1)".format(l,l))
exec("Ec{}2f = np.zeros([Ng,Ng,Ng])".format(l))
exec("Epsc{}2f = np.zeros([Ng,Ng,Ng]).reshape(-1,1)".format(l))
#
#
Epsi12=Epsn12+2*(S11*((2*mu0+lambda0)*(Epsc12-Epsn12)+lambda0*\
(Epsc22-Epsn22+Epsc32-Epsn32))+S12*(\
(2*mu0+lambda0)*(Epsc22-Epsn22)+\
lambda0*(Epsc12-Epsn12+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc32-Epsn32)+\
lambda0*(Epsc22-Epsn22+Epsc12-Epsn12) ) )
Epsi22=Epsn22+2*(S11*((2*mu0+lambda0)*(Epsc22-Epsn22)+lambda0*\
(Epsc12-Epsn12+Epsc32-Epsn32))+S12*(\
(2*mu0+lambda0)*(Epsc12-Epsn12)+\
lambda0*(Epsc22-Epsn22+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc32-Epsn32)+\
lambda0*(Epsc22-Epsn22+Epsc12-Epsn12) ) )
Epsi32=Epsn32+2*(S11*((2*mu0+lambda0)*(Epsc32-Epsn32)+lambda0*\
(Epsc12-Epsn12+Epsc22-Epsn22))+S12*(\
(2*mu0+lambda0)*(Epsc12-Epsn12)+\
lambda0*(Epsc22-Epsn22+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc22-Epsn22)+\
lambda0*(Epsc32-Epsn32+Epsc12-Epsn12) ) )
Epsi42=Epsn42+8*S44*mu0*(Epsc42-Epsn42)
Epsi52=Epsn52+8*S44*mu0*(Epsc52-Epsn52)
Epsi62=Epsn62+8*S44*mu0*(Epsc62-Epsn62)
#
for l in range(1,7):
exec("Epsn{}2=Epsi{}2".format(l,l))
exec("Epsi{}2 = Ec{}2f.reshape(-1,1)".format(l,l))
# Evaluation of th error err
err=(1/dep)*np.sqrt(np.abs( np.mean( ( (Epsn12-Epsc12)**2+ \
(Epsn22-Epsc22)**2+(Epsn32-Epsc32)**2+\
2*(Epsn42-Epsc42)**2+2*(Epsn52-Epsc52)**2+\
2*(Epsn62-Epsc62)**2 ) ) ) )
#
# increment the counter
# kcompt+=1
# print('err',err)
#
# Error < errmax
exec("sig=np.mean(2*mulocal*Epsn{}2)".format(nk))
G12=sig/(2*dep)
Ghs[lp]=G12.real
Ghl[lp]=G12.imag
res=[Temph[lp],Ghs[lp],Ghl[lp]]
str_q = str(res)[1 : -1]
file_out.write(str_q+'\n')
# t1 = time.clock() - t0
# print("Time elapsed: ", t1)
##############################################################################
# Shearing test 12
nk=6
#
# initialization step \epsilon=Epsn2
exec("Epsn{}2=dep+Epsn{}2".format(nk,nk))
#
# Error initialization
err=1
#
# Itteration scheme while err>=errmax :
# initalize a counter
kcompt=0
while err>errmax:
Taun12=2*(mulocal-mu0)*Epsn12+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun22=2*(mulocal-mu0)*Epsn22+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun32=2*(mulocal-mu0)*Epsn32+\
(llocal-lambda0)*(Epsn12+Epsn22+Epsn32)
Taun42=2*(mulocal-mu0)*Epsn42
Taun52=2*(mulocal-mu0)*Epsn52
Taun62=2*(mulocal-mu0)*Epsn62
#
# Calculate the FFT of the polarization field
for l in range(1,7):
exec("Tn{}2=Taun{}2.reshape((Ng,Ng,Ng))".format(l,l))
exec("Taun{}2f=fft(Tn{}2).reshape(-1,1)".format(l,l))
exec("Tn{}2 = np.zeros([Ng,Ng,Ng])".format(l))
exec("Taun{}2 = Tn{}2.reshape(-1,1)".format(l,l))
#
# Calculate \epsilon_comp
Epsc12f=-(G114*Taun12f+G124*Taun22f+G134*Taun32f+\
2*(G144*Taun42f+G154*Taun52f+G164*Taun62f))
Epsc22f=-(G124*Taun12f+G224*Taun22f+G234*Taun32f+\
2*(G244*Taun42f+G254*Taun52f+G264*Taun62f))
Epsc32f=-(G134*Taun12f+G234*Taun22f+G334*Taun32f+\
2*(G344*Taun42f+G354*Taun52f+G364*Taun62f))
Epsc42f=-(G144*Taun12f+G244*Taun22f+G344*Taun32f+\
2*(G444*Taun42f+G454*Taun52f+G464*Taun62f))
Epsc52f=-(G154*Taun12f+G254*Taun22f+G354*Taun32f+\
2*(G454*Taun42f+G554*Taun52f+G564*Taun62f))
Epsc62f=-(G164*Taun12f+G264*Taun22f+G364*Taun32f+\
2*(G464*Taun42f+G564*Taun52f+G664*Taun62f))
# Calculate the strain at step i+1 in the Fourier space for [0,0,0]
# which is the average value of the applied strain field
for l in range(1,7):
exec("Ec{}2f=Epsc{}2f.reshape((Ng,Ng,Ng))".format(l,l))
exec("Ec{}2f[Ng//2,Ng//2,Ng//2]=0.0".format(l))
#
exec("Ec{}2f[Ng//2,Ng//2,Ng//2]=dep*Ng**ndim".format(nk))
#
# Inverse FFT then calculate \Epsn2
for l in range(1,7):
exec("Epsc{}2=ifft(Ec{}2f).reshape(-1,1)".format(l,l))
exec("Ec{}2f = np.zeros([Ng,Ng,Ng])".format(l))
exec("Epsc{}2f = np.zeros([Ng,Ng,Ng]).reshape(-1,1)".format(l))
#
#
Epsi12=Epsn12+2*(S11*((2*mu0+lambda0)*(Epsc12-Epsn12)+lambda0*\
(Epsc22-Epsn22+Epsc32-Epsn32))+S12*(\
(2*mu0+lambda0)*(Epsc22-Epsn22)+\
lambda0*(Epsc12-Epsn12+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc32-Epsn32)+\
lambda0*(Epsc22-Epsn22+Epsc12-Epsn12) ) )
Epsi22=Epsn22+2*(S11*((2*mu0+lambda0)*(Epsc22-Epsn22)+lambda0*\
(Epsc12-Epsn12+Epsc32-Epsn32))+S12*(\
(2*mu0+lambda0)*(Epsc12-Epsn12)+\
lambda0*(Epsc22-Epsn22+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc32-Epsn32)+\
lambda0*(Epsc22-Epsn22+Epsc12-Epsn12) ) )
Epsi32=Epsn32+2*(S11*((2*mu0+lambda0)*(Epsc32-Epsn32)+lambda0*\
(Epsc12-Epsn12+Epsc22-Epsn22))+S12*(\
(2*mu0+lambda0)*(Epsc12-Epsn12)+\
lambda0*(Epsc22-Epsn22+Epsc32-Epsn32)+\
(2*mu0+lambda0)*(Epsc22-Epsn22)+\
lambda0*(Epsc32-Epsn32+Epsc12-Epsn12) ) )
Epsi42=Epsn42+8*S44*mu0*(Epsc42-Epsn42)
Epsi52=Epsn52+8*S44*mu0*(Epsc52-Epsn52)
Epsi62=Epsn62+8*S44*mu0*(Epsc62-Epsn62)
#
for l in range(1,7):
exec("Epsn{}2=Epsi{}2".format(l,l))
exec("Epsi{}2 = Ec{}2f.reshape(-1,1)".format(l,l))
# Evaluation of th error err
err=(1/dep)*np.sqrt(np.abs( np.mean( ( (Epsn12-Epsc12)**2+ \
(Epsn22-Epsc22)**2+(Epsn32-Epsc32)**2+\
2*(Epsn42-Epsc42)**2+2*(Epsn52-Epsc52)**2+\
2*(Epsn62-Epsc62)**2 ) ) ) )
#
# increment the counter
# kcompt+=1
# print('err',err)
#
# Error < errmax
exec("sig=np.mean(2*mulocal*Epsn{}2)".format(nk))
#
G12=sig/(2*dep)
Ghs[lp]=G12.real
Ghl[lp]=G12.imag
res=[Temph[lp],Ghs[lp],Ghl[lp]]
str_q = str(res)[1 : -1]
file_out.write(str_q+'\n')
# t1 = time.clock() - t0
# print("Time elapsed: ", t1)
file_out.close()
###############################################################################
I add a python code to generate lattice two-phase periodic microstructures, either random or worm-like microstructures or beads embedded in a matrix. For the matrix with beads, you need to provide with the position of the center (‘centers1.txt’) – The file as well as the python code below may be dowloaded here.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 3 17:50:29 2018
@author: Julie Diani
Generate three kinds of two-phase microstructures
-Spheres in a matrix
-Random microstructure
-Matrix with a polymer network built by self-avoiding random walk
"""
import numpy as np
import itertools
import time
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
# Cubic grid Ng x Ng x Ng
# Discretization number
Ng=128
##############################################################################
## Microstructure 1
## Spherical particles in a matrix
##############################################################################
# List of the sphere centers
centers=np.loadtxt('centers1.txt')
nbcenter=len(centers)
# Target volume fraction
vf=0.54
# Radius of the beads in pixels
ray=np.int((vf*Ng**3/(nbcenter*(4/3)*np.pi))**(1/3))
ray2=ray**2
# microstructure
micro = int(np.zeros(1))*np.mgrid[:Ng, :Ng, :Ng][0]
micenters = int(np.zeros(1))*np.mgrid[:Ng, :Ng, :Ng][0]
# define the centers in pixel
for i in range(nbcenter):
p=np.int(centers[i,0]*Ng)
q=np.int(centers[i,1]*Ng)
l=np.int(centers[i,2]*Ng)
micenters[p,q,l]=1
for k1 in range(2*ray+1):
kp=(p-ray+k1)
for k2 in range(2*ray+1):
kq=(q-ray+k2)
for k3 in range(2*ray+1):
kl=(l-ray+k3)
long=(kp-p)**2+(kq-q)**2+(kl-l)**2
if (long<=ray2):
kkp=kp % Ng
kkq=kq % Ng
kkl=kl % Ng
micro[kkp,kkq,kkl]=1
# Calculate the actual volume fraction
vfreel=np.sum(micro)/Ng**3
display(np.sum(micro)/Ng**3)
# setup some generic data
for no in range(2):
fm=no;
x, y = np.mgrid[:Ng, :Ng]
Z = micro[x,fm,y]
# mask out the negative and positive values, respectively
fig = plt.figure(frameon=False)
# Represent the microstructure
pos_neg_clipped = plt.imshow(Z, cmap='binary',
interpolation='none')
# Add minorticks on the colorbar to make it easy to read the
# values off the colorbar.
exec("plt.savefig('M21_Ng512_fY{}.png', dpi=300)".format(fm))
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt')
fig = plt.figure()
# ims is a list of lists, each row is a list of artists to draw in the
# current frame; here we are just animating one artist, the image, in
# each frame
ims = []
for i in range(Ng):
im = plt.imshow(micro[i,x,y], cmap='binary',animated=True)
ims.append([im])
ani = animation.ArtistAnimation(fig, ims, interval=70, blit=True)
#ani = animation.ArtistAnimation(fig, ims, interval=300, blit=True,
# repeat_delay=1000)
exec("ani.save('dynamic_micro21_{}.mp4',writer='ffmpeg')".format(Ng))
plt.show()
##############################################################################
## Microstructure 2
## Random
##############################################################################
# Space dimension
ndim=3
# Lattice dimension Ngxndim
Ng=32
Nmg2=(Ng-1)**2
Nmg3=(Ng-1)**3
Ng3=Ng**3
# Target volume fraction of MA
vf=0.52
nbpart=int(vf*Ng3) # number of pixel for phase 2
###########################
# Random microstructure
###########################
# initialization microstructure
micro = int(np.zeros(1))*np.mgrid[:Ng, :Ng, :Ng][0]
#
while (np.sum(micro)<nbpart):
p=random.randint(0, Ng-2)
q=random.randint(0, Ng-2)
l=random.randint(0, Ng-2)
micro[p,q,l]=1
# for the symmetry of the microstructure
if p==0:
micro[Ng-1,q,l]=1
if p==0 and q==0:
micro[Ng-1,Ng-1,l]=1
if p==0 and l==0:
micro[Ng-1,q,Ng-1]=1
if p==0 and q==0 and l==0:
micro[Ng-1,Ng-1,Ng-1]=1
if q==0:
micro[p,Ng-1,l]=1
if q==0 and l==0:
micro[p,Ng-1,Ng-1]=1
if l==0:
micro[p,q,Ng-1]=1
# Restructuration of the microstructure
micro2=micro.reshape(-1,1)
# Actual volume fraction
volfrac=np.sum(micro2)/Ng3
print('volfrac=',volfrac)
## Display of a microstructure
for no in [0,Ng-1]:
fm=no;
x, y = np.mgrid[:Ng, :Ng]
Z = micro[fm,x,y]
# mask out the negative and positive values, respectively
fig = plt.figure(frameon=False)
# Represent the microstructure
pos_neg_clipped = plt.imshow(Z, cmap='binary',
interpolation='none')
x, y, z = np.indices((Ng,Ng,Ng))
# combine the objects into a single boolean array
cube1 = micro[x,y,z]<1
cube2 = micro[x,y,z]>0.5
# combine the objects into a single boolean array
voxels = cube1 | cube2
# set the colors of each object
colors = np.empty(voxels.shape, dtype=object)
colors[cube1] = 'black'
colors[cube2] = 'white'
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.voxels(voxels, facecolors=colors, edgecolor='black', linewidth=0.1)
ax.axis('off')
ax.set_aspect('equal')
plt.savefig('Micro_random.eps', dpi=300)
plt.savefig('Micro_random.png', dpi=300)
plt.show()
##############################################################################
## Microstructure 3
## worm-like microstructure one phase is built by random walk
##############################################################################
# initialization microstructure
micro = int(np.zeros(1))*np.mgrid[:Ng, :Ng, :Ng][0]
p=random.randint(1,Ng-2)
q=random.randint(1,Ng-2)
l=random.randint(1,Ng-2)
micro[p,q,l]=1
while (np.sum(micro)<nbpart):
z=random.randint(1,3)
if z==1:
if p==Ng-2:
p=0
micro[p,q,l]=1
micro[Ng-1,q,l]=1
if q==0:
micro[p,Ng-1,l]=1
micro[Ng-1,Ng-1,l]=1
if l==0:
micro[p,q,Ng-1]=1
micro[Ng-1,q,Ng-1]=1
if q==0 and l==0:
micro[p,Ng-1,Ng-1]=1
micro[Ng-1,Ng-1,Ng-1]=1
else:
p+=1
micro[p,q,l]=1
if q==0:
micro[p,Ng-1,l]=1
if l==0:
micro[p,q,Ng-1]=1
if q==0 and l==0:
micro[p,Ng-1,Ng-1]=1
#
if z==2:
if q==Ng-2:
q=0
micro[p,q,l]=1
micro[p,Ng-1,l]=1
if p==0:
micro[Ng-1,q,l]=1
micro[Ng-1,Ng-1,l]=1
if l==0:
micro[p,q,Ng-1]=1
micro[p,Ng-1,Ng-1]=1
if p==0 and l==0:
micro[Ng-1,q,Ng-1]=1
micro[Ng-1,Ng-1,Ng-1]=1
else:
q+=1
micro[p,q,l]=1
if p==0:
micro[Ng-1,q,l]=1
if l==0:
micro[p,q,Ng-1]=1
if p==0 and l==0:
micro[Ng-1,q,Ng-1]=1
#
if z==3:
if l==Ng-2:
l=0
micro[p,q,l]=1
micro[p,q,Ng-1]=1
if p==0:
micro[Ng-1,q,l]=1
micro[Ng-1,q,Ng-1]=1
if q==0:
micro[p,Ng-1,l]=1
micro[p,Ng-1,Ng-1]=1
if p==0 and q==0:
micro[Ng-1,Ng-1,l]=1
micro[Ng-1,Ng-1,Ng-1]=1
else:
l+=1
micro[p,q,l]=1
if p==0:
micro[Ng-1,q,l]=1
if q==0:
micro[p,Ng-1,l]=1
if q==0 and p==0:
micro[Ng-1,Ng-1,l]=1
# Restructuration of the microstructure
micro2=micro.reshape(-1,1)
# Actual volume fraction
volfrac=np.sum(micro2)/Ng3
print('volfrac=',volfrac)
## Display of a microstructure
for no in [0,Ng-1]:
fm=no;
x, y = np.mgrid[:Ng, :Ng]
Z = micro[fm,x,y]
# mask out the negative and positive values, respectively
fig = plt.figure(frameon=False)
# Represent the microstructure
pos_neg_clipped = plt.imshow(Z, cmap='binary',
interpolation='none')
x, y, z = np.indices((Ng,Ng,Ng))
# combine the objects into a single boolean array
cube1 = micro[x,y,z]<1
cube2 = micro[x,y,z]>0.5
# combine the objects into a single boolean array
voxels = cube1 | cube2
# set the colors of each object
colors = np.empty(voxels.shape, dtype=object)
colors[cube1] = 'black'
colors[cube2] = 'white'
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.voxels(voxels, facecolors=colors, edgecolor='black', linewidth=0.1)
ax.axis('off')
ax.set_aspect('equal')
plt.savefig('Micro_worm.eps', dpi=300)
plt.savefig('Micro_worm.png', dpi=300)
plt.show() 