Contents  
Functions  
Index
 
 PDF |
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.4) on how to compute a response.
Two types of time integration algorithm are possible : the Newmark schemes 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 info,q0 entry.
The solvers selected by the string TimeOpt.Method are
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 alpha t0 deltaT Nstep Nf] |
| 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 forst 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 | (only for newton) : default is 0, 0 if modified Newton (no update in Newton iterations), 1 if update in Newton iteration |
| 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 |
| 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 ouput (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 preference getpref('OpenFEM','THRESHOLD',1e-6); |
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 ans DOF fields as in a FEM result data structure (section 4.4). 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 determines 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
model=fe_case(model,'remove','Point load 1');
% Setting initial conditions
q0=struct('def',[def.def(:,end) def.v(:,end)],'DOF',def.DOF);
model=stack_set(model,'info','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 a eventual stack entry.
You may define the time dependance of a load using curves as illustrated in section 7.9.
You may a specify the time steps by giving the 'TimeVector'
TimeOpt=struct('Method','Newmark','Opt',[.25 .5 ],...
'TimeVector',linspace(0,100e-4,101));
This is usefull if you want to use non constant time steps. There is no current implementation for self adaptive time steps.
To illustrate the ouput options, we use the example of a 2D propagation
model=fe_time('demo 2d');
TimeOpt=fe_time('TimeOpt Newmark .25 .5 0 1e-4 50 10');
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);
You may select specific output time step using TimeOpt.OutputFcn as a vector
TimeOpt.OutputFcn=[11e-4 12e-4]; model=stack_set(model,'info','TimeOpt',TimeOpt); def=fe_time(model);
or as a string to evaluate. The ouput is the out local variable in the fe_timefunction 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);
This example illustrates how to display the result (see feplot) and make a movie
cf=feplot(model,def);
fecom('colordataa');
cf.ua.clim=[0 2e-6];fecom(';view2;animtime;ch20;scd1e-2;');
st=fullfile(getpref('SDT','tempdir'),'test.avi');
fecom(['animavi ' st]);
endNote that you must choose the Anim:Time option in the feplotGUI.
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 SaveTimesas 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);';
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 Nf]
beta and gamma are the standard Newmark parameters [33], t0 the initial time, deltaT the fixed time step, Nstep the number of steps and Nf the optional number of time step of the input force.
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 100 1 100e-4'));
model = fe_case(model,'DOFLoad','Point load 1',data);
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
The time discontinuous Galerkin is a very accurate time solver [46] [47] 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
For the non linear Newmark scheme, TimeOpt has the same form as for the linear scheme (method tt 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 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, ... |
For non linear static problems, the Newton solver is used. TimeOpt has a similar form as with the tt NLNewmark method but no parameter Opt is used. 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, ... |
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,'DofLoad','right',struct('DOF',2.01,'def',1));
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);try;fecom('showdefarrow');end
For the α-HHT scheme, TimeOpt has the form
TimeOpt=struct('Method','hht','Opt',Opt)
where TimeOpt.Opt is defined by
[beta gamma t0 deltaT Nstep Nf]
beta and gamma are the standard Newmark parameters [33], t0 the initial time, deltaT the fixed time step, Nstep the number of steps and Nf the optional number of time step of the input force.
This is a simple 1D example plotting the propagation of the velocity field using a Newmark implicit algorithm :
model=fe_time('demo bar');
TimeOpt=struct('Method','hht','Opt',[.25 .5 3e-4 1e-4 100]);
TimeOpt.NeedUVA=[1 1 0];
def=fe_time(TimeOpt,model);
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 |
| newmarkinterp | Newmark interpolation (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])
The storelaststep command makes a deep copy of the displacement, celerity and acceleration fields (stored in each column of the variable uva in the preallocated variables u, v and a following the syntax:
Warning : this command modifies the variable last within a given function this may modify other identical constants in the same m-file.
of_time('storelaststep',uva,u,v,a);
The newmarkinterp command is used by fe_time when the user gives a TimeVector in the command using a Newmark scheme. Given an acceleration vector a1 at time t1 and the uva matrix containing in each column, displacement, celerity and acceleration at the preceding time step t0, it interpolates according to Newmark scheme (see Geradin p.371 eq. 7.3.9) the displacement at time t1. The low level call of newmarkinterp is
of_time ('newmarkinterp', out, beta,gamma,uva,a1, t0,t1)
Inside fe_time, at time tc=t(j1+1) using timestep dt, values are t0=tc-dt and t1=tc.
The out data structure must be preallocated and is a modified input containing the following fields :
| OutInd | Output indice, must be given |
| cur | [Step dt], must be given |
| def | must be preallocated |
See also
©1991-2008 by SDTools