# 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)