Capstone Logbook

W2

Capstone B W2 logbook

• Capstone Project logbook (@2025-W31-D7 23:54) (todo!())

This Week

Summary

• 2D lap joint -> 3D lap joint (still isotropic adherends)

Monday

Meeting w supervisor.

Thursday

Creating a 3D lap joint model (currently with isotropic adherends). Got better at seeding the mesh edges. Initially the mesh was way too fine.

I later fixed this to have less elements.

Here it also shows where the adhesive has yielded. This is with an eccentric load downwards and to the left. Next step is to make the adherends composite materials.

What I learnt

• Don’t create a mesh with over a million elements.
• Modelling the 3D case is very similar it’s just that you extrude the original sketch
• You can easily seed the edges by setting the view to partially transparent and dragging, as opposed to trying to click them
• It’s really hard to set the correct force in an Abaqus/Explicit model without damage modelling because if it’s too much end exceeds the adhesives ultimate strength it wants to just extend forever.

Todo

• Composite adherends
• Discuss with team what the loading shape should be

Link to original

W3

Capstone B W3 logbook

• Capstone Project logbook (@2025-W32-D7 23:54)
  

  Todo

  • Composite adherends
  • Discuss with team what the loading shape should be

  Link to original

This week

Monday

Meeting w/ Adrian. He said we should work towards having a theme of repair. This makes me think that moving towards a plug repair would be good. Ethan will have to change the numbers around for some of the single axis loading hand calcs but that should be fine.

Wednesday

Got composite adherend for the lap joint model. This one is a different model because it uses shells for the composites instead of the 3d elements. The composite follows the Hashin Criterion and shows when the fiber or matrix fails either from compression or tension (as shown below, where it’s showing the fibre tensile failure criterion). The composite doesn’t have damage evolution but that’s okay. The adhesive is modelled with cohesive elements.

Also, the force comes from the displacement of the end instead of applying a direct edge load. This is way better because it doesn’t lead to infinite stretching. In the photo above, the adhesive isn’t deleting even though I applied a damage criterion based on max stress, which is addressed later.

The adhesive damage and deletion is also working. I realised that the issue was that I hadn’t applied a mesh stacking direction nor a material orientation to the cohesive elements.

Thursday

I took the learning from the lap joint model and used it to make a scarf joint. First I modelled the damaged surface (aluminum) with a shell:

The part in the center is weakened by making the young’s modulus only 10% of what it was before. I will have to confirm with Jalisha to make sure this is the best way to do it. Then I made an adhesive layer (cohesive elements) and tied that to the surface. The cohesive elements need definitions for traction based separation. I tied a composite shell to the adhesive layer underneath to act as the patch.
Here is the right edge being pulled down:

And the right edge being pulled out:

It makes sense that the adhesive fails in the middle as that is where the material is weakest.

Right now the material properties just use magic numbers. Ethan will be filling in the material data eventually but I can keep working with the values I have now and change them later.

Friday

Created a document that we can put the material properties in. It’s just looking at what Abaqus needs and representing it in a way that allows us to easily reference it later.

Working on the completion plan diagram. I’ve been making it in Microsoft Visio because it has teams integration so it’s the easiest to work with.

I’ve discussed the plan with Ethan and he seems happy with it for the most part it just comes down to getting the details from the scarf repair side, as well as finalising who will do what in the second half of the semester.

What I learnt

• One way to define adhesives is to use cohesive elements that have a traction definition. You must define the mesh stack direction and material orientation. Then you have to define the elasticity in terms of traction. For damage, you’ll also have to define damage criterion. For deletion, you’ll have to specify damage progression as a sub-property.
• When you’re running jobs so many times, just make a coarser mesh initially and then worry about making it finer later.
• Microsoft Visio kinda sucks. The arrows are very unintuitive. I would’ve preferred excalidraw but no Teams integration.

Todo

• Confer with Jalisha to get the actual dimensions for the damage
• Try to run the jobs with python

Link to original

W4

Capstone B W4 logbook

• Capstone Project logbook (@2025-W33-D7 23:54) (todo!())

Todo

• Confer with Jalisha to get the actual dimensions for the damage
• Try to run the jobs with python

Link to original

This Week

Monday

Meeting with Adrian. He’s happy with the fact that we’re assigning/micromanaging tasks more.
Jalisha and I finalised that the crack length should be 3mm

Source
We have yet to finalise material properties but that can come later as it’s just a matter of changing numbers in Abaqus.
We have also decided to try to do as much of the report as possible before the poster because then making the poster is just a matter of summarising existing information :3

Tuesday

We have confirmed that the modelled surface will be 5 x 50 x 15 (depth) mm for the sake of having consistent analyses.

Wednesday

I added some information on the modelling of the patch repair to the completion plan introduction.

Thursday

Finalising completion timeline by adding in the tasks that Jalisha, Sam and Mitch will do for the scarf repair.

Sunday

Proofreading completion plan and adding patch repair info to the Completed Tasks part.

I can do a displacement load based off a normal co-ordinate system but not for a cylindrical system (which I wanted to use so that I can apply a torsional displacement). It gives an error for having elements on the rotational axis itself right now so maybe a rotation would have to be define through a script or some other method (like just rotating the corner nodes).

Submitted completion plan!

Todo

• Touch base w/ Ethan on material properties
• Get the jobs running via Python scripts
• Find out how to apply a torsional displacement BC

