Attaching my code which I can't manage to upload :

# Open Odb session

myOdb = jobName + '.odb'

odb_Mehcanical = session.openOdb(name = myOdb)

# Definning region selections for stress, strain, and IVOL data extraction

myAssembly = odb_Mehcanical.rootAssembly

myRegion = myAssembly.elementSets[' ALL ELEMENTS']

# Number of steps in the job

numOfSteps = len(odb_Mehcanical.steps)

# Initiating Stiffness matrix

myStiffMatrix = [ [ 0 for i in range(6) ] for j in range(6) ]

# RVE's volume

totalVolume = xlength*ylength*zlength

# Iterating through steps

for i in range(numOfSteps):

stepIndex = i # i ranges from [0 1 2 3 4 5], total of 6 steps.

print 'Step Index = ' + str(stepIndex)

print ' '

# Current step

myStep = odb_Mehcanical.steps.values()[stepIndex]

# Path for field stresses calculated at integration point

# frames[-1] select the last deformed frame.

myStress = myStep.frames[-1].fieldOutputs['S'].getSubset(position=INTEGRATION_POINT, region = myRegion)

# Path for field strain calculated at integration point

myStrain = myStep.frames[-1].fieldOutputs['E'].getSubset(position=INTEGRATION_POINT, region = myRegion)

# Path for field Volume cahnge of integration point

myVolume = myStep.frames[-1].fieldOutputs['IVOL'].getSubset(position=INTEGRATION_POINT, region = myRegion)

# Number of elements in the model

myModel = p # p was defined earlier as the model's meshed assembled part.

len1 = len(myModel.elements[:])

# Checking how many integration points need to iterate

if (p.elements[0].type == C3D10):

total_Length = len1*4 # for each integration point in C3D10 element.

else:

total_Length = len1 # in case of linear element (C3D4)

# Indexs, reference.

# S11 = 0 # S22 = 1 # S33 = 2 # S12 = 3 # S13 = 4 # S23 = 5

# E11 = 0 # E22 = 1 # E33 = 2 # E12 = 3 # E13 = 4 # E23 = 5

# Initiating a

effectiveStrain = 0

# Iterating through stress components for each displacement (i).

# Requires less iterations each time as the tensor is symmetric.

for j in range(i, 6):

myStress_temp = 0

myStrain_temp = 0

myStrainCalc = 0

# iterating through integration points

for k in range(total_Length):

# Stress

myStress_temp = myStress.values[k].data[j] # k for integration point; j for stress component

if (j == i):

# Strain

myStrain_temp = myStrain.values[k].data[j] # k for integration point; j for strain component

# Volume

myVol = myVolume.values[k].data # k for integration point;

# Initial Stiffness Tensor(matrix) precalculations

# This line Sums the product of stress(per)element times volumeChange(per)element

myStiffMatrix[i][j] += myStress_temp*myVol

# Strain

# This line Sums the product of Strain/element times volumeChange/element

myStrainCalc += myStrain_temp*myVol

# Effective Strain, should be calculated only in the direction of the applied displacement.

if (j == i):

effectiveStrain = (1/totalVolume) * myStrainCalc

print 'MyStrainCalc = ' + str(myStrainCalc)

print ' '

print 'Stiffness Tensor column equals: ' + str(j)

print 'Step number is ' + str(i)

print myStep

print 'Effective Strain is ' + str(effectiveStrain)

print ' '

# Calculating the Stiffness Tensor through homogenization formula

myStiffMatrix[i][j] = (1/totalVolume) * myStiffMatrix[i][j]

myStiffMatrix[i][j] = myStiffMatrix[i][j] / (effectiveStrain)

# Due to the tensor's symmetry:

if (i != j) and (i < 5):

myStiffMatrix[j][i] = myStiffMatrix[i][j]

print myStiffMatrix

# Changing data type to allow matrix type calculations

stiff_mat = mat(myStiffMatrix)

# Calculating the Compliance Tensor

comp = stiff_mat**(-1)

print comp

myStiffMatrix = mat(myStiffMatrix)

# Averaging the first 3 diagonal components of the Compliance tensor.

effElasticModulus = (1/comp[0,0] + 1/comp[1,1] + 1/comp[2,2]) / 3

# Averaging the poisson's ratio components from the Compliance tensor.

effpRatio = (-1) * effElasticModulus * ((comp[0,1] + comp[0,2] + comp[1,2]) / 3)

myEffectiveProperties = [effElasticModulus, effpRatio]

print ' '

print 'Effective Elastic Modulus = ' + str(effElasticModulus)

print 'Effective Poisson''s ratio = ' + str(effpRatio)

Dear Gilad,

The main problem is related to macro strain calculations.

When there are voids or cracks inside the RVE, the macro strain components shouldn't be calculated by averaging the local strains over the volume of the RVE. Because, The strain components inside the voids are unknown and there is displacement discontinuity near the edges of the voids. In contrast, the stress in the voids = 0. And, normally, the value we obtain using this technique is smaller than the real value and by using this technique, higher number voids mean lower macro strain. Therefore, the effective stiffness increases with the increase of voids. So, in your case, the correct value of macro strain can be calculated from the average displacements of the RVE's corners (Macro strain in x for example = delta Ux/Lx). Here, Ux is the change in displacement of the corners and Lx is the initial length in the x-direction.

In the absence of voids, both methods converge to the same value (calculating the macro strain using the average method or from RVE's corners displacement). You can check that.

Good luck.

Mahmoud

Dear Gilad,

You can check the homogenization method and the FE model used in this article for RVEs that contain voids or zero stiffness elements:

"Efects of Defects on Efective Material Properties of Triaxial Braided

Textile Composite"

https://link.springer.com/article/10.1007/s42405-019-00244-8

Dear Mahmoud Ali ,

Thank you very much for your answer.

Your answer makes sense to me.

Probing the displacements at the corner nodes(on the displaced face) will just give me the prescribed initial displacement. That is, if I have a 400*400*400 RVE, with prescribed displacement of 1 unit along the x axis, while maintaining periodicity, and preventing rigid body motion by fully constraining one corner node, I should get a macro strain that is equal to 1/400 = 0.0025. In that case, why should I probe the displacements and not just use a fixed value?

Thanks again,

- Gilad

Dear Gilad Nave

Under the BCs you mentioned, you can use the prescribed displacement to calculate the macro strain. Then you solve the boundary value problem to get the macro stress. The Applied value is equal to the average value. Keep this in mind, if the applied value is either stress or strain (displacement).

If you need more info, please let me know.

mahmoud

Hello, young researchers

I see that you are in the stage I was in 1995 when I began to deal with problems of homogenization. In 2005 I did a monograph on the subject , where You can found answers to all questions You had mention (and even more) .It might be hardly available, so I send You pdf. Please , do the refence Gilad Nave Mahmoud Ali Mohsen Barzegar

pdf Gilad Nave Mahmoud Ali Mohsen Barzegar

Dear Aleksander Urbanski,

Thank you for sharing your monograph. It looks very informative and fun to read, and the illustrations are great!

I appreciate your help.

Sincerely,

Gilad