Contents     Functions              PDF Index

## fe_time,of_time

Purpose

Computation of time and non linear responses.

Syntax

```  def=fe_time(model)
def=fe_time(TimeOpt,model)
[def,model,opt]=fe_time(TimeOpt,model)
model=fe_time('TimeOpt...',model)
TimeOpt=fe_time('TimeOpt...')
```

Description

fe_time groups static non-linear and transient solvers to compute the response of a FE model given initial conditions, boundary conditions, load case (section 4.2.3) and time parameters. Note that you may access to the fe_time commands graphically with the simulate tab of the feplot GUI. See tutorial (section 4.5) on how to compute a response.

### Solvers and options

Three types of time integration algorithm are possible: the Newmark schemes, the Theta-method, and the time Discontinuous Galerkin method. Implicit and explicit methods are implemented for the Newmark scheme, depending on the Newmark coefficients β and γ, and non linear problems are supported.

The parameters of a simulation are stored in a time option data structure TimeOpt given as input argument or in a model.Stack entry info,TimeOpt. Initial conditions are stored as a curve,q0 entry.

The solvers selected by the string TimeOpt.Method are

• newmark, linear Newmark
• NLNewmark, non linear Newmark (with Newton iterations)
• staticNewton, static Newton
• Theta, Theta-Method (linear)
• Euler method for first order time integration.
• dg, Discontinuous Galerkin

Here is a simple example to illustrate the common use of this function.

``` model=fe_time('demo bar'); % build the model

% set the time options in model.Stack
model=fe_time('TimeOpt Newmark .25 .5 0 1e-4 100',model);

def=fe_time(model); % compute the response
```

fe_time can also be called with TimeOpt as the first argument. This is often more convenient when the user modifies options fields by hand

``` def=fe_time(TimeOpt,model);
```

The TimeOpt data structure has fields

 Method selection of the solver Opt if any, for example for Newmark[beta gamma t0 deltaT Nstep] OutInd DOF output indices (see 2D example). This selection is based on the state DOFs which can be found using fe_case(model,'GettDof'). MaxIter maximum number of iterations. nf optional value of the first residual norm. IterInit,IterEnd see the non linear solvers. Jacobian string to be evaluated to generate a factored jacobian matrix in matrix or ofact object ki. The default jacobian matrix is 'ki=ofact(model.K{3}+2/dt*model.K{2} +4/(dt*dt)*model.K{1});' for the dynamic case and'ki=ofact(model.K{3});' for the static case. JacobianUpdate controls the update of Jacobian in Newton and quasi-Newton loops. Use 1 for updates and 0 for a fixed Jacobian (default). Residual Callback evaluated for residual. The default residual is method dependent. InitAcceleration optional field to be evaluated to initialize the acceleration field. OutputFcn string to be evaluated for post-processing or time vector containing the output time steps FinalCleanupFcn string to be evaluated for final post-processing of the output time steps
 c_u, c_v, c_a optional observation matrices for displacement, velocity and acceleration outputs. lab_u, lab_v, lab_a optional cell array containing labels describing each output (lines of observation matrices) NeedUVA [NeedU NeedV NeedA], if NeedU is equal to 1, output displacement, etc. OutputInit optional string to be evaluated to initialize the output (before the time loop) SaveTimes optional time vector, saves time steps on disk TimeVector optional value of computed time steps, if exists TimeVector is used instead of time parameters RelTol threshold for convergence tests. The default is the OpenFEM preferencegetpref('OpenFEM','THRESHOLD',1e-6);

### Input and output options

This section details the applicable input and the output options.

Initial conditions may be provided in a model.Stack entry of type info named q0 or in an input argument q0. q0 is a data structure containing def and DOF fields as in a FEM result data structure (section 4.5). If any, the second column gives the initial velocity. If q0 is empty, zero initial conditions are taken. In this example, a first simulation is used to determine the initial conditions of the final simulation.

``` model=fe_time('demo bar');
TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 100');
TimeOpt.NeedUVA=[1 1 0];
% first computation to determine initital conditions
def=fe_time(TimeOpt,model);

% no input force

% Setting initial conditions
q0=struct('def',[def.def(:,end) def.v(:,end)],'DOF',def.DOF);
model=stack_set(model,'curve','q0',q0);

def=fe_time(TimeOpt,model);
```