Link to original

W5

Capstone B W5 logbook

Todo

• Touch base w/ Ethan on material properties
• Get the jobs running via Python scripts
• Find out how to apply a torsional displacement BC

Link to original

This week

Monday

Meeting with Adrian. Only caught the end of it because I slept in. Overall he’s happy with the report (besides formatting)

Wednesday

Been working on python scripts for the model. Right now I’m getting a bit confused over what I can do in Abaqus/CAE vs the Abaqus kernel, but overall I have a VSCode environment working using the Abaqus API, which gives good typing information https://pypi.org/project/pyabaqus/. So far I’ve just been learning the environment and tinkering with existing scripts. The script below is modified to give more field values such as the plastic strain as well.

Max Mises Stress Script

VonMisesPy

 
import math
from odbAccess import openOdb, Odb
from abaqus.Odb.FieldValue import FieldValue
from sys import argv, exit
 
 
def rightTrim(input, suffix):
    if input.find(suffix) == -1:
        input = input + suffix
    return input
 
 
from dataclasses import dataclass
from typing import Any, Callable, Generic, Optional, Tuple, TypeVar
 
 
@dataclass
class StepFrame:
    stepName: str
    step: Optional[int]
    frame: int
 
    def __str__(self) -> str:
        return f'"{self.stepName}" (S{self.step}) F{self.frame}'
 
    def __hash__(self):
        return hash((self.stepName, self.frame))
 
 
@dataclass
class DataPoint:
    fieldValue: FieldValue
    element: Optional[int]
    node: Optional[int]
 
    def __str__(self) -> str:
        return f"{self.fieldValue} in element #{self.element}"
 
 
FieldComparisonCallable = Callable[[FieldValue], float]
 
 
def getMaxFieldValue(
    odb: Odb,
    elsetName: Optional[str],
    fieldOutputName: str,
    fieldComparisonKey: FieldComparisonCallable,
) -> dict[StepFrame, Optional[DataPoint]]:
    """Print max mises location and value given odbName
    and elset(optional)
    """
    assembly = odb.rootAssembly
 
    # process the element subset
    if elsetName:
        if elsetName in assembly.elementSets:
            elemset = assembly.elementSets[elsetName]
        else:
            raise ValueError(
                "An assembly level elset named %s does"
                "not exist in the output database %s" % (elsetName, odbName)
            )
    else:
        elemset = None
 
    isFieldPresent = False
    maxFieldKey = -float("inf")
    maxFieldValues = {}
 
    for stepName, step in odb.steps.items():
        print("Processing Step:", step.number)
        for frame in step.frames:
            stepFrame = StepFrame(stepName, step.number, frame.incrementNumber)
            frameMaxValue: Optional[DataPoint] = None
            frameFields = frame.fieldOutputs
            if fieldOutputName in frameFields:
                isFieldPresent = 1
                fieldOutput = frameFields[fieldOutputName]
                if elemset:
                    fieldOutput = fieldOutput.getSubset(region=elemset)
                fieldValueArray = fieldOutput.values
                if fieldValueArray is None:
                    continue
                for fieldValue in fieldValueArray:
                    fieldKey = fieldComparisonKey(fieldValue)
                    if fieldKey > maxFieldKey:
                        maxFieldKey = fieldKey
                        frameMaxValue = DataPoint(
                            fieldValue=fieldValue,
                            element=fieldValue.elementLabel,
                            node=fieldValue.nodeLabel,
                        )
            maxFieldValues[stepFrame] = frameMaxValue
        # end for frame
 
    if not isFieldPresent:
        raise ValueError("No field present in element set")
 
    return maxFieldValues
 
 
if __name__ == "__main__":
    odbName = None
    elsetName = None
    argList = argv
    argc = len(argList)
    i = 0
    while i < argc:
        if argList[i][:2] == "-o":
            i += 1
            name = argList[i]
            odbName = rightTrim(name, ".odb")
        elif argList[i][:2] == "-e":
            i += 1
            elsetName = argList[i]
        elif argList[i][:2] == "-h":
            print(__doc__)
            exit(0)
        i += 1
    if not (odbName):
        print(" **ERROR** output database name is not provided")
        print(__doc__)
        exit(1)
 
    odb = openOdb(odbName)
 
    fieldComparisonKey: FieldComparisonCallable = (
        lambda fv: fv.mises if fv.mises is not None else math.nan
    )
 
    try:
        maxMisesValues = getMaxFieldValue(odb, elsetName, "S", fieldComparisonKey)
        for stepFrame, maxMisesPoint in maxMisesValues.items():
            value = None
            if maxMisesPoint is not None:
                value = fieldComparisonKey(maxMisesPoint.fieldValue)
            print(f"at {stepFrame}: Max Von Mises Stress = {value}")
    except Exception as e:
        print(e)
 
    fieldComparisonKey: FieldComparisonCallable = lambda fv: fv.data[0]
 
    try:
        maxMisesValues = getMaxFieldValue(odb, elsetName, "PE", fieldComparisonKey)
        for stepFrame, maxMisesPoint in maxMisesValues.items():
            value = None
            if maxMisesPoint is not None:
                value = maxMisesPoint.fieldValue.data
            print(f"at {stepFrame}: Max Plastic Strain = {value}")
    except Exception as e:
        print(e)
 
    odb.close()
 

Link to original

