b2composite_consistent_damage_model_mcmd_2006.H

//------------------------------------------------------------------------
// b2composite_consistent_damage_model_mcmd_2006.H --
//
//
// written by Mathias Doreille
//
// Copyright (c) 2006-2012,2016
//     SMR Engineering & Development SA
//     2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef __B2COMPOSITE_CONSISTENT_DAMAGE_MODEL_MCMD_2006_H__
#define __B2COMPOSITE_CONSISTENT_DAMAGE_MODEL_MCMD_2006_H__
 
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <memory>
 
#include "utils/b2exception.H"
#include "utils/b2type.H"
 
namespace b2000 {
// Implement the model described in the following NASA Technical
// Memorandum:
//
// Maimí P., Camanho P.P., Mayugo J.A., Dávila C.G., 2006.
// A thermodynamically consistent damage model for advanced composites.
// NASA/TM-2006-214282
//
// Available at:
// http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20060008658_2006009390.pdf
//
// G : G12, G13, G23
// X : XT, XC, YT, YC, SL
// Gf : G1+, G1-, G2+, G2-, G6
// r : r1+, r1-, r2+, r2- ; must be set at 1 (no damage) at begening
//
// The viscous regularisation was applied to all the damages
// and not only to the longitudinal direction has in the Memorandum.
//
class CompositeConsistentDamageModelMCMD2006 {
public:
    CompositeConsistentDamageModelMCMD2006() : damage_data(nullptr) {}
 
    ~CompositeConsistentDamageModelMCMD2006() { delete damage_data; }
 
    void get_stress(
          const double E[2], const double P, const double G[3], const double X[5],
          const double Gf[5], const double area, const double b, const double time,
          const double delta_time, const double viscosity, const double strain[6],
          const EquilibriumSolution equilibrium_solution, double stress[6],
          double d_stress_d_strain[21], int& elasticity_domain);
 
    const double* get_damages() const {
        if (damage_data) {
            return damage_data->r;
        } else {
            static const double undamage_r[4] = {1.0, 1.0, 1.0, 1.0};
            return undamage_r;
        }
    }
 
    const double* get_A() const {
        if (damage_data) {
            return damage_data->A;
        } else {
            static const double undamage_A[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
            return undamage_A;
        }
    }
 
private:
    struct DamageData {
        DamageData() {
            std::fill_n(r, 4, 1.0);
            std::fill_n(A, 5, -1.0);
            std::fill_n(warnings, 5, false);
        }
        double r[4];
        double A[5];
        bool warnings[5];
    };
    DamageData* damage_data;
 
    template <typename T>
    static double simpson_integration(double a, double b, int n, const T& f) {
        n *= 2;
        double r = (f(a) + f(b)) / 2;
        const double h = (b - a) / n;
        a += h;
        for (int i = 1; i != n; ++i, a += h) { r += (1 + i % 2) * f(a); }
        return 2 * r * h / 3;
    }
 
    template <typename D_M>
    class G_M {
    public:
        G_M(const D_M& d_M_, double E, double P, double e_stress)
            : d_M(d_M_), P2(P * P), s_E(e_stress * e_stress / (2 * E)), A(0) {}
        double operator()(double r) const {
            double tmp = (1 - P2) / (1 - (1 - d_M.value(r, A)) * P2);
            return tmp * tmp * s_E * d_M.d_value_d_r(r, A);
        }
        double value(double A_) {
            A = A_;
            const double ln_inv_K = std::log(1 / 20.0);
            const double s = -1 / A * ln_inv_K;
            return simpson_integration(1, 1 + s, 100, *this);
        }
 
    private:
        const D_M& d_M;
        double P2;
        double s_E;
        double A;
    };
 
    // Method that compute the adjustment parameter by secant method
    template <typename D_M>
    static double compute_adjustment_parameter(
          double E, double X, double Gf, double characteristic_length, G_M<D_M>& g_M, double tol) {
        X *= X;
        double lA_0;
        double lA_1;
        double lg_0;
        double lg_1;
        Gf = Gf / characteristic_length;
        {
            double A_1 = (2 * characteristic_length * X) / (2 * E * Gf - characteristic_length * X);
            if (A_1 <= 0) {
                Exception() << "The element size are too big for this material" << THROW;
            }
            double g_1 = g_M.value(A_1);
            double A_0 = 0.5 * A_1;
            double g_0 = g_M.value(A_0);
 
            for (int nbiter = 0; nbiter != 50 && g_1 > Gf; ++nbiter) {
                A_1 *= 10;
                g_1 = std::max(g_M.value(A_1), 1e-200);
            }
            lA_1 = std::log(A_1);
            lg_1 = std::log(g_1);
 
            for (int nbiter = 0; nbiter != 50 && g_0 < Gf; ++nbiter) {
                A_0 *= 0.1;
                g_0 = g_M.value(A_0);
            }
            lA_0 = std::log(A_0);
            lg_0 = std::log(g_0);
        }
        const double lGc = std::log(Gf);
        for (int nbiter = 0; nbiter != 20; ++nbiter) {
            double lA_2 = lA_1 - (lg_1 - lGc) * (lA_1 - lA_0) / (lg_1 - lg_0);
            lg_0 = lg_1;
            lA_2 = std::min(lA_2, 200.0);
            double A = std::exp(lA_2);
            double g_1;
            for (int nbiter1 = 0; nbiter1 != 10; ++nbiter1) {
                g_1 = g_M.value(A);
                if (g_1 == 0) {
                    lA_2 /= 2;
                    A = std::exp(lA_2);
                } else {
                    break;
                }
            }
            if (g_1 != g_1) { Exception() << THROW; }
            if (std::abs(g_1 - Gf) <= tol) { return A; }
            lg_1 = std::log(g_1);
            if (std::abs(lg_1 - lg_0) <= tol) { return A; }
            lA_0 = lA_1;
            lA_1 = lA_2;
        }
        Exception() << THROW;
        return 0;
    }
 
    class D1minus {
    public:
        D1minus(const double Apm_, const double A1plus, const double rplus)
            : Apm(Apm_), f(1 - Apm + Apm / rplus * std::exp(A1plus * (1 - rplus))) {}
        double value(const double r, const double A) const {
            return 1 - std::exp(A * (1 - r)) * f / r;
        }
        double d_value_d_r(const double r, const double A) const {
            return f * std::exp(A * (1 - r)) * (1 + r * A) / (r * r);
        }
 
    private:
        const double Apm;
        const double f;
    };
 
    class D2plus {
    public:
        D2plus(const double g_) : g(g_) {}
        double value(const double r, const double A) const {
            const double f2plus = (g - 1 + std::sqrt((1 - g) * (1 - g) + 4 * g * r * r)) / (2 * g);
            return 1 - std::exp(A * (1 - f2plus)) / f2plus;
        }
        double d_value_d_r(const double r, const double A) const {
            const double s = std::sqrt((1 - g) * (1 - g) + 4 * g * r * r);
            return 4 * (g * r * std::exp(A * (1 - s)) * (g * (2 + A) - A * (1 - s)))
                   / (s * (g - 1 + s));
        }
 
    private:
        double g;
    };
 
    class D2minus {
    public:
        double value(const double r, const double A) const { return 1 - std::exp(A * (1 - r)) / r; }
        double d_value_d_r(const double r, const double A) const {
            return std::exp(A * (1 - r)) * (1 / (r * r) + A);
        }
    };
};
}  // namespace b2000
 
#endif