I am trying to implement hydraulic ram end stops/limit stops, this being one component within a wider hydraulic model. Each component contributes residual equations such as force equilibrium, mass continuity, etc. These equations are then assembled as a system of non-linear equations and solved using Newtons method with backtracking. In addition, some of the residual equations define an implicit first-order time integration in order for the model to propagate through time.
Essentially, I am looking for a clean, efficient, and robust method of limiting the ram stroke to stay between a lower and upper bound, and for a resulting pressure increase to occur if pushed up against one of the bounds. I am not concerned about conservation of energy/impact.
For context, I have written a simple representation of a single-acting hydraulic ram (i.e. only one fluid port) as follows:
Residual 0 - force equilibrium:
$R_0 = F_{ram} - X_{ramPressure}^t * AreaPiston = 0$
Residual 1 - ram stroke first-order time integration:
$R_1 = X_{stroke}^t - X_{stroke}^{t-1} - dt * X_{ramFlowrate}^t / AreaPiston = 0$
where $X^t_{(*)}$ represents a state coefficient at time t to be determined by the solver, all other parameters are pre-defined.
Things I have tried so far:
1. Penalty method
Inspired by contact problems in finite-element analysis, I have tried adding an end-stop stiffness multiplied by the end stop penetration to the force equilibrium as follows:
$R_0 = F_{ram} - X_{ramPressure}^t * AreaPiston + K_{endStop} * penetration = 0$
My problem with this approach has been similar to that quoted for FE contact problems: I require such a large stiffness for the model to accurately respect the end stop that the Jacobian becomes ill-conditioned and impacts the solver performance. Note, with this I also need to implement a piecewise smoothing function that transitions from end stop not-active to active as penetration approaches zero, I have fitted a 3rd order polynomial for this purpose.
2. Constraint on ram stroke
I have tried adding an additional equation that requires that the ram stroke does not exceed the end stop. Such an equation can be represented simply as follows:
$R_2 = X^t_{stroke} - minimum(X^t_{stroke}, endStop) = 0$
Being that I add another equation, I also add another state which I define to be a force on the ram to be determined by the solver, thus the force equilibrium becomes:
$R_0 = F_{ram} - X_{ramPressure}^t * AreaPiston + X^t_{endstopForce} = 0$
This approach can work well with a smooth implementation of $R_2$. However, the problem I am having is that when the ram is not at an end stop, $R2 = 0$, and thus the Jacobian of the system becomes singular. Whilst I can use a matrix factorisation such as QR with pivoting to handle this, my requirements prefer the faster LU factorisation which requires a non-singular matrix.
Thoughts:
- I know I can use state scaling to avoid having such large coefficients in my Jacobian. However, I think I would need to dynamically change the scaling factors depending on if the end stop is active or not which may be tricky to integrate with the solver.
- $R_2$ is essentially a piecewise function in which it makes the Jacobian singular in one region, and non-singular in another region. I'm wondering are there solvers/methods to avoid piecewise functions causing a system of equations to go singular?
- Can this problem be considered as a constraint on the system of equations or bound on a state? I am reluctant to move from solving a system of equations to formulating/solving as an optimisation as I worry it won't be as efficient.
- Is this a situation in which a Lagrange or augmented Lagrange method might help me?
- Could semi-smooth newton solvers avoid the need for adding piecewise smoothing functions between the end stop being active/not active?
Thank you for taking the time to read, I look forward to your responses.