Todo

• Touch base w/ Ethan on material properties
• Get the jobs running via Python scripts
• Find out how to apply a torsional displacement BC

Link to original

(Literally the same as last week. Have been focusing on structures studio)

Link to original

W6

Capstone B W6 logbook

Todo

• Touch base w/ Ethan on material properties
• Get the jobs running via Python scripts
• Find out how to apply a torsional displacement BC

Link to original

This week

Monday

Meeting with Adrian. Not much to note last week was mostly about completing the structures assignment.

Thursday

I found out how to define a boundary condition requiring torsion using analytical fields. The field defines the distribution, and that gets scaled by the constraint you pick. So, first we start with the rotation matrices:

$$\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)=\left(\begin{array}{cc}\cos (\theta )& -\sin (\theta )\\ \sin (\theta )& \cos (\theta )\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)-\left(\begin{array}{c}x\\ y\end{array}\right)$$ $$\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)=\left(\begin{array}{cc}\cos \theta -1& -\sin \theta \\ \sin \theta & \cos \theta -1\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$$

In Abaqus we define the two equations separately, one for horizontal and one for vertical

$$\begin{array}{rl}\Delta x& =x(\cos \theta -1)-y\sin \theta \\ \Delta y& =x(\sin \theta )+y(\cos \theta -1)\end{array}$$

Now we define the fields while picking a $\theta$ to use. A key difference is that due to our co-ordinate system being ZY the equations look a little different


And then we apply these analytical fields to the stretch nodes:

There was more boundary conditions added, such as setting the node rotation on torsion to $\theta$ as well as fixing the other DOF (U1, UR2, UR3).

There was some difficulty, namely with making sure co-ordinate systems and displacements were correct. One funny thing I couldn’t understand was why my whole model was rotating, and that was because I had temporarily disabled the wall BC. Unfortunately wasted half an hour on that, thinking it was a buckling issue (as I was getting negative eigenvalue errors).

Eventually, the torsion worked out alright:

Currently exact adhesive values are still being found by other group members

Since this is a displacement load, it’s not the actual torsional moment being applied. That can be done by looking at the reaction forces and moments. In fact, you can define a derived field output based on the reaction force RF and reaction moments RM at a node.

Friday

Remodelling

For Friday, the main focus was to remodel the mesh to actually be more indicative of a patch repair. The shape of the surface was modified to be a square and a 10x20 patch was created to sit over the 3mm crack.

Adaptive meshing

In order to remove errors in stress calculation, I remeshed the surface and patch to make the patch corner stresses more realistic, as well as better model the stress around the crack. The adhesive could not be remeshed as Abaqus only supports remeshing tri and quad elements, but the patch remeshing and stresses gave insight into how the adhesive mesh would experience stresses.

^ Patch before and after remeshing

During the remeshing process, it was found that the sharp corners of the initial mesh were carrying unrealistic stresses and so these corners were rounded and the remeshing was done again. Another change was changing the loading from twisting (which would be non-symmetrical) to pure tension. This would still have similar stress concentrations to the torsional case but would make sure the remeshing would be relatively symmetrical.

What I learnt

• You can create basically any positional constraints by using an analytical field.
• You can operate on field values to create new field values (RF, U, and RM to give total rotational moment)
• You can’t really use an adaptivity process to decrease the number of elements, only increase.]

Todo

• In Python:
  • Run a single job
  • Put at least one result we want in a different object database so we can compare stuff from multiple runs
• Follow up on material values w/ Ethan

Link to original

W7

Capstone B W7 logbook

Todo

• In Python:
  • Run a single job
  • Put at least one result we want in a different object database so we can compare stuff from multiple runs
• Follow up on material values w/ Ethan

Link to original

This Week

Monday

I inputted the new material values for the aluminium, adhesive, and patch. The adhesive is now using significantly stronger values, which makes sense as an adhesive needs to be strong.

After doing this, I’m focusing on making scripts that can help me run tests and extract information. A significant portion of this time was just getting a simple working script for extracting odb information easily. I am making many mistakes but I am also learning many things, which is a plus.

After finally setting up a good way to get position (COORD), displacement (U/UR), and force (RF/RM) values for the wall and stretch positions, I can now use the nodal values at the stretch location to get specific useful values:

# for each frame:
tension = -sum(val.data[0] for val in fieldOutputs["RF"])
stretch = statistics.mean(val.data[0] for val in fieldOutputs["U"])
elongation = stretch / GEOMETRY["Length"]

From this I can then start getting values of crack stress vs applied tensile stress:

What I learnt

• You can also specify a COORD field output to get the integration point original positions
  
• With FEA, some field outputs are at nodes, whereas some are at elements. The sets for these are different. You cannot get displacements (U) on element sets; you must use node sets.
• You can use the Abaqus/CAE kernel line to easily sandbox things before putting them in code

