This python method evaluates the average stresses and strains, given the odb path and step name.

rve_homogenization_post_processing.py

"""
    Suggested to run this script from terminal,
    `abaqus cae noGUI="/path/to/this/file"`
    
"""
from visualization import *
import numpy as np
import os
import sys
import inspect

#
def print_output(msg,terminal=True):
    if terminal:
        print >> sys.__stdout__, msg
    else:
        print msg
#
#
def array_pretty_printing(arr, name=""):
    if type(arr)==np.ndarray:
        print_output("{}::".format(name))
        if arr.ndim==2:
            for a_row in arr:
                print_output("    ".join("{0:^10.6E}".format(k) for k in a_row))
        elif arr.ndim==1:
            print_output("    ".join("{0:^10.6E}".format(k) for k in arr))
    elif type(arr) in (int, float):
        print_output("{0}:: {1}".format(name, arr))
#
#
def get_volume_averaged_values(odb_path, step_name, print_res=False):
    #
    odb = openOdb(path=odb_path, readOnly=True)
    last_frame = odb.steps[step_name].frames[-1]
    #
    ip_vol = last_frame.fieldOutputs["IVOL"].getSubset(position=INTEGRATION_POINT).values    
    ip_stresses = last_frame.fieldOutputs["S"].getSubset(position=INTEGRATION_POINT).values
    ip_strains = last_frame.fieldOutputs["E"].getSubset(position=INTEGRATION_POINT).values
    #
    cum_wtd_stresses = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    cum_wtd_strains = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    total_strain_energy = 0.0
    cum_vol = 0.0
    for i in range(len(ip_vol)):
        ith_ip_vol = ip_vol[i].data
        cum_wtd_stresses += (ip_stresses[i].data * ith_ip_vol)
        cum_wtd_strains += (ip_strains[i].data * ith_ip_vol)
        total_strain_energy += (ip_strains[i].data * ip_stresses[i].data * ith_ip_vol * 0.5)
        cum_vol += ith_ip_vol
        #
    vol_avg_stress = cum_wtd_stresses/cum_vol
    vol_avg_strain = cum_wtd_strains/cum_vol
    total_strain_energy = sum(total_strain_energy)
    if print_res:
        print_output("\nVolume averaged quantites where the total volume is {0}".format(cum_vol))
        array_pretty_printing(vol_avg_stress, "Averaged Stress")
        array_pretty_printing(vol_avg_strain, "Averaged Strain")
        array_pretty_printing(total_strain_energy, "Microscopic Strain Energy")
    return vol_avg_stress, vol_avg_strain, total_strain_energy
#
#
def get_voigtSEvec_from_abaqusSEvec(arr, ax_map=None):
    """
    Get the stress(S)/Strain(E) vector in Voigt notation, 
    given the corresponding vector generated by ABAQUS.
    Result depends on the position of material axes(1,2,3) relative to the abaqus global CSYS(x,y,z).
    ABAQUS notation of stress/ strain vector => 11, 22, 33, 23, 31, 12
    ABAQUS notation of stress/ strain vector => xx, yy, zz, xy, xz, yz
    """
    assert len(arr)==6
    if ax_map==None:
        ax_map={"x":2, "y":3, "z":1}
   
    for (csys_ax, mat_ax) in ax_map.items():
        if mat_ax==1:
            one=csys_ax
        elif mat_ax==2:
            two=csys_ax
        elif mat_ax==3:
            three=csys_ax
        else:
            raise Warning("material axis can be 1, 2 or 3")
    #
    voigt_notation = [one+one, two+two, three+three, two+three, three+one, one+two]
    abaqus_notation = ["xx", "yy", "zz", "xy", "xz", "yz"]
    notation_mapping = []
    for i in voigt_notation:
        if i in abaqus_notation:
            notation_mapping.append(abaqus_notation.index(i))
        elif i[::-1] in abaqus_notation:
            notation_mapping.append(abaqus_notation.index(i[::-1]))
        else:
            raise Warning("Component {0} is not in abaqus notation {1}".format(i, abaqus_notation))
    print(notation_mapping)
    #
    return np.array([arr[i] for i in notation_mapping])
#
#
def get_Cijkl_col(Savg_abq, Eavg_abq, col_num):
    Savg_Voigt = (get_voigtSEvec_from_abaqusSEvec(Savg_abq)).reshape(-1,1)
    Eavg_Voigt = (get_voigtSEvec_from_abaqusSEvec(Eavg_abq)).reshape(1,-1)
    return (Savg_Voigt/Eavg_Voigt)[:, col_num-1]