An alternative call is possible using input arguments

``` def=fe_time(TimeOpt,model,Case,q0)
```

In this case, it is the input argument q0 which is used instead of an eventual stack entry.

You may define the time dependence of a load using curves as illustrated in section 7.9.

You may specify the time steps by giving the 'TimeVector'

``` TimeOpt=struct('Method','Newmark','Opt',[.25 .5 ],...
'TimeVector',linspace(0,100e-4,101));
```

This is useful if you want to use non constant time steps. There is no current implementation for self adaptive time steps.

To illustrate the output options, we use the example of a 2D propagation. Note that this example also features a time dependent DofLoad excitation (see fe_case) defined by a curve, (see fe_curve), here named Point load 1.

``` model=fe_time('demo 2d');
TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 50');
```

You may specify specific output by selecting DOF indices as below

``` i1=fe_case(model,'GettDof'); i2=feutil('findnode y==0',model)+.02;
TimeOpt.OutInd=fe_c(i1,i2,'ind');
model=stack_set(model,'info','TimeOpt',TimeOpt);
def=fe_time(model); % Don't animate this (only bottom line)
```

You may select specific output time step using TimeOpt.OutputFcn as a vector

``` TimeOpt.OutputFcn=[11e-4 12e-4];
TimeOpt=feutil('rmfield',TimeOpt','OutInd');
model=stack_set(model,'info','TimeOpt',TimeOpt);
def=fe_time(model); % only two time steps saved
```

or as a string to evaluate. The output is the out local variable in the fe_time function and the current step is j1+1. In this example the default output function (for TimeOpt.NeedUVA=[1 1 1]) is used but specified for illustration

``` TimeOpt.OutputFcn=['out.def(:,j1+1)=u;' ...
'out.v(:,j1+1)=v;out.a(:,j1+1)=a;'];
model=stack_set(model,'info','TimeOpt',TimeOpt);
def=fe_time(model); % full deformation saved
```

This example illustrates how to display the result (see feplot) and make a movie

```  cf=feplot(model,def);
fecom('ColorDataEvalA');
fecom(cf,'SetProp sel(1).fsProp','FaceAlpha',1,'EdgeAlpha',0.1);
cf.ua.clim=[0 2e-6];fecom(';view2;AnimTime;ch20;scd1e-2;');
st=fullfile(getpref('SDT','tempdir'),'test.avi');
fecom(['animavi ' st]);fprintf('\nGenerated movie %s\n',st);
```

Note that you must choose the Anim:Time option in the feplot GUI.

You may want to select outputs using observations matrix

``` model=fe_time('demo bar'); Case=fe_case('gett',model);
i1=feutil('findnode x>30',model);
TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 100');
TimeOpt.c_u=fe_c(Case.DOF,i1+.01);          % observation matrix
TimeOpt.lab_u=fe_c(Case.DOF,i1+.01,'dofs'); % labels

def=fe_time(TimeOpt,model);
```

If you want to specialize the output time and function you can specify the SaveTimes as a time vector indicating at which time the SaveFcn string will be evaluated. A typical TimeOpt would contain

``` TimeOpt.SaveTimes=[0:Ts:TotalTime];
TimeOpt.SaveFcn='My_function(''Output'',u,v,a,opt,out,j1,t);';
```

### Cleanup

The field FinalCleanupFcn of the TimeOpt can be used to specify what is done just after the time integration.
fe_simul provides a generic clean up function which can be called using
opt.FinalCleanupFcn='fe_simul(''fe_timeCleanUp'',model)';
Some command options can be used.
-cf i stores the result of time integration in the stack of iiplot or feplot figure number i. i=-1 can be specified to use current iiplot figure and i=-2 for current feplot figure. Displacements are stored in curve,disp entry of the stack. Velocities and accelerations (if any) are respectively stored in the curve,vel and curve,acc stack entries. If command option -reset is present, existent stack entries (disp, vel, acc, etc. ...) are lost whereas if not stack entries name are incremented (disp(1), disp(2), etc. ...).
Command option '-ExitFcn"AnotherCleanUpFcn"' can be used to call an other clean up function just after fe_simul('fe_timeCleanUp') is performed.

### newmark