>>> step.frames[0].fieldOutputs["COORD"].values[0].data
array([-5. , -0.075 , 9.799999], dtype=float32)
>>> step.frames[0].fieldOutputs["COORD"].getSubset(wallSet)
TypeError: arg1; found OdbSet, expecting UNDEFINED_POSITION, NODAL, ELEMENT_NODAL, INTEGRATION_POINT, ELEMENT_FACE_INTEGRATION_POINT, WHOLE_ELEMENT, WHOLE_REGION, WHOLE_MODEL, CENTROID, SURFACE_INTEGRATION_POINT, SURFACE_NODAL, ELEMENT_FACE or WHOLE_PART_INSTANCE
>>> step.frames[0].fieldOutputs["COORD"].getSubset()
openOdb(r'C:/Users/.../Twist.odb').steps['Step-1'].frames[0].fieldOutputs['COORD'].getSubset()
>>> step.frames[0].fieldOutputs["COORD"].getSubset(region=wallSet)
openOdb(r'C:/Users/.../Twist.odb').steps['Step-1'].frames[0].fieldOutputs['COORD'].getSubset(region=openOdb(r'C:/Users/.../Twist.odb').rootAssembly.elementSets['WALL'])
>>> wallOutputs = step.frames[0].fieldOutputs["COORD"].getSubset(region=wallSet)
>>> wallOutputs.values
openOdb(r'C:/Users/.../Twist.odb').steps['Step-1'].frames[0].fieldOutputs['COORD'].getSubset(region=openOdb(r'C:/Users/.../Twist.odb').rootAssembly.elementSets['WALL']).values
>>> wallOutputs.values[0].data
array([-58.564137, 0. , -47.995377], dtype=float32)
>>> wallOutputs = step.frames[0].fieldOutputs["U"].getSubset(region=wallSet)
>>> wallOutputs.values[0].data
IndexError: Sequence index out of range

C:\Users\S3943498\...\Lapjoint3d>abaqus python abaqus_plugins/loadinginfo.py -odb Job-Patch.odb
Processing Step  1
"Step-1" F0 LoadingCase(tension=0.0, stretch=0.0, elongation=0.0, crackStress=0.0)
"Step-1" F1 LoadingCase(tension=14032.019073486328, stretch=0.1, elongation=0.0016666666915019354, crackStress=70.94763946533203)
"Step-1" F2 LoadingCase(tension=28032.499877929688, stretch=0.2, elongation=0.0033333333830038708, crackStress=141.6964569091797)
"Step-1" F3 LoadingCase(tension=42001.510009765625, stretch=0.3, elongation=0.005000000198682149, crackStress=212.3720245361328)
"Step-1" F4 LoadingCase(tension=55933.23010253906, stretch=0.4, elongation=0.0066666667660077415, crackStress=283.9537048339844)
"Step-1" F5 LoadingCase(tension=69818.291015625, stretch=0.5, elongation=0.008333333333333333, crackStress=358.57122802734375)
"Step-1" F6 LoadingCase(tension=83657.84155273438, stretch=0.6, elongation=0.010000000397364298, crackStress=436.7843933105469)
"Step-1" F7 LoadingCase(tension=97455.65112304688, stretch=0.7, elongation=0.011666666467984517, crackStress=518.5008544921875)
"Step-1" F8 LoadingCase(tension=111215.55419921875, stretch=0.8, elongation=0.013333333532015483, crackStress=603.2343139648438)
"Step-1" F9 LoadingCase(tension=124942.49658203125, stretch=0.9, elongation=0.014999999602635702, crackStress=689.9479370117188)
"Step-1" F10 LoadingCase(tension=138631.0654296875, stretch=1.0, elongation=0.016666666666666666, crackStress=778.7988891601562)
"Step-1" F11 LoadingCase(tension=152280.49853515625, stretch=1.1, elongation=0.01833333373069763, crackStress=871.8869018554688)
"Step-1" F12 LoadingCase(tension=165906.837890625, stretch=1.2, elongation=0.020000000794728596, crackStress=964.7281494140625)
"Step-1" F13 LoadingCase(tension=179509.77685546875, stretch=1.3, elongation=0.02166666587193807, crackStress=1056.3543701171875)
"Step-1" F14 LoadingCase(tension=193085.2802734375, stretch=1.4, elongation=0.023333332935969033, crackStress=1148.4420166015625)
"Step-1" F15 LoadingCase(tension=206632.98095703125, stretch=1.5, elongation=0.025, crackStress=1241.2662353515625)
"Step-1" F16 LoadingCase(tension=220152.28076171875, stretch=1.6, elongation=0.026666667064030966, crackStress=1334.7576904296875)
"Step-1" F17 LoadingCase(tension=233644.70556640625, stretch=1.7, elongation=0.02833333412806193, crackStress=1425.341796875)
"Step-1" F18 LoadingCase(tension=247105.5390625, stretch=1.8, elongation=0.029999999205271403, crackStress=1518.4952392578125)
"Step-1" F19 LoadingCase(tension=260536.919921875, stretch=1.9, elongation=0.03166666626930237, crackStress=1609.6663818359375)
"Step-1" F20 LoadingCase(tension=273938.232421875, stretch=2.0, elongation=0.03333333333333333, crackStress=1701.3314208984375)
Data written to Job-Patch.odb-data.csv

Link to original

W13

Capstone B W13 logbook

Small gap between updates. Many things were done, such as:

• Developing Python code to create the input files.
• Developing Python code that could run all the jobs in a batched way
• Developing python code to read the ODBs and get key information
• Use that information in MATLAB to then get more data.
• Develop Python code to get the contour images for me (Why do something repetitive for 2 hours if I can just spend a mere 6 hours automating it)
• Actual analysis of results

I ended up creating 434 jobs, which took over ~10 hours to run

I’m simply going to embed some files here and some in the pdf I’ll submit, such as the matlab code.

