Question about getting trial acceleration

[jiejingz]

Hello Everyone,

I was wondering if there is a build-in function that can get trail accelerations?

For instance, in elastic material, we use get strain
in viscous material, we use get strain rate
I'm trying to find if there is anything that can get the "rate of strain rate"

Appreciate your help!

[fmk]

is this from the interp[reter or the code .. the UniaxialMaterial class does have a getStrainRate method .. if not implementd in a material we can put it ine

[jiejingz]

Hello fmk, thanks for the reply.

This is from the code I presume.

What I am trying to to do is to basically create a damping material that provide damping force proportional to the acceleration.

In the source code for uniaxial viscous material, where the damping force is proportional to the velocity, we can use getStrainrate method as our trial velocity
but for the uniaxial material I am trying to create, I can't really do this because I don't have/or don't know of a method that can directly get me the trial acceleration?

Is there an efficient way to do this? Thank you！

[fmk]

you need to store the rate in your material .. only a few elements provide the strain rate to the material .. if you need it for a material that will be in an element that does not supply the rate you will have to compute it yourself using the strains and the global delta T variable

[jiejingz]

Hello fmk,

Here is what I have done so far, I tried to create a viscous material using the strains and the global delta T, but it seems ops_Dt is 0 as my displacements and velocity are either #IND or #QNAN.

Also, if I replace ops_Dt with a finite number (e.g. 0.05), it gives reasonable results, but it is not very accurate and reliable because damping force depend on that number.

here is my cpp code (i used it to build a dll):


#include <elementAPI.h>
#include <math.h>

#include "Material_1.h"
#include <Vector.h>
#include <Channel.h>


#include <OPS_Globals.h>

#ifdef _USRDLL
#define OPS_Export extern "C" _declspec(dllexport)
#elif _MACOSX
#define OPS_Export extern "C" __attribute__((visibility("default")))
#else
#define OPS_Export extern "C"
#endif

static int numMaterial_1 = 0;

OPS_Export void *
OPS_Material_1()
{
// print out some KUDO's
if (numMaterial_1 == 0) {
opserr << "Material_1 unaxial material - Written by fmk UC Berkeley Copyright 2008 - Use at your Own Peril\n";
numMaterial_1 = 1;
}

// Pointer to a uniaxial material that will be returned
UniaxialMaterial *theMaterial = 0;

//
// parse the input line for the material parameters
//

int iData[1];
double dData[3];
int numData;
numData = 1;
if (OPS_GetIntInput(&numData, iData) != 0) {
opserr << "WARNING invalid uniaxialMaterial Material_1 tag" << endln;
return 0;
}

numData = 3;
if (OPS_GetDoubleInput(&numData, dData) != 0) {
opserr << "WARNING invalid c,a,and minV";
return 0;
}

//
// create a new material
//

theMaterial = new Material_1(iData[0], dData[0], dData[1], dData[2]);

if (theMaterial == 0) {
opserr << "WARNING could not create uniaxialMaterial of type QuadraticPPcpp\n";
return 0;
}

// return the material
return theMaterial;
}

Material_1::Material_1(int tag, double c, double a, double minV)
:UniaxialMaterial(tag, MAT_TAG_Viscous),
trialRate(0.0), C(c), Alpha(a), minVel(minV)
{
if (Alpha < 0.0) {
opserr << "Material_1::Material_1 -- Alpha < 0.0, setting to 1.0\n";
Alpha = 1.0;
}

minVel = fabs(minVel);
if (minVel == 0.0) {
opserr << "Material_1::Material_1 -- minVel == 0.0, setting to 1.0e-21\n";
minVel = 1.0e-21;
}
}

Material_1::Material_1()
:UniaxialMaterial(0, MAT_TAG_Viscous),
trialRate(0.0), trialStrain(0.0), C(0.0), Alpha(0.0), minVel(1e-11),
trialStress(0.0), trialTangent(0.0), commitStrain(0.0), commitStress(0.0), commitTangent(0.0)
{

}

Material_1::~Material_1()
{
// does nothing
}

int
Material_1::setTrialStrain(double strain, double strainRate)
{
trialStrain = strain;
//trialRate = strainRate;
trialRate = (trialStrain-commitStrain)/ops_Dt;

//printf
return 0;
}

double
Material_1::getStress(void)
{
double stress = 0.0;
double absRate = fabs(trialRate);

if (absRate > minVel)
stress = C*pow(absRate, Alpha);
else
stress = C*pow(minVel, Alpha);

//stress = C*pow(absRate, Alpha);

if (trialRate < 0.0)
return -stress;
else
return stress;
}

double
Material_1::getTangent(void)
{
return 0.0;
}

double
Material_1::getInitialTangent(void)
{
return 0.0;
}

double
Material_1::getDampTangent(void)
{
double absRate = fabs(trialRate);

if (absRate < minVel)
return Alpha*C*pow(minVel, Alpha - 1.0);
else
return Alpha*C*pow(absRate, Alpha - 1.0);
}


double
Material_1::getStrain(void)
{
return 0.0;
}

double
Material_1::getStrainRate(void)
{
return trialRate;
}

int
Material_1::commitState(void)
{

commitStrain = trialStrain;

return 0;
}

int
Material_1::revertToLastCommit(void)
{
trialStrain = commitStrain;
return 0;
}

int
Material_1::revertToStart(void)
{
trialRate = 0.0;

return 0;
}

UniaxialMaterial *
Material_1::getCopy(void)
{
Material_1 *theCopy = new Material_1(this->getTag(), C, Alpha);

theCopy->trialRate = trialRate;

return theCopy;
}

int
Material_1::sendSelf(int cTag, Channel &theChannel)
{
int res = 0;
static Vector data(5);
data(0) = this->getTag();
data(1) = C;
data(2) = Alpha;
data(3) = trialRate;
data(4) = minVel;
res = theChannel.sendVector(this->getDbTag(), cTag, data);
if (res < 0)
opserr << "Material_1::sendSelf() - failed to send data\n";

return res;
}

int
Material_1::recvSelf(int cTag, Channel &theChannel,
FEM_ObjectBroker &theBroker)
{
int res = 0;
static Vector data(5);
res = theChannel.recvVector(this->getDbTag(), cTag, data);

if (res < 0) {
opserr << "Material_1::recvSelf() - failed to receive data\n";
C = 0;
this->setTag(0);
}
else {
this->setTag((int)data(0));
C = data(1);
Alpha = data(2);
trialRate = data(3);
minVel = data(4);
}

return res;
}

void
Material_1::Print(OPS_Stream &s, int flag)
{
s << "Viscous tag: " << this->getTag() << endln;
s << " C: " << C << endln;
s << " Alpha: " << Alpha << endln;
s << " minVel: " << minVel << endln;
}

[jiejingz]

I guess for now my primary concern is accessing ops_Dt with a dll

For some reason I believe ops_Dt is giving me 0 all the time. Is there a way around this?

[fmk]

the only way around it till i can put it in the source code is to keep a variable in the material to keep track of the last committed time. through the ops_TheActiveDomain pointer you can obtain the current time and the dT can be calculated if you have stored the committed time.

if that does not work let me know and i will build something.

[jiejingz]

Dear fmk,

Thanks for your reply. I think this method works well, but there is something I do not quite understand.

In implementing a mass damper material, here is my main code, where C is the mass of the mass damper, and I left getTangent, getInitialTangent, and getDampedTangeng methods to return 0;

_________________________________________________________________________________
Material_5::setTrialStrain(double strain, double strainRate)
{
double dt;
currentTime = ops_TheActiveDomain->getCurrentTime();
dt = currentTime - commitTime;
trialStrain = strain;
trialStrainRate = strainRate;

double nextStrain = 0;

//nextStrain = trialStrain + trialStrainRate*dt;
//Acceleration = (commitStrain - 2 * trialStrain + nextStrain)/pow(dt,2);

Acceleration = (trialStrainRate - commitStrainRate)/dt;

if (Acceleration >= 0)
{
trialStress = fabs(C*Acceleration);
}
else
{
trialStress = -fabs(C*Acceleration);
}

return 0;
}
____________________________________________________________________________________

However, unless I specify a value proportional to C for the damped tangent determined through trial and error, e.g.:
____________________________________________________________________________________
Material_5::getDampTangent(void)
{

return 39.9*C;
}
____________________________________________________________________________________

The single dof system with a mass, a spring (elastic zero length), and this mass damper (Material 5 zero length element) will produce sinusoidal increasing displacements and never reach steady state.

Why is this the case? What is the implication of the damped tangent? The coefficient is not at all intuitive in this case.

Many Thanks!!

[fmk]

the problem could be the calculation of acceleration .. did you look at the theory for the integrator you are using to see how the acceleration changes with regards to the disp and vel changes and try and make your computed accel agree with this approximation.

[jiejingz]

Hello Frank,

By aligning the acceleration calculation to newmark's method (also my integration method),the code worked pretty well. I just have two more questions regarding Opensees:

(1) For the above to work, I need to set GetDampedTangent method to return a value proportional to the inverse of dt, but couldn't figure out the exact reason, is this method related to numerical damping?

(2) Is there a pointer that allows me to access the inputs of the integrator (e.g. in tcl, I would state "integrator Newmark 0.5 0.25", is there a way to access the 0.5 and 0.25" in the c++ code?

Many thanks,

Jiejing Zhou

[fmk]

1) I would say it has more to do with the fact that the stress is related to accelerations. C/dt for the damp tangent would make sense. I think an approach without the 1/dt would be like creating a Viscous material in which stress is C*V, the tangent stiffness was C and the damping tangent was 0.

2) No.

[Gholamreza]

Hello jiejingz,
If you are done with your code, Could you share that with us ?