Import source code

main.py

@@ -0,0 +1,45 @@
+# TP 1 Pasticité - Jovian Hersemeule
+print("PLARU starting ...")
+
+import scipy as sc
+from plaruCompute import radialReturn
+from plaruPlot import fullPlot
+
+# Parameters
+E = 1.0
+sigY = 2.0
+H = 0.25
+
+# test
+# size = 50
+# deps = 0.1 * sc.ones(size)
+
+# a)
+# d1 = 0.1 * sc.ones(25)
+# d2 = -0.1 * sc.ones(21)
+# deps = sc.concatenate((d1, d2))
+
+# b)
+# d1 = 0.1 * sc.ones(25)
+# d2 = -0.1 * sc.ones(45)
+# d3 = 0.1 * sc.ones(20)
+# deps = sc.concatenate((d1, d2, d3))
+
+# c) 1
+d1 = 0.1 * sc.ones(15)
+d2 = -0.1 * sc.ones(15)
+d3 = 0.1 * sc.ones(25)
+d4 = -0.1 * sc.ones(21)
+d5 = 0.1 * sc.ones(35)
+d6 = -0.1 * sc.ones(24)
+deps = sc.concatenate((d1, d2, d3, d4, d5, d6))
+
+# c) 2
+# deps = sc.concatenate((d5, d6, d3, d4, d1, d2))
+
+
+# Compute results
+(eps, sig, epsP) = radialReturn(deps, E, sigY, H = H)
+
+# Plot results
+fullPlot(eps, sig, epsP, deps)

plaruCompute.py

@@ -0,0 +1,44 @@
+# Implements plasticity algorithm
+
+import scipy as sc
+
+def radialReturn(deps, E, sigY, H = 0):
+    """
+    Simulate a traction experience.
+    :param deps: A iteration vector for epsilon.
+    :param E: Hook factor.
+    :param sigY: Max elastic stress.
+    :param H: Ecrouissage.
+    :return: Three vectors (eps, sig, epsP)
+    """
+    # Announcement
+    print("Computing results ...")
+
+    # Prepare results
+    n = len(deps) + 1
+    eps = sc.zeros(n)
+    sig = sc.zeros(n)
+    epsP = sc.zeros(n)
+
+    # Algorithm loop
+    for k in range(0, n - 1):
+        # Add deformation
+        eps[k + 1] = eps[k] + deps[k]
+
+        # Plastic prediction
+        sig_pred = E*(eps[k + 1] - epsP[k])
+        f_pred = abs(sig_pred - H*epsP[k]) - sigY
+
+        # Plasticity test
+        if f_pred < 0.0 :
+            # Elastic domain
+            epsP[k + 1] = epsP[k]
+            sig[k + 1] = sig_pred
+        else:
+            # Radial return
+            lda = f_pred / (E + H)
+            epsP[k + 1] = epsP[k] + sc.sign(sig_pred - H*epsP[k])*lda
+            sig[k + 1] = E*(eps[k + 1] - epsP[k + 1])
+
+    # End of computation
+    return (eps, sig, epsP)

plaruPlot.py

@@ -0,0 +1,54 @@
+import matplotlib.pyplot as plt
+import scipy as sc
+
+def fullPlot(eps, sig, epsP, deps):
+    print("Ploting results ...")
+
+    fig, axs = plt.subplots(1, 4, figsize=(20,9))
+    plt.subplots_adjust(wspace=0.5)
+    plt.rcParams.update({'font.size': 24})
+
+    # Draw axis
+    for k in range(0, 4):
+        axs[k].grid()
+        axs[k].axvline(color = 'black')
+        axs[k].axhline(color = 'black')
+
+    # Plot 1
+    axs[0].title.set_text("$\\varepsilon \\to \sigma$")
+    axs[0].set_xlabel("$\\varepsilon$")
+    axs[0].set_ylabel("$\sigma \, (MPa)$")
+    axs[0].plot(eps, sig)
+
+    # Plot 2
+    axs[1].title.set_text("$\\varepsilon \\to \\varepsilon_p$")
+    axs[1].set_xlabel("$\\varepsilon$")
+    axs[1].set_ylabel("$\\varepsilon_p$")
+    axs[1].plot(eps, epsP)
+
+    # Plot 3
+    axs[2].title.set_text("$\\varepsilon_p \\to \sigma$")
+    axs[2].set_xlabel("$\\varepsilon_p$")
+    axs[2].set_ylabel("$\sigma \, (MPa)$")
+    axs[2].plot(epsP, sig)
+
+    # Plot 4
+    axs[3].title.set_text("$k \\to d\\varepsilon_p$")
+    axs[3].set_xlabel("k")
+    axs[3].set_ylabel("$d\\varepsilon_p$")
+    axs[3].plot(deps, marker='+')
+
+    # Better scale
+    xadd = 0.2
+    yadd = 0.2
+    for k in range(0, 4):
+        lft, rig = axs[k].get_xlim()
+        lft -= xadd
+        axs[k].set_xlim(lft, rig)
+
+        lft, rig = axs[k].get_ylim()
+        lft -= yadd
+        rig += yadd
+        axs[k].set_ylim(lft, rig)
+
+    plt.show()