PatchAnalysisPy

from odbAccess import *
from abaqus import *
import os
# from loadinginfo import createLoadingData
from typing import Any
import numpy as np
import datetime
import traceback
 
mdb: Mdb = openMdb("../patch.cae")
model = mdb.models["Model-1"]
# set loading
 
 
def jobWithLoading(mdb, model, dx, dz, dr, with_patch):
    # suppressing/enabling
    def setFeatureStatus(parent, objects, enable):
        if enable:
            return parent.resumeFeatures(objects)
        else:
            return parent.suppressFeatures(objects)
    def setStatus(parent, objects, enable):
        for object in objects:
            if enable: parent[object].resume()
            else: parent[object].suppress()
    setFeatureStatus(model.rootAssembly, ("Adhesive-Layer-1", "Patch-1"), with_patch)
    setStatus(model.constraints, ("Adhesive-Patch-1", "Surface-Adhesive"), with_patch)
    setStatus(model.fieldOutputRequests, ("Adhesive", "Composite"), with_patch)
    model.ExpressionField("RotationY", "Z*sin(%f) + Y*(cos(%f)-1)" % (dr, dr))
    model.ExpressionField(
        "RotationZ", "Z * (cos(%f)-1) + Y*sin(%f) + %f" % (dr, dr, dz)
    )
    model.boundaryConditions["StretchR"].setValues(u1=dx, ur1=-dr) # must be negative rotation
    patchText = "" if with_patch else "less"
    name = ("Patch%s%+04.2fX%+04.2fZ%+05.2fR" % (patchText, dx, dz, dr * 180 / pi)).replace(".", "_")
    print(name)
    return mdb.Job(name, model), name
 
 
def writeInputsInRange(mdb, model, dxs, dzs, drs, patch_statuses):
    jobs = {}
    dxmax = max(1e-6,max(dxs))
    dzmax = max(1e-6,max(dzs))
    drmax = max(1e-6,max(drs))
    print(dxmax, dzmax, drmax)
    for with_patch in patch_statuses:
        for dx in dxs:
            for dz in dzs:
                for dr in drs:
                    if dx != dxmax and dz != dzmax and dr != drmax:
                        continue
                    job, name = jobWithLoading(mdb, model, dx, dz, dr, with_patch)
                    job.writeInput()
                    job = mdb.JobFromInputFile(name, name + ".inp")
                    jobs[name] = (job, dx, dz, dr, with_patch)
                    print(jobs[name])
    return jobs
 
 
 
def fieldComponent(fieldOutputs, fieldOutput, region, getComponent):
    return np.array(np.nan) if not region else np.array([
        getComponent(value)
        for value in fieldOutputs[fieldOutput].getSubset(region=region).values
    ])
 
 
def volumeMean(fieldOutputs, fieldOutput, region, getComponent):
    fo = fieldComponent(fieldOutputs, fieldOutput, region, getComponent)
    return np.mean(fo)
    # print("FO:" + str(fo.shape))
    # vol = fieldComponent(fieldOutputs, "EVOL", region, lambda v: v.data)
    # print("VOL:" + str(vol.shape))
    # return np.average(fo, weights=vol)
 
 
