ab: add code

Baldit_2021_JMBBM_Dorfmann_and_Ogden_2004_script.py

+# Copyright University of Lorraine - LEM3
+# Contributor(s) : 
+#         Adrien Baldit
+# Contact: adrien.baldit@univ-lorraine.fr
+# 
+# This script is a computer program whose purpose is to produce 
+# biomechanical and bioengineering data processing.
+#
+# This script is governed under French law and abiding by the rules 
+# of distribution of free script.  You can  use,  modify and/ or 
+# redistribute the script under the terms of honor
+#
+# As a counterpart to the access to the source code and  rights to copy,
+# modify and redistribute granted by the honor, users are provided only
+# with a limited warranty  and the script's author,  the holder of the
+# economic rights,  and the successive licensors  have only  limited
+# liability. 
+#
+# In this respect, the user's attention is drawn to the risks associated
+# with loading,  using,  modifying and/or developing or reproducing the
+# script by the user in light of its specific status of free software,
+# that may mean  that it is complicated to manipulate,  and  that  also
+# therefore means  that it is reserved for developers  and  experienced
+# professionals having in-depth computer knowledge. Users are therefore
+# encouraged to load and test the script's suitability as regards their
+# requirements in conditions enabling the security of their systems and/or 
+# data to be ensured and,  more generally, to use and operate it in the 
+# same conditions as regards security. 
+# 
+# The fact that you are presently reading this means that you have 
+# had knowledge of the rules and accepted them.
+
+#!/usr/bin/python
+# -*- coding:utf8 -*-
+
+######################
+#     Libraries import
+######################
+
+import numpy as np
+import pylab as pl
+import os
+
+###############
+# Pylab parameter
+##################
+#~ pl.rc('font', size=18)          # controls default text sizes
+pl.rc('axes', titlesize=25)     # fontsize of the axes title
+pl.rc('axes', labelsize=25)    # fontsize of the x and y labels
+pl.rc('xtick', labelsize=20)    # fontsize of the tick labels
+pl.rc('ytick', labelsize=20)    # fontsize of the tick labels
+pl.rc('legend', fontsize=20)    # legend fontsize
+pl.rc('lines', linewidth=2)    # legend fontsize
+pl.rc('lines', markersize=10)    # legend fontsize
+pl.rc('figure', titlesize=18)  # fontsize of the figure title
+pl.rc('text', usetex=True)
+
+# Interactive plotting
+pl.ion()
+
+######################
+#       Parameters
+######################
+
+dict_data = {
+            "Dorfmann and Ogden 2004" : {
+                "alpha_i" : [-1.011467,4.2047799,-4.398598],\
+                "mu_i MPa" : [-1.52838,0.222564,-1.13418e-3],\
+                "r" : 1.25,\
+                "m" : 0.965,\
+                "q" : 0.4,\
+                "maximum stretch" : 3.,\
+                },\
+            "Wharton's jelly" : {
+                "alpha_i" : [0.606683,7.770503,6.736465],\
+                "mu_i MPa" : [-0.63024,0.052257,0.034789],\
+                "r" : 0.865844,\
+                "m" : 0.006835,\
+                "q" : 0.162428,\
+                "maximum stretch" : 1.078783,\
+                },\
+            }
+
+########################################
+#  Dorfmann and Ogden 2004 Model
+########################################
+
+def dorfmann_ogden_uniaxial(stretch_list,mu_i,alpha_i, r, m,coef) :
+    """ A. Dorfmann and R.W. Ogden. A constitutive model for the Mullins effect with permanent set in particle-reinforced rubber. International Journal of Solids and Structures 41 (2004) 1855-1878
+    The stretch_list is the same for loading and unloading phases"""
+    # Stretch list maximum
+    stretch_max = max(stretch_list)
+    # initialisation of the unloading strain energy density
+    Wsul = pl.zeros(len(stretch_list))
+    # initialisation of the strain enery density for the maximum stretch
+    Wsm = 0
+    # mu parameter depending on Ogden model parameters
+    mu = 0
+    # initialisation of the load stress vector
+    tl = pl.zeros(len(stretch_list))
+    # loop over Ogden's model parameters
+    for j in pl.arange(len(mu_i)) :
+        # extending load stress vector
+        tl += mu_i[j]*(stretch_list**(alpha_i[j]-1.)-stretch_list**(-alpha_i[j]/2.-1.))
+        # extending loading strain energy density
+        Wsul += mu_i[j]/alpha_i[j]*(stretch_list**alpha_i[j]+2.*stretch_list**(-alpha_i[j]/2.)-3.)
+        # extending strain enery density for the maximum stretch
+        Wsm += mu_i[j]/alpha_i[j]*(stretch_max**alpha_i[j]+2.*stretch_max**(-alpha_i[j]/2.)-3.)
+        # extending mu parameter
+        mu += 0.5*mu_i[j]*alpha_i[j]
+    # computation of the alpha parameter
+    alpha = 0.3+0.16/mu*Wsm
+    # computation of the nu parameters
+    nu_1 = coef*mu*(1.-1/3.5*pl.tanh((stretch_max-1.)/0.1))
+    nu_2 = coef*mu
+    # compuatation of the eta parameters
+    eta_1 = 1.-1./r*pl.tanh((Wsm-Wsul)/(m*mu))
+    eta_2 = pl.tanh((Wsul/Wsm)**alpha)/pl.tanh(1.)
+    # unload stress vector
+    tul = eta_1*tl + (1.-eta_2)*(nu_1*stretch_list-nu_2*stretch_list**(-2.))
+    return tl, tul
+
+########################################
+#  Plot container
+########################################
+fig = pl.figure(figsize=(16,8))
+# subplot parameters
+fig.subplots_adjust(left = 0.07, bottom = 0.1, right = 0.97, top = 0.94, wspace = 0.22, hspace = 0.03)
+
+########################################
+# Loop over dataset
+########################################
+
+for ind,data_type in enumerate(dict_data.keys()) :
+    # Data
+    data = dict_data[data_type]
+    # Ogden model parameters
+    alpha_i = pl.array(data["alpha_i"])
+    mu_i = pl.array(data["mu_i MPa"])
+    # control of parameter consitency
+    if len(mu_i) != len(alpha_i):
+        print("ERROR in Ogden model parameters")
+    else :
+        mu = (0.5*mu_i*alpha_i).sum()
+        print("mu = {} MPa".format(mu))
+    # Dorfmann & Ogden parameters
+    r = data["r"]
+    m = data["m"]
+    coef = data["q"]
+
+    # Maximum stretch applied on the sample
+    stretch_max = data["maximum stretch"]
+
+    # Stretch sampling
+    Dstretch = (stretch_max-1.)/100.
+    stretch_list_load = pl.arange(1.,stretch_max+Dstretch,Dstretch)
+
+    # stress computation
+    tl,tul = dorfmann_ogden_uniaxial(stretch_list = stretch_list_load,\
+                                        mu_i = mu_i,\
+                                        alpha_i = alpha_i,\
+                                        r = r,\
+                                        m = m,\
+                                        coef = coef)
+
+    # subfigure
+    ax = fig.add_subplot(121+ind)
+    # title
+    ax.set_title(data_type)
+    # plot
+    ax.plot(stretch_list_load, tl,"-bx",label="loading "+"$\mu = {}\ MPa$".format(np.round(mu,2)),markevery = 2)
+    ax.plot(stretch_list_load, tul,"-bo",label="unloading",markevery = 2)
+    # axis label
+    ax.set_ylabel(r"$Nominal\ stress\ [MPa]$")
+    ax.set_xlabel(r"$stretch\ [-]$")
+    ax.legend()
+    ax.grid()
+# save figure
+fig.savefig(os.getcwd()+os.sep+"Baldit_2021_JMBBM_Dorfmann_and_Ogden_2004.png", dpi=150)