For the Newmark scheme, TimeOpt has the form

``` TimeOpt=struct('Method','Newmark','Opt',Opt)
```

where TimeOpt.Opt is defined by

```   [beta gamma t0 deltaT Nstep]
```

beta and gamma are the standard Newmark parameters [36] ([0 0.5] for explicit and default at [.25 .5] for implicit), t0 the initial time, deltaT the fixed time step, Nstep the number of steps.

The default residual is r = (ft(j1,:)*fc'-v'*c-u'*k)'; (notice the sign change when compared to NLNewmark).
This is a simple 1D example plotting the propagation of the velocity field using a Newmark implicit algorithm. Rayleigh damping is declared using the info,Rayleigh case entry.

``` model=fe_time('demo bar');
data=struct('DOF',2.01,'def',1e6,...
'curve',fe_curve('test ricker 10e-4 1'));
TimeOpt=struct('Method','Newmark','Opt',[.25 .5 3e-4 1e-4 100],...
'NeedUVA',[1 1 0]);
def=fe_time(TimeOpt,model);

% plotting velocity (propagation of the signal)
def_v=def;def_v.def=def_v.v; def_v.DOF=def.DOF+.01;
feplot(model,def_v);
if sp_util('issdt'); fecom(';view2;animtime;ch30;scd3');
else; fecom(';view2;scaledef3'); end
```

### dg

The time discontinuous Galerkin is a very accurate time solver [49] [50] but it is much more time consuming than the Newmark schemes. No damping and no non linearities are supported for Discontinuous Galerkin method.

The options are [unused unused t0 deltaT Nstep Nf], deltaT is the fixed time step, Nstep the number of steps and Nf the optional number of time step of the input force.

This is the same 1D example but using the Discontinuous Galerkin method:

``` model=fe_time('demo bar');
TimeOpt=fe_time('TimeOpt DG Inf Inf 0. 1e-4 100');
TimeOpt.NeedUVA=[1 1 0];
def=fe_time(TimeOpt,model);

def_v=def;def_v.def=def_v.v; def_v.DOF=def.DOF+.01;
feplot(model,def_v);
if sp_util('issdt'); fecom(';view2;animtime;ch30;scd3'); ...
else; fecom(';view2;scaledef3'); end
```

### NLNewmark

For the non linear Newmark scheme, TimeOpt has the same form as for the linear scheme (method Newmark). Additional fields can be specified in the TimeOpt data structure

 Jacobian string to be evaluated to generate a factored jacobian matrix in matrix or ofact object ki. The default jacobian matrix is 'ki=ofact(model.K{3}+2/dt*model.K{2}' +4/(dt*dt)*model.K{1});' Residual Defines the residual used for the Newton iterations of each type step. It is typically a call to an external function. The default residual is 'r = model.K{1}*a+model.K{2}*v+model.K{3}*u-fc;' where fc is the current external load, obtained using (ft(j1,:)*fc')' at each time step. IterInit evaluated when entering the correction iterations. This can be used to initialize tolerances, change mode in a co-simulation scheme, etc. IterEnd evaluated when exiting the correction iterations. This can be used to save specific data, ... IterFcn Correction iteration algorithm function, available are iterNewton (default when ommitted) or iterNewton_Sec. Details of the implementation are given in the staticNewton below. MaxIterSec for iterNewton_Sec applications (see staticNewton). ResSec for iterNewton_Sec applications (see staticNewton).

### staticNewton

For non linear static problems, the Newton solver iterNewton is used. TimeOpt has a similar form as with the NLNewmark method but no parameter Opt is used.

An increment control algorithm iterNewton_Sec can be used when convergence is difficult or slow (as it happens for systems showing high stiffness variations). The Newton increment Δ q is then the first step of a line search algorithm to optimize the corective displacement increment ρ Δ q, ρ∈ R in the iteration. This optimum is found using the secant iteration method. Only a few optimization iterations are needed since this does not control the mechanical equilibrium but only the relevace of the Newton increment. Each secant iteration requires two residual computations, which can be costly, but more efficient when a large number of standard iterations (matrix inversion) is required to obtain convergence.