def getStatsForFrame(frame, assembly, stepLoading):
    hasPatch = stepLoading[3] != 0 # means we don't have to worry about weird float issues
    frameLoading = stepLoading[0:3] * frame.frameValue
    fo = frame.fieldOutputs
    # instances
    instances = assembly.instances
    crack_proximity = instances["SURFACE-TOP"].elementSets["CRACK-ANALYSIS"]
    crack = instances["SURFACE-TOP"].elementSets["SURFACE-WEAK"]
    surfaceElem = instances["SURFACE-TOP"].elementSets["PATCH-NO-CRACK"]
    surfaceNode = instances["SURFACE-TOP"].nodeSets["PATCH-NO-CRACK"]
    adhesive = None if ~hasPatch else instances["ADHESIVE-LAYER-1"]
    patch = None if ~hasPatch else instances["PATCH-1"]
    stretch = assembly.nodeSets["STRETCH"]
    # surface strains n rotation
    sE11 =  volumeMean(fo, "LE", surfaceElem, lambda v: v.data[0])
    sE12 =  volumeMean(fo, "LE", surfaceElem, lambda v: v.data[3])
    sUR =  np.mean(fieldComponent(fo, "UR", surfaceNode, lambda v: v.data[0]))
    patchStressMax = np.max(fieldComponent(fo, "S", patch, lambda v: v.mises))
    surfaceStressMax = np.max(fieldComponent(fo, "S", surfaceElem, lambda v: v.mises))
    crackStressMax = np.max(fieldComponent(fo, "S", crack_proximity, lambda v: v.mises))
    crackStrainE11 = np.max(fieldComponent(fo, "LE", crack, lambda v: v.data[0]))
    crackStrainEMax = np.max(fieldComponent(fo, "LE", crack, lambda v: v.maxInPlanePrincipal))
    crackStrainE22 = np.max(fieldComponent(fo, "LE", crack, lambda v: v.data[1]))
    crackStrainE12 = np.max(fieldComponent(fo, "LE", crack, lambda v: v.data[3]))
    adhesiveDegradationMax = np.max(fieldComponent(fo, "SDEG", adhesive, lambda v: v.data))
    adhesiveDegradationAvg= np.mean(fieldComponent(fo, "SDEG", adhesive, lambda v: v.data))
    sRFX = fieldComponent(fo, "RF", stretch, lambda v: v.data[0])
    sRFY = fieldComponent(fo, "RF", stretch, lambda v: v.data[1])
    sRFZ = fieldComponent(fo, "RF", stretch, lambda v: v.data[2])
    sPY =  fieldComponent(fo, "COORD", stretch, lambda v: v.data[1]) + fieldComponent(fo, "U", stretch, lambda v: v.data[1])
    sPZ =  fieldComponent(fo, "COORD", stretch, lambda v: v.data[2]) + fieldComponent(fo, "U", stretch, lambda v: v.data[2])
    stretchRMX = sPZ * sRFY - sPY * sRFZ
    # hashin it
    hashinFC = np.max(fieldComponent(fo, "HSNFCCRT", patch, lambda v: v.data))
    hashinFT = np.max(fieldComponent(fo, "HSNFTCRT", patch, lambda v: v.data))
    hashinMC = np.max(fieldComponent(fo, "HSNMCCRT", patch, lambda v: v.data))
    hashinMT = np.max(fieldComponent(fo, "HSNMTCRT", patch, lambda v: v.data))
    hashinModes = [hashinFC, hashinFT, hashinMC, hashinMT]
    hashinMax = max(hashinModes)
    return {
        "t"  : frame.frameValue,
        "hasPatch": hasPatch,
        "dx" : frameLoading[0],
        "dz" : frameLoading[1],
        "dr" : frameLoading[2], 
        "E11": sE11,
        "E12": sE12,
        "UR1": sUR,
        "RFX": np.sum(sRFX),
        "RFY": np.sum(sRFY),
        "RFZ": np.sum(sRFZ),
        "RMX": np.sum(stretchRMX),
        "AdhDMax": adhesiveDegradationMax,
        "AdhDAvg": adhesiveDegradationAvg,
        "PtchSMax": patchStressMax,
        "SurfSMax": surfaceStressMax,
        "CrckSMax": crackStressMax,
        "CrckE11": crackStrainE11,
        "CrckE12": crackStrainE12,
        "CrckE22": crackStrainE22,
        "CrckEMax": crackStrainEMax,
        "HshnMax": hashinMax,
        "HshnType": hashinModes.index(hashinMax),
    }
 
 
def getStatsForStep(step, assembly, stepLoading):
    return[getStatsForFrame(frame, assembly, stepLoading) for frame in step.frames]
 
 
def dictToCsvString(data: list[dict[str,Any]]):
    return "\n".join([", ".join([f"{value}" for value in data[0].keys()])] + \
              [", ".join([f"{value:+.8e}" for value in row.values()]) for row in data])
 
 
# jobs = writeInputsInRange(
#     mdb, model, np.linspace(0, 0.7, 8), np.linspace(0, 1, 8), np.linspace(0, pi/18, 8), [True, False]
# )
jobs = writeInputsInRange(
    mdb, model, np.linspace(0, 0.5, 9), np.linspace(0, 0.9, 9), np.linspace(0, pi/18, 9), [True, False]
)
 
 
# uniaxial cases:
if(False):
    jobs = {}
    jobs.update(writeInputsInRange(mdb, model, [1], [0], [0], [False, True]))
    jobs.update(writeInputsInRange(mdb, model, [0], [1], [0], [False, True]))
    jobs.update(writeInputsInRange(mdb, model, [0], [0], [pi/12], [False, True]))
 
 
def writeOdbStats(jobName, job, batchFilename):
    count_errored = 0
    odb = None
    try:
        odb = openOdb("./" + jobName + ".odb", readOnly=True)
        step = odb.steps["Step-1"]
        loading = np.array(jobs[jobName][1:5])
        csvstring = dictToCsvString(getStatsForStep(step, odb.rootAssembly, loading))
    except Exception:
        count_errored+=1
        extxt = traceback.format_exc()
        print(extxt)
        if odb: odb.close()
        return 1
        raise
    with open("../csvs2/" + jobName + "-data.csv", "w") as file:
        file.write(csvstring)
    if batchFilename:
        with open("../csvs2/" + batchFilename + ".csv", "a") as file:
            file.write("\n"+csvstring)
    with open("../csvs2/collated-data.csv", "a") as file:
        file.write("\n"+csvstring)
    if odb: odb.close()
    return 0
 
 
n = 6
count_total = len(jobs)
count_completed = 0
count_errored = 0
pending_jobs = list(jobs.keys())
stats = {}
batchFilename = "Data-" + datetime.datetime.now().strftime("%Y-%m-%dT%H_%M_%S")
while pending_jobs:
    jobsToComplete: list[tuple[str,Any]] = []
    jobsAdded = 0
    while pending_jobs and jobsAdded < 6:
        jobName = pending_jobs.pop()
        if os.path.exists("../csvs2/" + jobName + "-data.csv"):
            print(f"Skipped submitting {jobName}. data exists")
            continue # don't redo existing jobs
        # todo rewrite
        job = jobs[jobName][0]
        print(f"submitting {jobName}")
        job.submit()
        jobsToComplete.append((jobName,job))
        jobsAdded+=1
    for jobName, job in jobsToComplete:
        job.waitForCompletion()
        count_completed+=1
        odb = None
        count_errored += writeOdbStats(jobName, jobs[jobName], batchFilename)
    print(f"Jobs completed {count_completed}/{count_total} ({count_errored} errors)")
 

