Dorfmann and Ogden 2004 script

# 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+"Dorfmann_and_Ogden_2004_model.pdf", dpi=None, facecolor='w', edgecolor='w', orientation='portrait', format = 'pdf', transparent=False, bbox_inches=None, pad_inches=0.1)