Fields can be specified in the TimeOpt data structure

 Jacobian string to be evaluated to generate a factored jacobian matrix in matrix or ofact object ki. The default jacobian matrix is 'ki=ofact(model.K{3});' Residual Defines the residual used for the Newton iterations of each type step. It is typically a call to an external function. The default residual is 'r = model.K{3}*u-fc;' IterInit evaluated when entering the Newton solver. This can be used to initialize tolerances, change mode in a co-simulation scheme, etc. IterEnd evaluated when exiting the Newton solver. This can be used to save specific data, ... IterFcn Correction iteration function. It is a string naming an fe_timeiteration subfunction available, or an external function. When performing the time simulation initialization, the string will be replaced by the function handle (e.g. @iterNewton). The available iteration algorithms in fe_timeare either iterNewton, which is defaulted when the IterFcn field is ommitted and performs basic Newton iterations, either iterNewton_Sec which is the Newton increment control algorithm. MaxIterSec Maximum secant iterations for the iterNewton_Sec iteration algorithm. The default is 3 when ommitted. ResSec Residual evaluation for the secant iterations of the iterNewton_Sec iteration algorithm. When ommited, fe_timetries to interpret the Residual field. The function must fill in the secant residual evaluation r1 which two columns will contain the residual for solution rho(1)*dq and rho(2)*dq. The default ResSec field will be then 'r1(:,1) = model.K{3}*(u-rho(1)*dq)-fc; r1(:,2) = model.K{3}*(u-rho(2)*dq)-fc;'.

Below is a trivial demonstration of formats and call for non-linear bar statically loaded in 5 steps.

```mdl=femesh('testbar1');

mdl=fe_case(mdl,'FixDof','left',1);
mdl=fe_case(mdl,'setcurve','right', ...
fe_curve('testramp 5 1')); % 5 steps gradual load
mdl=stack_set(mdl,'info','TimeOpt', ...
struct('Opt',[],'Method','staticnewton',...
'Jacobian','ki=basic_jacobian(model,ki,0.,0.,opt.Opt);',...
'Residual','r = model.K{3}*u-fc-model.K{3}*u.^1.5;'));

def=fe_time(mdl);feplot(mdl,def);fecom('showdefarrow');
```

### numerical damping for Newmark, HHT-alpha schemes

You may want to use numerical damping in a time integration scheme, the first possibility is to tune the Newmark parameters using a coefficient α such that β = (1+α)2/4 and γ=1/2+α. This is known to implement too much damping at low frequencies and is very depending on the time step [36].

A better way to implement numerical damping is to use the HHT-α method which applies the Newmark time integration scheme to a modified residual balancing the forces with the previous time step.

For the HHT-α scheme, TimeOpt has the form

``` TimeOpt=struct('Method','nlnewmark','Opt',Opt,...
'HHTalpha',alpha)
```

where TimeOpt.Opt is defined by

```   [beta gamma t0 deltaT Nstep]
```

beta and gamma are the standard Newmark parameters  [36] with numerical damping, t0 the initial time, deltaT the fixed time step, Nstep the number of steps.

The automatic TimeOpt generation call takes the form [alpha unused t0 deltaT Nstep] and will compute the corresponding β, γ parameters.

This is a simple 1D example plotting the propagation of the velocity field using the HHT-α implicit algorithm:

``` model=fe_time('demo bar');
TimeOpt=fe_time('TimeOpt hht .05 Inf 3e-4 1e-4 100');
TimeOpt.NeedUVA=[1 1 0];
def=fe_time(TimeOpt,model);
```

The call

``` TimeOpt=fe_time('TimeOpt hht .05 Inf 3e-4 1e-4 100');
```

is strictly equivalent to

```TimeOpt=struct('Method','nlnewmark',...
'Opt',[.275625 .55 3e-4 1e-4 100],...
'HHTalpha',.05);
```

### Theta

The θ-method is a velocity based solver, whose formulation is given for example in [51, 52]. It considers the acceleration as a distribution, thus relaxing discontinuity problems in non-smooth dynamics. Only a linear implementation is provided in fe_time. The user is nevertheless free to implement a non-linear iteration, through his own IterFcn.

This method takes only one integration parameter for its scheme, θ set by default at 0.5. Any values between 0.5 and 1 can be used, but numerical damping occurs for θ>0.5.

The TimeOpt.Opt takes the form [theta unused t0 deltaT Nstep].