Link to original

AbaqusMacrosPy

# -*- coding: mbcs -*-
# Do not delete the following import lines
from abaqus import *
from abaqusConstants import *
import __main__
 
FRAME_FOR_LOADING = {
    "+0_50X+0_00Z+0_00R": 8,
    "+0_00X+0_90Z+0_00R": 9,
    "+0_00X+0_00Z+10_00R": 17,
    "+0_12X+0_90Z+2_50R": 8,
    "+0_50X+0_00Z+10_00R": 9
    }
 
def createNextViewport():
    maxNum = max(map(lambda vp: int(vp.split(" ")[1]), session.viewports.keys()))
    session.Viewport(f"Viewport: {maxNum+1}")
 
def manageViewportCount(desiredCount):
    print(f"Want {desiredCount} viewports and have {len(session.viewports)}")
    sortedKeys = sorted(session.viewports.keys(), key = lambda vp: int(vp.split(" ")[1]))
    if desiredCount < len(session.viewports) and False: # actually don't ever delete viewports
        session.viewports[sortedKeys[0]].makeCurrent()
        print("Removing excess viewports")
        for i, key in enumerate(sortedKeys):
            if i >= desiredCount:
                del session.viewports[key]
    elif desiredCount > len(session.viewports):
        print("Adding needed viewports")
        for i in range(len(session.viewports), desiredCount):
            createNextViewport()
    # re-sort
    sortedKeys = sorted(session.viewports.keys(), key = lambda vp: int(vp.split(" ")[1]))
    vps = {i:session.viewports[i] for i in sortedKeys}
    return vps
 
def ShowFieldOutput(viewportSettings):
    
    vps = manageViewportCount(len(viewportSettings))
    usedViewports = [];
                     
    for settings, viewport in zip(viewportSettings,vps.values()):
        viewport.view.setProjection(projection=PARALLEL)
        viewport.odbDisplay.contourOptions.setValues(maxAutoCompute=ON, minAutoCompute=ON)
        viewport.makeCurrent()
        viewport.odbDisplay.contourOptions.setValues(
        colorByMatchingPlies=ON if settings[3] else OFF)
        if settings[1]:
            viewport.odbDisplay.setPrimaryVariable(
                variableLabel=settings[0],
                outputPosition=INTEGRATION_POINT,
                refinement=settings[1]
                )
        else:
            viewport.odbDisplay.setPrimaryVariable(
                variableLabel=settings[0],
                outputPosition=INTEGRATION_POINT
                )
        
        viewport.odbDisplay.basicOptions.setValues(
            sectionResults=USE_ENVELOPE,
            envelopeCriteria=settings[2]
            ) 
        viewport.odbDisplay.display.setValues(plotState=(CONTOURS_ON_UNDEF,))
        usedViewports.append(viewport)
        
    
    # get mins and maxes
    legendGroups = set([settings[4] for settings in viewportSettings])
    for lg in legendGroups:
        if lg is None: continue
        connectedVPsIx = [i for i,settings in enumerate(viewportSettings) if settings[4] == lg]
        connectedVPs = [usedViewports[i] for i in connectedVPsIx]
        lgMin = 1e9
        lgMax = -1e9
        for vp in connectedVPs:
            lgMax = max(lgMax, vp.odbDisplay.contourOptions.autoMaxValue)
            lgMin = min(lgMax, vp.odbDisplay.contourOptions.autoMinValue)
        # now that we have min/max set to that
        for vp in connectedVPs:
            vp.odbDisplay.contourOptions.setValues(
                maxAutoCompute=OFF, maxValue=lgMax, minAutoCompute=OFF, minValue=lgMin)
            
    
    return usedViewports
        
 
def FOStress():
    viewportSettings = [
        ("S", (INVARIANT, "Mises"), MAX_VALUE, False, 0),
        ("S", (INVARIANT, "Mises"), MIN_VALUE, False, 0),
        ("S", (COMPONENT, "S11"), MAX_VALUE, False, 1),
        ("S", (COMPONENT, "S11"), MIN_VALUE, False, 1),
        ("S", (COMPONENT, "S22"), MAX_VALUE, False, 2),
        ("S", (COMPONENT, "S22"), MIN_VALUE, False, 2),
        ("S", (COMPONENT, "S12"), MAX_VALUE, False, 3),
        ("S", (COMPONENT, "S12"), MIN_VALUE, False, 3)
        ]
    return ShowFieldOutput(viewportSettings)
 
    
def FOHashin():
    viewportSettings = [
        ("HSNFCCRT", None, MAX_VALUE, False, None),
        ("HSNFCCRT", None, MIN_VALUE, True, None),
        ("HSNFTCRT", None, MAX_VALUE, False, None),
        ("HSNFTCRT", None, MIN_VALUE, True, None),
        ("HSNMCCRT", None, MAX_VALUE, False, None),
        ("HSNMCCRT", None, MIN_VALUE, True, None),
        ("HSNMTCRT", None, MAX_VALUE, False, None),
        ("HSNMTCRT", None, MIN_VALUE, True, None)
        ]
    return ShowFieldOutput(viewportSettings)
    
    
