Taylor bar example

This example is inspired by the “Taylor Bar” example on the LS-DYNA Examples site. It shows how to use DYNA-Lib to create a keyword file for LS-DYNA and solve it within a Pythonic environment.

Perform required imports

Import required packages, including those for the keywords, deck, and solver.

# sphinx_gallery_thumbnail_number = 1

import datetime
import os
import tempfile
import time

import matplotlib.pyplot as plt
import pandas as pd

from ansys.dyna import examples
from ansys.dyna.keywords import Deck
from ansys.dyna.keywords import keywords as kwd
from ansys.dyna.run import run_dyna

Import DPF-Core

Import the PyDPF-Core package.

try:
    from ansys.dpf import core as dpf

    HAS_DPF = True
except ImportError:
    HAS_DPF = False

Create a deck and keywords

Create a deck, which is the container for all the keywords. Then, create and append individual keywords to the deck.

def create_input_deck(**kwargs):
    initial_velocity = kwargs.get("initial_velocity")
    wd = kwargs.get("wd")
    if not all((initial_velocity, wd)):
        raise Exception("Missing input!")

    deck = Deck()
    deck.title = "Taylor-Bar Velocity - %s - Unit: t-mm-s" % initial_velocity

    # Import mesh
    mesh = examples.taylor_bar_mesh
    include1 = kwd.Include(filename=mesh)

    # Define material
    mat_1 = kwd.Mat003(mid=1)
    mat_1.ro = 7.85000e-9
    mat_1.e = 150000.0
    mat_1.pr = 0.34
    mat_1.sigy = 390.0
    mat_1.etan = 90.0

    # Define section
    sec_1 = kwd.SectionSolid(secid=1)
    sec_1.elform = 1

    # Define part
    part_1 = kwd.Part()
    part_1.parts = pd.DataFrame({"pid": [1], "mid": [mat_1.mid], "secid": [sec_1.secid]})

    # Define coordinate system
    cs_1 = kwd.DefineCoordinateSystem(cid=1)
    cs_1.xl = 1.0
    cs_1.yp = 1.0

    # Define initial velocity
    init_vel = kwd.InitialVelocityGeneration()
    init_vel.id = part_1.parts["pid"][0]
    init_vel.styp = 2
    init_vel.vy = initial_velocity
    init_vel.icid = cs_1.cid

    # Define box for node set
    box_1 = kwd.DefineBox(boxid=1, xmn=-500, xmx=500, ymn=39.0, ymx=40.1, zmn=-500, zmx=500)

    # Create node set
    set_node_1 = kwd.SetNodeGeneral()
    set_node_1.sid = 1
    set_node_1.option = "BOX"
    set_node_1.e1 = box_1.boxid

    # Define rigid wall
    rw = kwd.RigidwallPlanar(id=1)
    rw.nsid = set_node_1.sid
    rw.yt = box_1.ymx
    rw.yh = box_1.ymn

    # Define control termination
    control_term = kwd.ControlTermination(endtim=8.00000e-5, dtmin=0.001)

    # Define database cards
    deck_dt_out = 8.00000e-8
    deck_glstat = kwd.DatabaseGlstat(dt=deck_dt_out, binary=3)
    deck_matsum = kwd.DatabaseMatsum(dt=deck_dt_out, binary=3)
    deck_nodout = kwd.DatabaseNodout(dt=deck_dt_out, binary=3)
    deck_elout = kwd.DatabaseElout(dt=deck_dt_out, binary=3)
    deck_rwforc = kwd.DatabaseRwforc(dt=deck_dt_out, binary=3)
    deck_d3plot = kwd.DatabaseBinaryD3Plot(dt=4.00000e-6)

    # Define deck history node
    deck_hist_node_1 = kwd.DatabaseHistoryNodeSet()
    deck_hist_node_1.id1 = set_node_1.sid

    # Append all cards to input deck
    deck.extend(
        [
            deck_glstat,
            deck_matsum,
            deck_nodout,
            deck_elout,
            deck_rwforc,
            deck_d3plot,
            set_node_1,
            control_term,
            rw,
            box_1,
            init_vel,
            cs_1,
            part_1,
            mat_1,
            sec_1,
            deck_hist_node_1,
            include1,
        ]
    )
    deck_string = deck.write()

    # Create LS-DYNA input deck
    os.makedirs(wd, exist_ok=True)
    file_handle = open(os.path.join(wd, "input.k"), "w")
    file_handle.write(deck_string)
    file_handle.close()

    return deck

Define the Dyna solver function

def run(**kwargs):
    wd = kwargs.get("wd")
    inputfile = os.path.join(wd, "input.k")
    run_dyna(inputfile)
    assert os.path.isfile(os.path.join(wd, "d3plot"))
    assert os.path.isfile(os.path.join(wd, "lsrun.out.txt"))

Define the DPF output function

def get_global_ke(**kwargs):
    wd = kwargs.get("wd")
    ds = dpf.DataSources()
    ds.set_result_file_path(os.path.join(wd, "d3plot"), "d3plot")
    model = dpf.Model(ds)

    ts = [set for set in range(1, model.metadata.time_freq_support.n_sets + 1)]
    ke_op = model.results.global_kinetic_energy(time_scoping=ts)

    time = model.metadata.time_freq_support.time_frequencies.data_as_list
    ke = []
    for t in ts:
        ke.append(ke_op.outputs.global_kinetic_energy.get_data()[t - 1].data[0])

    return time, ke


# Define base working directory
root_out_dir = tempfile.gettempdir()

stamp = datetime.datetime.fromtimestamp(time.time()).strftime("%Y-%m-%d_%H_%M_%S")
base_wd = os.path.join(root_out_dir, "PyDyna.%s" % stamp)
os.mkdir(base_wd)

View the model

etc etc TODO Need to support solids for this to work

deck_for_graphic = create_input_deck(initial_velocity=300e3, wd="run")
# print(deck_for_graphic.expand())
deck_for_graphic.plot(cwd="run")

Run a parametric solve

etc etc

color = ["b", "r", "g", "y"]
# Specify different velocities in mm/s
initial_velocities = [275.0e3, 300.0e3, 325.0e3, 350.0e3]

for index, initial_velocity in enumerate(initial_velocities):
    # Create a folder for each parameter
    wd = os.path.join(base_wd, "tb_vel_%s" % initial_velocity)
    os.mkdir(wd)
    # Create LS-Dyna input deck
    create_input_deck(initial_velocity=initial_velocity, wd=wd)
    # Run Solver
    try:
        run(wd=wd)
        # Run PyDPF Post
        if HAS_DPF:
            time, ke = get_global_ke(initial_velocity=initial_velocity, wd=wd)
            # Add series to the plot
            plt.plot(time, ke, color[index], label="KE at vel. %s mm/s" % initial_velocity)
    except Exception as e:
        print(e)

        # sphinx_gallery_defer_figures

plt.xlabel("Time (s)")
plt.ylabel("Energy (mJ)")

Text(0, 0.5, 'Energy (mJ)')

Generate graphical output

etc etc

plt.show()