This is a simple 1D example plotting the propagation of the velocity field using the θ-Method:

```model=fe_time('demo bar');
TimeOpt=fe_time('TimeOpt theta .5 0 3e-4 100');
def=fe_time(TimeOpt,model);
```

### Euler

This method can be used to integrate first order problem of the form Mq + Kq=F. One can use it to solve transient heat diffusion equation (see p_heat).

Integration scheme is of the form qn+1=qn+(1−θ)h qnh qn+1
θ can be define in opt.Opt(1). Explicit Euler (θ=0) is not implemented at this time. Best accuracy is obtained with θ=1/2 (Crank-Nicolson).

### of_time

The of_time function is a low level function dealing with CPU and/or memory consuming steps of a time integration.

The commands are

 lininterp linear interpolation. storelaststep pre-allocated saving of a time step. interp Time scheme interpolations (low level call).

The lininterp command which syntax is

out = of_time ('lininterp',table,val,last) ,

computes val containing the interpolated values given an input table which first column contains the abscissa and the following the values of each function. Due to performance requirements, the abscissa must be in ascending order. The variable last contains [i1 xi si], the starting index (beginning at 0), the first abscisse and coordinate. The following example shows the example of 2 curves to interpolate:

```out=of_time('lininterp',[0 0 1;1 1 2;2 2 4],linspace(0,2,10)',[0 0 0])
```

Warning : this command modifies the variable last within a given function this may modify other identical constants in the same m-file. To avoid any problems, this variable should be generated using zeros (the Matlab function) to assure its memory allocation independence.

The storelaststep command makes a deep copy of the displacement, velocity and acceleration fields (stored in each column of the variable uva in the preallocated variables u, v and a following the syntax:

of_time('storelaststep',uva,u,v,a);

The interp command is used by fe_time when the user gives a TimeVector in the command. For time t1 and the uva matrix containing in each column, displacement, velocity and acceleration at the preceding time step t0, the interpolation strategy depends on the arguments given to of_time.

Warning : this command modifies out.def at very low level, out.def thus cannot be initialized by simple numerical values, but by a non trivial command (use zeros(1) instead of 0 for example) to ensure the unicity of this data in memory.

For a Newmark or HHT-alpha scheme, the low level call command is

```of_time ('interp', out, beta,gamma,uva,a, t0,t1,model.FNL);
```

where beta and gamma are the coefficients of the Newmark scheme, first two values of opt.Opt.

Thus the displacement (u1) and velocity (v1) at time t1 will be computed from the displacement (u0), velocity (v0), acceleration (a0) stored in uva, the new acceleration a (a1), and the time step (h=t1−t0) as

v1 = v0 + h (1− γ) a0 + h γ a1
u1 = u0 + h v0 + h2 (
 1 2
− β) a0 + h2 β a1
(9.4)

For the Theta-Method scheme, the low level command is

```of_time ('interp', out, opt.Opt(1),[],uva,v, t0,t1,model.FNL);
```

Thus the displacement (u1) at time t1 will be computed from the displacement (u0), velocity (v0), stored in uva, the new velocity v (v1), and the time step (h=t1−t0) as

 u1 = u0 + h (1−θ) v0 + h θ v1     (9.5)

For the staticnewton method, it is possible to use the same storage strategy (since it is optimized for performance), using

```of_time ('interp', out, [],[], [],u, t0,t1,model.FNL);
```

In this case no interpolation is performed.

Please note that this low-level call uses the internal variables of fe_time at the state where is is evaluated. It is then usefull to know that inside fe_time:

• at time tc=t(j1+1) using time step dt, values are t0=tc-dt and t1=tc.
• uva is generally stored in Case.uva.
• the current acceleration, velocity or displacement values when interpolation is performed are always a, v, and u.
• The out data structure must be preallocated and is modified by low level C calls. Expected fields are

 def must be preallocated with size length(OutInd) x length(data) data column vector of output times OutInd int32 vector of output indices, must be given cur [Step dt], must be given FNL possibly preallocated data structure to store non-linear loads. FNL.def must be length(model.FNL) by size(out.data,1) (or possibly size(out.FNL.data,1), in this case fieldnames must be def,DOF,data,cur)
• non linear loads in model.FNL are never interpolated.