def FOAdhesive():
    viewportSettings = [
        ("SDEG", None, MAX_VALUE, False, None),
        ("S", (INVARIANT, "Mises"), MAX_VALUE, False, None)
        ]
    return ShowFieldOutput(viewportSettings)
    
    
def ViewSettings():
    import connectorBehavior
    viewport = session.viewports.values()[0]
    viewport.viewportAnnotationOptions.setValues(title=OFF, 
        state=OFF, compass=OFF)
    viewport.viewportAnnotationOptions.setValues(
        triadPosition=(95, 90))
    viewport.viewportAnnotationOptions.setValues(
        legendFont='-*-verdana-medium-r-normal-*-*-70-*-*-p-*-*-*')
 
 
def StressSurface():
    session.viewports[session.currentViewportName].odbDisplay.contourOptions.setValues(
        maxAutoCompute=OFF, maxValue=345, minAutoCompute=OFF, minValue=-345)
        
        
def StressPatch():
    session.viewports[session.currentViewportName].odbDisplay.contourOptions.setValues(
        maxAutoCompute=OFF, maxValue=400, minAutoCompute=OFF, minValue=-400)
 
 
def ContourColors():
    session.viewports[session.currentViewportName].odbDisplay.contourOptions.setValues(
        outsideLimitsAboveColor='#FF00F8', outsideLimitsBelowColor='#6C008A')
 
 
def SaveToImage(imageInfo, viewports):
    size = (3840, 2010)# (2560, 1340)
    session.printOptions.setValues(vpDecorations=OFF, reduceColors=False)
    session.pngOptions.setValues(imageSize = size)
    odbName = session.viewports.values()[0].displayedObject.name.split("/")[-1].strip(".odb");
    loading = odbName.strip("Patch").strip("less");
    hasPatch = not "less" in odbName
    fileName = loading + "-" + imageInfo
    if imageInfo == "SurfaceStress":
        fileName +=  "-Patch" if hasPatch else "-Patchless"
    filePatch = 'C:/Users/S3943498/OneDrive - RMIT University/Documents/Capstone personal/Patch/images/' + fileName
    session.printToFile(
        fileName=filePatch, 
        format=PNG, canvasObjects=tuple(viewports))
    print(f"Saving {imageInfo} to ./images/{fileName}")
 
 
def ShowSurface():
    import displayGroupOdbToolset as dgo
    leaf = dgo.LeafFromElementSets(elementSets=('SURFACE-TOP.REMESH', ))
    session.viewports[session.currentViewportName].odbDisplay.displayGroup.replace(leaf=leaf)
    session.viewports[session.currentViewportName].view.setValues(session.views['User-1'])
 
 
 
def ShowPatch():
    import displayGroupOdbToolset as dgo
    leaf = dgo.LeafFromPartInstance(partInstanceName=('PATCH-1', ))
    session.viewports[session.currentViewportName].odbDisplay.displayGroup.replace(leaf=leaf)
    session.viewports[session.currentViewportName].view.setValues(session.views['User-3'])
 
 
def ShowAdhesive():
    import displayGroupOdbToolset as dgo
    leaf = dgo.LeafFromPartInstance(partInstanceName=('ADHESIVE-LAYER-1', ))
    session.viewports[session.currentViewportName].odbDisplay.displayGroup.replace(leaf=leaf)
    session.viewports[session.currentViewportName].view.setValues(session.views['User-2'])
 
    
def APatchStress():
    ShowPatch()
    SaveToImage("PatchStress", FOStress())
    
    
def APatchHashin():
    ShowPatch()
    SaveToImage("PatchHashin", FOHashin())
    
    
def ASurfaceStress():
    ShowSurface()
    SaveToImage("SurfaceStress", FOStress())
    
def AAdhesiveAll():
    ShowAdhesive()
    SaveToImage("AdhesiveAll", FOAdhesive())
 
def AAAGetImagesForOdb(odb):
    odbName = odb.name.split("/")[-1].strip(".odb")
    loading = odbName.strip("Patch").strip("less")
    # get frame for this odb
    if loading in FRAME_FOR_LOADING:
        frameIx = FRAME_FOR_LOADING[loading]
        frame = odb.steps['Step-1'].frames[frameIx]
    else:
        frameIx = len(odb.steps['Step-1'].frames)-1
        
    print(odbName + " : " + loading + ", Frame " + str(frameIx))
    for viewport in session.viewports.values():
        viewport.setValues(displayedObject = odb)
        viewport.odbDisplay.setFrame(step='Step-1', frame=frameIx) 
    manageViewportCount(8)
    if "less" in odbName:
        # patchless, only do surface
        ASurfaceStress()
    else:
        ASurfaceStress()
        APatchHashin()
        APatchStress() # have to do stress before changing element set
        AAdhesiveAll()
    
def AAAGetImagesForViewportOdb():
    odb = session.viewports[session.currentViewportName].displayedObject
    AAAGetImagesForOdb(odb)
 
 
def AAAAGetImagesForLoadedOdbs():
    for odb in session.odbs.values():
        AAAGetImagesForOdb(odb)
 
 
def ZZImports():
    import section
    import regionToolset
    import displayGroupMdbToolset as dgm
    import part
    import material
    import assembly
    import step
    import interaction
    import load
    import mesh
    import optimization
    import job
    import sketch
    import visualization
    import xyPlot
    import displayGroupOdbToolset as dgo
    import connectorBehavior
    pass
 

Link to original

Link to original