#
#
def get_engineering_const_from_compliance_matrix(s_matrix):
    assert s_matrix.shape==(6,6)
    s11, s12, s13, s14, s15, s16 = s_matrix[0, :]
    s22, s23, s24, s25, s26 = s_matrix[1, 1:]
    s33, s34, s35, s36 = s_matrix[2, 2:]
    s44, s45, s46 = s_matrix[3, 3:]
    s55, s56 = s_matrix[4, 4:]
    s66 = s_matrix[5, 5]
    results = {}
    results.update({
    "E11"   : 1.0/s11,
    "E22"   : 1.0/s22,
    "E33"   : 1.0/s33,
    "nu_12" : -s12/s11,
    "nu_21" : -s12/s22,
    "nu_23" : -s23/s22,
    "nu_32" : -s23/s33,
    "nu_13" : -s13/s11,
    "nu_31" : -s13/s33,
    "G23"   : 1.0/s44,
    "G13"   : 1.0/s55,
    "G12"   : 1.0/s66,
    })
    return results
#
#
def main(odb_paths, step_id = "STEP-1", print_vol_avg_res=False, csys_voigt_map=None, msys_voigt_map=None):
    if csys_voigt_map is None:
        csys_voigt_map = {"XX": 2, "YY": 3, "ZZ": 1, "YZ": 5, "ZX": 6, "XY": 4,}
    if msys_voigt_map is None:
        msys_voigt_map = {"22": 2, "33": 3, "11": 1, "31": 5, "12": 6, "23": 4,}
    #
    C_ijkl_columns = {}
    for (a_odb_id, a_odb_path) in odb_paths.items():
        print_output("Working on {0} case".format(a_odb_id))    
        column_index = csys_voigt_map[a_odb_id[1:]]
        a_avg_stress, a_avg_strain, a_Senergy = get_volume_averaged_values(a_odb_path, step_id, print_vol_avg_res)
        C_ijkl_columns[column_index] = get_Cijkl_col(a_avg_stress, a_avg_strain, column_index)
    #
    eff_C_ijkl = np.column_stack(tuple(C_ijkl_columns[i] for i in range(1, 7)))
    eff_S_ijkl = np.linalg.inv(eff_C_ijkl)
    return eff_C_ijkl, eff_S_ijkl


# ===================================
#               MAIN 
# ===================================
#
# getting the path of the file
CFD = os.path.dirname(inspect.getfile(lambda: None))
t0 = time.time()
odbs_dir = CFD # replace with directory of choice
results_dir = CFD # replace with directory of choice
Odb_Paths = {
    "EXX": os.path.join(odbs_dir, "E22.odb"), # Transverse Tensile
    "EYY": os.path.join(odbs_dir, "E33.odb"), # Transverse Tensile
    "EZZ": os.path.join(odbs_dir, "E11.odb"), # Axial Tensile
    "GYZ": os.path.join(odbs_dir, "G31.odb"), # In-plane Shear
    "GZX": os.path.join(odbs_dir, "G12.odb"), # In-plane Shear
    "GXY": os.path.join(odbs_dir, "G23.odb"), # Transverse Shear
}
#
C_ijkl, S_ijkl = main(Odb_Paths, "STEP-1", True)
#
array_pretty_printing(C_ijkl, "stiffness_matrix")
array_pretty_printing(S_ijkl, "compliance_matrix")
#
CS_ijkl = np.dstack((C_ijkl, S_ijkl))
np.save(os.path.join(results_dir, "{0}_CSijkl.npy".format(a_case_id)), CS_ijkl)
#
eng_const = get_engineering_const_from_compliance_matrix(S_ijkl)
#
print_output("Engineering constants : ")
for (k,v) in eng_const.items():
    print_output("{0} :: {1}".format(k,v))       
# Export to a text file
with open(os.path.join(results_dir, "{0}_results.txt"), "w") as file1:
    for (k,v) in eng_const.items():
        file1.write("{0} :: {1}\n".format(k,v))
#
end_card = "Finished the case {0} in {1:4.3f} seconds".format(a_case_id, time.time()-t0)
print_output("\n\n{0}\n{1}\n{0}\n".format(len(end_card)*"~", end_card))